| Brian Silverman | 72890c2 | 2015-09-19 14:37:37 -0400 | [diff] [blame^] | 1 | MatrixXf A = MatrixXf::Random(4,4); |
| 2 | MatrixXf B = MatrixXf::Random(4,4); |
| 3 | RealQZ<MatrixXf> qz(4); // preallocate space for 4x4 matrices |
| 4 | qz.compute(A,B); // A = Q S Z, B = Q T Z |
| 5 | |
| 6 | // print original matrices and result of decomposition |
| 7 | cout << "A:\n" << A << "\n" << "B:\n" << B << "\n"; |
| 8 | cout << "S:\n" << qz.matrixS() << "\n" << "T:\n" << qz.matrixT() << "\n"; |
| 9 | cout << "Q:\n" << qz.matrixQ() << "\n" << "Z:\n" << qz.matrixZ() << "\n"; |
| 10 | |
| 11 | // verify precision |
| 12 | cout << "\nErrors:" |
| 13 | << "\n|A-QSZ|: " << (A-qz.matrixQ()*qz.matrixS()*qz.matrixZ()).norm() |
| 14 | << ", |B-QTZ|: " << (B-qz.matrixQ()*qz.matrixT()*qz.matrixZ()).norm() |
| 15 | << "\n|QQ* - I|: " << (qz.matrixQ()*qz.matrixQ().adjoint() - MatrixXf::Identity(4,4)).norm() |
| 16 | << ", |ZZ* - I|: " << (qz.matrixZ()*qz.matrixZ().adjoint() - MatrixXf::Identity(4,4)).norm() |
| 17 | << "\n"; |