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