MatrixXf A = MatrixXf::Random(4, 4);
MatrixXf B = MatrixXf::Random(4, 4);
RealQZ<MatrixXf> qz(4);  // preallocate space for 4x4 matrices
qz.compute(A, B);        // A = Q S Z,  B = Q T Z

// print original matrices and result of decomposition
cout << "A:\n"
     << A << "\n"
     << "B:\n"
     << B << "\n";
cout << "S:\n"
     << qz.matrixS() << "\n"
     << "T:\n"
     << qz.matrixT() << "\n";
cout << "Q:\n"
     << qz.matrixQ() << "\n"
     << "Z:\n"
     << qz.matrixZ() << "\n";

// verify precision
cout << "\nErrors:"
     << "\n|A-QSZ|: " << (A - qz.matrixQ() * qz.matrixS() * qz.matrixZ()).norm()
     << ", |B-QTZ|: " << (B - qz.matrixQ() * qz.matrixT() * qz.matrixZ()).norm()
     << "\n|QQ* - I|: " << (qz.matrixQ() * qz.matrixQ().adjoint() - MatrixXf::Identity(4, 4)).norm()
     << ", |ZZ* - I|: " << (qz.matrixZ() * qz.matrixZ().adjoint() - MatrixXf::Identity(4, 4)).norm() << "\n";