320 lines
9.0 KiB
C++
320 lines
9.0 KiB
C++
//
|
|
// Copyright (c) 2010 Athanasios Iliopoulos
|
|
//
|
|
// Distributed under the Boost Software License, Version 1.0. (See
|
|
// accompanying file LICENSE_1_0.txt or copy at
|
|
// http://www.boost.org/LICENSE_1_0.txt)
|
|
//
|
|
|
|
#include <boost/numeric/ublas/assignment.hpp>
|
|
#include <boost/numeric/ublas/vector.hpp>
|
|
#include <boost/numeric/ublas/vector_proxy.hpp>
|
|
#include <boost/numeric/ublas/matrix_proxy.hpp>
|
|
#include <boost/numeric/ublas/vector_sparse.hpp>
|
|
#include <boost/numeric/ublas/matrix_sparse.hpp>
|
|
#include <boost/numeric/ublas/io.hpp>
|
|
#include <boost/numeric/ublas/matrix.hpp>
|
|
|
|
using namespace boost::numeric::ublas;
|
|
|
|
int main() {
|
|
// Simple vector fill
|
|
vector<double> a(3);
|
|
a <<= 0, 1, 2;
|
|
std::cout << a << std::endl;
|
|
// [ 0 1 2]
|
|
|
|
// Vector from vector
|
|
vector<double> b(7);
|
|
b <<= a, 10, a;
|
|
std::cout << b << std::endl;
|
|
// [ 0 1 2 10 0 1 2]
|
|
|
|
// Simple matrix fill
|
|
matrix<double> A(3,3);
|
|
A <<= 0, 1, 2,
|
|
3, 4, 5,
|
|
6, 7, 8;
|
|
std::cout << A << std::endl;
|
|
// [ 0 1 2 ]
|
|
// [ 3 4 5 ]
|
|
// [ 6 7 8 ]
|
|
|
|
// Matrix from vector
|
|
A <<= 0, 1, 2,
|
|
3, 4, 5,
|
|
a;
|
|
std::cout << A << std::endl;
|
|
// [ 0 1 2 ]
|
|
// [ 3 4 5 ]
|
|
// [ 0 1 2 ]
|
|
|
|
// Matrix from vector - column assignment
|
|
A <<= move(0,2), traverse_policy::by_column(),
|
|
a;
|
|
std::cout << A << std::endl;
|
|
// [ 0 1 0 ]
|
|
// [ 3 4 1 ]
|
|
// [ 0 1 2 ]
|
|
|
|
// Another matrix from vector example (watch the wraping);
|
|
vector<double> c(9); c <<= 1, 2, 3, 4, 5, 6, 7, 8, 9;
|
|
A <<= c;
|
|
std::cout << A << std::endl;
|
|
// [ 1 2 3 ]
|
|
// [ 4 5 6 ]
|
|
// [ 7 8 9 ]
|
|
|
|
// If for performance(Benchmarks are not definite about that) or consistency reasons you need to disable wraping:
|
|
static next_row_manip endr; //This can be defined globally
|
|
A <<= traverse_policy::by_row_no_wrap(),
|
|
1, 2, 3, endr,
|
|
4, 5, 6, endr,
|
|
7, 8, 9, endr;
|
|
// [ 1 2 3 ]
|
|
// [ 4 5 6 ]
|
|
// [ 7 8 9 ]
|
|
// If by default you need to disable wraping define
|
|
// BOOST_UBLAS_DEFAULT_NO_WRAP_POLICY, in the compilation options,
|
|
// so that you avoid typing the "traverse_policy::by_row_no_wrap()".
|
|
|
|
// Plus and minus assign:
|
|
A <<= fill_policy::index_plus_assign(),
|
|
3,2,1;
|
|
std::cout << A << std::endl;
|
|
// [ 4 4 4 ]
|
|
// [ 4 5 6 ]
|
|
// [ 7 8 9 ]
|
|
|
|
// Matrix from proxy
|
|
A <<= 0, 1, 2,
|
|
project(b, range(3,6)),
|
|
a;
|
|
std::cout << A << std::endl;
|
|
// [ 0 1 2 ]
|
|
// [10 0 1 ]
|
|
// [ 6 7 8 ]
|
|
|
|
// Matrix from matrix
|
|
matrix<double> B(6,6);
|
|
B <<= A, A,
|
|
A, A;
|
|
std::cout << B << std::endl;
|
|
// [ A A ]
|
|
// [ A A ]
|
|
|
|
// Matrix range (vector is similar)
|
|
B = zero_matrix<double>(6,6);
|
|
matrix_range<matrix<double> > mrB (B, range (1, 4), range (1, 4));
|
|
mrB <<= 1,2,3,4,5,6,7,8,9;
|
|
std::cout << B << std::endl;
|
|
// [ 0 0 0 0 0 0]
|
|
// [ 0 1 2 3 0 0]
|
|
// [ 0 4 5 6 0 0]
|
|
// [ 0 0 0 0 0 0]
|
|
// [ 0 0 0 0 0 0]
|
|
// [ 0 0 0 0 0 0]
|
|
|
|
// Horizontal concatenation can be achieved using this trick:
|
|
matrix<double> BH(3,9);
|
|
BH <<= A, A, A;
|
|
std::cout << BH << std::endl;
|
|
// [ A A A]
|
|
|
|
// Vertical concatenation can be achieved using this trick:
|
|
matrix<double> BV(9,3);
|
|
BV <<= A,
|
|
A,
|
|
A;
|
|
std::cout << BV << std::endl;
|
|
// [ A ]
|
|
// [ A ]
|
|
// [ A ]
|
|
|
|
// Watch the difference when assigning matrices for different traverse policies:
|
|
matrix<double> BR(9,9, 0);
|
|
BR <<= traverse_policy::by_row(), // This is the default, so this might as well be omitted.
|
|
A, A, A;
|
|
std::cout << BR << std::endl;
|
|
// [ A A A]
|
|
// [ 0 0 0]
|
|
// [ 0 0 0]
|
|
|
|
matrix<double> BC(9,9, 0);
|
|
BC <<= traverse_policy::by_column(),
|
|
A, A, A;
|
|
std::cout << BC << std::endl;
|
|
// [ A 0 0]
|
|
// [ A 0 0]
|
|
// [ A 0 0]
|
|
|
|
// The following will throw a run-time exception in debug mode (matrix mid-assignment wrap is not allowed) :
|
|
// matrix<double> C(7,7);
|
|
// C <<= A, A, A;
|
|
|
|
// Matrix from matrix with index manipulators
|
|
matrix<double> C(6,6,0);
|
|
C <<= A, move(3,0), A;
|
|
// [ A 0 ]
|
|
// [ 0 A ]
|
|
|
|
// A faster way for to construct this dense matrix.
|
|
matrix<double> D(6,6);
|
|
D <<= A, zero_matrix<double>(3,3),
|
|
zero_matrix<double>(3,3), A;
|
|
// [ A 0 ]
|
|
// [ 0 A ]
|
|
|
|
// The next_row and next_column index manipulators:
|
|
// note: next_row and next_column functions return
|
|
// a next_row_manip and and next_column_manip object.
|
|
// This is the manipulator we used earlier when we disabled
|
|
// wrapping.
|
|
matrix<double> E(2,4,0);
|
|
E <<= 1, 2, next_row(),
|
|
3, 4, next_column(),5;
|
|
std::cout << E << std::endl;
|
|
// [ 1 2 0 5 ]
|
|
// [ 3 4 0 0 ]
|
|
|
|
// The begin1 (moves to the begining of the column) index manipulator, begin2 does the same for the row:
|
|
matrix<double> F(2,4,0);
|
|
F <<= 1, 2, next_row(),
|
|
3, 4, begin1(),5;
|
|
std::cout << F << std::endl;
|
|
// [ 1 2 5 0 ]
|
|
// [ 3 4 0 0 ]
|
|
|
|
// The move (relative) and move_to(absolute) index manipulators (probably the most useful manipulators):
|
|
matrix<double> G(2,4,0);
|
|
G <<= 1, 2, move(0,1), 3,
|
|
move_to(1,3), 4;
|
|
std::cout << G << std::endl;
|
|
// [ 1 2 0 3 ]
|
|
// [ 0 0 0 4 ]
|
|
|
|
// Static equivallents (faster) when sizes are known at compile time:
|
|
matrix<double> Gs(2,4,0);
|
|
Gs <<= 1, 2, move<0,1>(), 3,
|
|
move_to<1,3>(), 4;
|
|
std::cout << Gs << std::endl;
|
|
// [ 1 2 0 3 ]
|
|
// [ 0 0 0 4 ]
|
|
|
|
// Choice of traverse policy (default is "row by row" traverse):
|
|
|
|
matrix<double> H(2,4,0);
|
|
H <<= 1, 2, 3, 4,
|
|
5, 6, 7, 8;
|
|
std::cout << H << std::endl;
|
|
// [ 1 2 3 4 ]
|
|
// [ 5 6 7 8 ]
|
|
|
|
H <<= traverse_policy::by_column(),
|
|
1, 2, 3, 4,
|
|
5, 6, 7, 8;
|
|
std::cout << H << std::endl;
|
|
// [ 1 3 5 7 ]
|
|
// [ 2 4 6 8 ]
|
|
|
|
// traverse policy can be changed mid assignment if desired.
|
|
matrix<double> H1(4,4,0);
|
|
H1 <<= 1, 2, 3, traverse_policy::by_column(), 1, 2, 3;
|
|
|
|
std::cout << H << std::endl;
|
|
// [1 2 3 1]
|
|
// [0 0 0 2]
|
|
// [0 0 0 3]
|
|
// [0 0 0 0]
|
|
|
|
// note: fill_policy and traverse_policy are namespaces, so you can use them
|
|
// by a using statement.
|
|
|
|
// For compressed and coordinate matrix types a push_back or insert fill policy can be chosen for faster assginment:
|
|
compressed_matrix<double> I(2, 2);
|
|
I <<= fill_policy::sparse_push_back(),
|
|
0, 1, 2, 3;
|
|
std::cout << I << std::endl;
|
|
// [ 0 1 ]
|
|
// [ 2 3 ]
|
|
|
|
coordinate_matrix<double> J(2,2);
|
|
J<<=fill_policy::sparse_insert(),
|
|
1, 2, 3, 4;
|
|
std::cout << J << std::endl;
|
|
// [ 1 2 ]
|
|
// [ 3 4 ]
|
|
|
|
// A sparse matrix from another matrix works as with other types.
|
|
coordinate_matrix<double> K(3,3);
|
|
K<<=fill_policy::sparse_insert(),
|
|
J;
|
|
std::cout << K << std::endl;
|
|
// [ 1 2 0 ]
|
|
// [ 3 4 0 ]
|
|
// [ 0 0 0 ]
|
|
|
|
// Be careful this will not work:
|
|
//compressed_matrix<double> J2(4,4);
|
|
//J2<<=fill_policy::sparse_push_back(),
|
|
// J,J;
|
|
// That's because the second J2's elements
|
|
// are attempted to be assigned at positions
|
|
// that come before the elements already pushed.
|
|
// Unfortunatelly that's the only thing you can do in this case
|
|
// (or of course make a custom agorithm):
|
|
compressed_matrix<double> J2(4,4);
|
|
J2<<=fill_policy::sparse_push_back(),
|
|
J, fill_policy::sparse_insert(),
|
|
J;
|
|
|
|
std::cout << J2 << std::endl;
|
|
// [ J J ]
|
|
// [ 0 0 0 0 ]
|
|
// [ 0 0 0 0 ]
|
|
|
|
// A different traverse policy doesn't change the result, only they order it is been assigned.
|
|
coordinate_matrix<double> L(3,3);
|
|
L<<=fill_policy::sparse_insert(), traverse_policy::by_column(),
|
|
J;
|
|
std::cout << L << std::endl;
|
|
// (same as previous)
|
|
// [ 1 2 0 ]
|
|
// [ 3 4 0 ]
|
|
// [ 0 0 0 ]
|
|
|
|
typedef coordinate_matrix<double>::size_type cmst;
|
|
const cmst size = 30;
|
|
//typedef fill_policy::sparse_push_back spb;
|
|
// Although the above could have been used the following is may be faster if
|
|
// you use the policy often and for relatively small containers.
|
|
static fill_policy::sparse_push_back spb;
|
|
|
|
// A block diagonal sparse using a loop:
|
|
compressed_matrix<double> M(size, size, 4*15);
|
|
for (cmst i=0; i!=size; i+=J.size1())
|
|
M <<= spb, move_to(i,i), J;
|
|
|
|
|
|
// If typedef was used above the last expression should start
|
|
// with M <<= spb()...
|
|
|
|
// Displaying so that blocks can be easily seen:
|
|
for (unsigned int i=0; i!=M.size1(); i++) {
|
|
std::cout << M(i,0);
|
|
for (unsigned int j=1; j!=M.size2(); j++) std::cout << ", " << M(i,j);
|
|
std::cout << "\n";
|
|
}
|
|
// [ J 0 0 0 ... 0]
|
|
// [ 0 J 0 0 ... 0]
|
|
// [ 0 . . . ... 0]
|
|
// [ 0 0 ... 0 0 J]
|
|
|
|
|
|
// A "repeat" trasverser may by provided so that this becomes faster and an on-liner like:
|
|
// M <<= spb, repeat(0, size, J.size1(), 0, size, J.size1()), J;
|
|
// An alternate would be to create a :repeater" matrix and vector expression that can be used in other places as well. The latter is probably better,
|
|
return 0;
|
|
}
|
|
|