ceph/src/boost/libs/numeric/ublas/doc/samples/assignment_examples.cpp

   1 //
   2 //  Copyright (c) 2010 Athanasios Iliopoulos
   3 //
   4 //  Distributed under the Boost Software License, Version 1.0. (See
   5 //  accompanying file LICENSE_1_0.txt or copy at
   6 //  http://www.boost.org/LICENSE_1_0.txt)
   7 //
   8
   9 #include <boost/numeric/ublas/assignment.hpp>
  10 #include <boost/numeric/ublas/vector.hpp>
  11 #include <boost/numeric/ublas/vector_proxy.hpp>
  12 #include <boost/numeric/ublas/matrix_proxy.hpp>
  13 #include <boost/numeric/ublas/vector_sparse.hpp>
  14 #include <boost/numeric/ublas/matrix_sparse.hpp>
  15 #include <boost/numeric/ublas/io.hpp>
  16 #include <boost/numeric/ublas/matrix.hpp>
  17
  18 using namespace boost::numeric::ublas;
  19
  20 int main() {
  21         // Simple vector fill
  22     vector<double> a(3);
  23     a <<= 0, 1, 2;
  24     std::cout << a << std::endl;
  25     // [ 0 1 2]
  26
  27     // Vector from vector
  28     vector<double> b(7);
  29     b <<= a, 10, a;
  30     std::cout << b << std::endl;
  31     // [ 0 1 2 10 0 1 2]
  32
  33     // Simple matrix fill
  34     matrix<double> A(3,3);
  35     A <<= 0, 1, 2,
  36          3, 4, 5,
  37          6, 7, 8;
  38     std::cout << A << std::endl;
  39     // [ 0 1 2 ]
  40     // [ 3 4 5 ]
  41     // [ 6 7 8 ]
  42
  43     // Matrix from vector
  44     A <<= 0, 1, 2,
  45          3, 4, 5,
  46          a;
  47     std::cout << A << std::endl;
  48     // [ 0 1 2 ]
  49     // [ 3 4 5 ]
  50     // [ 0 1 2 ]
  51
  52     // Matrix from vector - column assignment
  53     A <<= move(0,2), traverse_policy::by_column(),
  54          a;
  55     std::cout << A << std::endl;
  56     // [ 0 1 0 ]
  57     // [ 3 4 1 ]
  58     // [ 0 1 2 ]
  59
  60     // Another matrix from vector example (watch the wraping);
  61     vector<double> c(9); c <<= 1, 2, 3, 4, 5, 6, 7, 8, 9;
  62     A <<= c;
  63     std::cout << A << std::endl;
  64     // [ 1 2 3 ]
  65     // [ 4 5 6 ]
  66     // [ 7 8 9 ]
  67
  68     // If for performance(Benchmarks are not definite about that) or consistency reasons you need to disable wraping:
  69     static next_row_manip endr; //This can be defined globally
  70     A <<= traverse_policy::by_row_no_wrap(),
  71             1, 2, 3, endr,
  72             4, 5, 6, endr,
  73             7, 8, 9, endr;
  74     // [ 1 2 3 ]
  75     // [ 4 5 6 ]
  76     // [ 7 8 9 ]
  77     // If by default you need to disable wraping define
  78     // BOOST_UBLAS_DEFAULT_NO_WRAP_POLICY, in the compilation options,
  79     // so that you avoid typing the "traverse_policy::by_row_no_wrap()".
  80
  81     //  Plus and minus assign:
  82     A <<= fill_policy::index_plus_assign(),
  83          3,2,1;
  84     std::cout << A << std::endl;
  85     // [ 4 4 4 ]
  86     // [ 4 5 6 ]
  87     // [ 7 8 9 ]
  88
  89     // Matrix from proxy
  90     A <<= 0, 1, 2,
  91          project(b, range(3,6)),
  92          a;
  93     std::cout << A << std::endl;
  94     // [ 0 1 2 ]
  95     // [10 0 1 ]
  96     // [ 6 7 8 ]
  97
  98     // Matrix from matrix
  99     matrix<double> B(6,6);
 100     B <<= A, A,
 101          A, A;
 102     std::cout << B << std::endl;
 103     // [ A A ]
 104     // [ A A ]
 105
 106     // Matrix range (vector is similar)
 107     B = zero_matrix<double>(6,6);
 108     matrix_range<matrix<double> > mrB (B, range (1, 4), range (1, 4));
 109     mrB <<= 1,2,3,4,5,6,7,8,9;
 110     std::cout << B << std::endl;
 111     // [ 0 0 0 0 0 0]
 112     // [ 0 1 2 3 0 0]
 113     // [ 0 4 5 6 0 0]
 114     // [ 0 0 0 0 0 0]
 115     // [ 0 0 0 0 0 0]
 116     // [ 0 0 0 0 0 0]
 117
 118     // Horizontal concatenation can be achieved using this trick:
 119     matrix<double> BH(3,9);
 120     BH <<= A, A, A;
 121     std::cout << BH << std::endl;
 122     // [ A A A]
 123
 124     // Vertical concatenation can be achieved using this trick:
 125     matrix<double> BV(9,3);
 126     BV <<= A,
 127           A,
 128           A;
 129     std::cout << BV << std::endl;
 130     // [ A ]
 131     // [ A ]
 132     // [ A ]
 133
 134     // Watch the difference when assigning matrices for different traverse policies:
 135     matrix<double> BR(9,9, 0);
 136     BR <<= traverse_policy::by_row(), // This is the default, so this might as well be omitted.
 137           A, A, A;
 138     std::cout << BR << std::endl;
 139     // [ A A A]
 140     // [ 0 0 0]
 141     // [ 0 0 0]
 142
 143     matrix<double> BC(9,9, 0);
 144     BC <<= traverse_policy::by_column(),
 145           A, A, A;
 146     std::cout << BC << std::endl;
 147     // [ A 0 0]
 148     // [ A 0 0]
 149     // [ A 0 0]
 150
 151     // The following will throw a run-time exception in debug mode (matrix mid-assignment wrap is not allowed) :
 152     // matrix<double> C(7,7);
 153     // C <<= A, A, A;
 154
 155     // Matrix from matrix with index manipulators
 156     matrix<double> C(6,6,0);
 157     C <<= A, move(3,0), A;
 158     // [ A 0 ]
 159     // [ 0 A ]
 160
 161     // A faster way for to construct this dense matrix.
 162     matrix<double> D(6,6);
 163     D <<= A, zero_matrix<double>(3,3),
 164          zero_matrix<double>(3,3), A;
 165     // [ A 0 ]
 166     // [ 0 A ]
 167
 168     // The next_row and next_column index manipulators:
 169     // note: next_row and next_column functions return
 170     // a next_row_manip and and next_column_manip object.
 171     // This is the manipulator we used earlier when we disabled
 172     // wrapping.
 173     matrix<double> E(2,4,0);
 174     E <<= 1, 2, next_row(),
 175          3, 4, next_column(),5;
 176     std::cout << E << std::endl;
 177     // [ 1 2 0 5 ]
 178     // [ 3 4 0 0 ]
 179
 180     // The begin1 (moves to the begining of the column) index manipulator, begin2 does the same for the row:
 181     matrix<double> F(2,4,0);
 182     F <<= 1, 2, next_row(),
 183          3, 4, begin1(),5;
 184     std::cout << F << std::endl;
 185     // [ 1 2 5 0 ]
 186     // [ 3 4 0 0 ]
 187
 188     // The move (relative) and move_to(absolute) index manipulators (probably the most useful manipulators):
 189     matrix<double> G(2,4,0);
 190     G <<= 1, 2, move(0,1), 3,
 191          move_to(1,3), 4;
 192     std::cout << G << std::endl;
 193     // [ 1 2 0 3 ]
 194     // [ 0 0 0 4 ]
 195
 196     // Static equivallents (faster) when sizes are known at compile time:
 197     matrix<double> Gs(2,4,0);
 198     Gs <<= 1, 2, move<0,1>(), 3,
 199          move_to<1,3>(), 4;
 200     std::cout << Gs << std::endl;
 201     // [ 1 2 0 3 ]
 202     // [ 0 0 0 4 ]
 203
 204     // Choice of traverse policy (default is "row by row" traverse):
 205
 206     matrix<double> H(2,4,0);
 207     H <<= 1, 2, 3, 4,
 208          5, 6, 7, 8;
 209     std::cout << H << std::endl;
 210     // [ 1 2 3 4 ]
 211     // [ 5 6 7 8 ]
 212
 213     H <<= traverse_policy::by_column(),
 214         1, 2, 3, 4,
 215         5, 6, 7, 8;
 216     std::cout << H << std::endl;
 217     // [ 1 3 5 7 ]
 218     // [ 2 4 6 8 ]
 219
 220     // traverse policy can be changed mid assignment if desired.
 221      matrix<double> H1(4,4,0);
 222      H1 <<= 1, 2, 3, traverse_policy::by_column(), 1, 2, 3;
 223
 224     std::cout << H << std::endl;
 225     // [1 2 3 1]
 226     // [0 0 0 2]
 227     // [0 0 0 3]
 228     // [0 0 0 0]
 229
 230     // note: fill_policy and traverse_policy are namespaces, so you can use them
 231     // by a using statement.
 232
 233     // For compressed and coordinate matrix types a push_back or insert fill policy can be chosen for faster assginment:
 234     compressed_matrix<double> I(2, 2);
 235     I <<=    fill_policy::sparse_push_back(),
 236             0, 1, 2, 3;
 237     std::cout << I << std::endl;
 238     // [ 0 1 ]
 239     // [ 2 3 ]
 240
 241     coordinate_matrix<double> J(2,2);
 242     J<<=fill_policy::sparse_insert(),
 243         1, 2, 3, 4;
 244     std::cout << J << std::endl;
 245     // [ 1 2 ]
 246     // [ 3 4 ]
 247
 248     // A sparse matrix from another matrix works as with other types.
 249     coordinate_matrix<double> K(3,3);
 250     K<<=fill_policy::sparse_insert(),
 251         J;
 252     std::cout << K << std::endl;
 253     // [ 1 2 0 ]
 254     // [ 3 4 0 ]
 255     // [ 0 0 0 ]
 256
 257     // Be careful this will not work:
 258     //compressed_matrix<double> J2(4,4);
 259     //J2<<=fill_policy::sparse_push_back(),
 260      //   J,J;
 261     // That's because the second J2's elements
 262     // are attempted to be assigned at positions
 263     // that come before the elements already pushed.
 264     // Unfortunatelly that's the only thing you can do in this case
 265     // (or of course make a custom agorithm):
 266     compressed_matrix<double> J2(4,4);
 267     J2<<=fill_policy::sparse_push_back(),
 268         J, fill_policy::sparse_insert(),
 269         J;
 270
 271     std::cout << J2 << std::endl;
 272     // [  J   J  ]
 273     // [ 0 0 0 0 ]
 274     // [ 0 0 0 0 ]
 275
 276     // A different traverse policy doesn't change the result, only they order it is been assigned.
 277     coordinate_matrix<double> L(3,3);
 278     L<<=fill_policy::sparse_insert(), traverse_policy::by_column(),
 279         J;
 280     std::cout << L << std::endl;
 281     // (same as previous)
 282     // [ 1 2 0 ]
 283     // [ 3 4 0 ]
 284     // [ 0 0 0 ]
 285
 286     typedef coordinate_matrix<double>::size_type cmst;
 287     const cmst size = 30;
 288     //typedef fill_policy::sparse_push_back spb;
 289     // Although the above could have been used the following is may be faster if
 290     //  you use the policy often and for relatively small containers.
 291     static fill_policy::sparse_push_back spb;
 292
 293     // A block diagonal sparse using a loop:
 294     compressed_matrix<double> M(size, size, 4*15);
 295     for (cmst i=0; i!=size; i+=J.size1())
 296         M <<= spb, move_to(i,i), J;
 297
 298
 299     // If typedef was used above the last expression should start
 300     // with M <<= spb()...
 301
 302     // Displaying so that blocks can be easily seen:
 303     for (unsigned int i=0; i!=M.size1(); i++) {
 304         std::cout << M(i,0);
 305         for (unsigned int j=1; j!=M.size2(); j++) std::cout << ", " << M(i,j);
 306         std::cout << "\n";
 307     }
 308     // [ J 0 0 0 ... 0]
 309     // [ 0 J 0 0 ... 0]
 310     // [ 0 . . . ... 0]
 311     // [ 0 0 ... 0 0 J]
 312
 313
 314     // A "repeat" trasverser may by provided so that this becomes faster and an on-liner like:
 315     // M <<= spb, repeat(0, size, J.size1(), 0, size, J.size1()), J;
 316     // 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,
 317     return 0;
 318 }
 319