Benoit Jacob | 61158b1 | 2007-09-07 08:18:21 +0000 | [diff] [blame^] | 1 | #include"../src/All.h" |
| 2 | |
| 3 | using namespace std; |
| 4 | using namespace Eigen; |
| 5 | |
| 6 | int main(int, char **) |
| 7 | { |
| 8 | Matrix<double,2,2> m; // 2x2 fixed-size matrix with uninitialized entries |
| 9 | m(0,0) = 1; |
| 10 | m(0,1) = 2; |
| 11 | m(1,0) = 3; |
| 12 | m(1,1) = 4; |
| 13 | cout << "Here is a 2x2 matrix m:" << endl << m << endl; |
| 14 | cout << "Let us now build a 4x4 matrix m2 by assembling together four 2x2 blocks." << endl; |
| 15 | MatrixX<double> m2(4, 4); // dynamic matrix with initial size 4x4 and uninitialized entries |
| 16 | // notice how we are mixing fixed-size and dynamic-size types. |
| 17 | |
| 18 | cout << "In the top-left block, we put the matrix m shown above." << endl; |
| 19 | // here we need to use .xpr() to allow write access to the block. |
| 20 | m2.xpr().block(0,1,0,1) = m; |
| 21 | cout << "In the bottom-left block, we put the matrix m*m, which is:" << endl << m*m << endl; |
| 22 | m2.xpr().block(2,3,0,1) = m * m; |
| 23 | cout << "In the top-right block, we put the matrix m+m, which is:" << endl << m+m << endl; |
| 24 | m2.xpr().block(0,1,2,3) = m + m; |
| 25 | cout << "In the bottom-right block, we put the matrix m-m, which is:" << endl << m-m << endl; |
| 26 | m2.xpr().block(2,3,2,3) = m - m; |
| 27 | cout << "Now the 4x4 matrix m2 is:" << endl << m2 << endl; |
| 28 | |
| 29 | // here we don't need to use .xpr() because we only need read access. |
| 30 | cout << "The central 2x2 block of m2 is:" << endl << m2.block(1,2,1,2) << endl; |
| 31 | cout << "Row 0 of m2, written as a column vector, is:" << endl << m2.row(0) << endl; |
| 32 | cout << "Column 1 of m2 is:" << endl << m2.col(1) << endl; |
| 33 | cout << "The matrix m2 with row 0 and column 1 removed is:" << endl << m2.minor(0,1) << endl; |
| 34 | |
| 35 | cout << endl << "Now let us study a tricky issue." << endl; |
| 36 | cout << "Recall that the matrix product m*m is:" << endl << m*m << endl; |
| 37 | cout << "We want to store that into m, i.e. do \"m = m * m;\"" << endl; |
| 38 | cout << "Here we must be very careful. For if we do \"m = m * m;\"," << endl |
| 39 | << "the matrix m becomes" << endl; |
| 40 | Matrix<double,2,2> m_save = m; |
| 41 | m = m * m; // the bogus operation |
| 42 | cout << m << "," << endl; |
| 43 | cout << "which is not what was wanted!" << endl |
| 44 | << "Explanation: because of the way expression templates work, the matrix m gets" << endl |
| 45 | << "overwritten _while_ the matrix product m * m is being computed." << endl |
| 46 | << "This is the counterpart of eliminating temporary objects!" << endl |
| 47 | << "Anyway, if you want to store m * m into m, you can do this:" << endl |
| 48 | << " m.alias() = m * m;" << endl; |
| 49 | m = m_save; |
| 50 | m.alias() = m * m; |
| 51 | cout << "And m is now:" << endl << m << endl << "as was expected." << endl; |
| 52 | cout << "To make your life easier, operator*= between matrices automatically" << endl |
| 53 | << "creates an alias. So a *= b is equivalent to a.alias() *= b," << endl |
| 54 | << "when a and b are matrices. So, coming back to the original matrix m," << endl |
| 55 | << "if we do m *= m, the matrix m becomes:" << endl; |
| 56 | m = m_save; |
| 57 | m *= m; |
| 58 | cout << m << endl << "as expected." << endl; |
| 59 | return 0; |
| 60 | } |