This macro shows several ways to invert a matrix .
Each method is a trade-off between accuracy of the inversion and speed. Which method to chose depends on "how well-behaved" the matrix is. This is best checked through a call to Condition(), available in each decomposition class. A second possibility (less preferred) would be to check the determinant
USAGE
This macro can be executed with Cling or ACLIC
- via the interpretor, do
- via ACLIC
root > .x invertMatrix.C+
Processing /mnt/build/workspace/root-makedoc-v610/rootspi/rdoc/src/v6-10-00-patches/tutorials/matrix/invertMatrix.C...
--------------------------------------------------------
Inversion results for a (6,6) matrix
For each inversion procedure we check the maximum size
of the off-diagonal elements of Inv(A) * A
--------------------------------------------------------
1. Use .InvertFast(&det)
Maximum off-diagonal = 8.31175e-05
Determinant = 5.3673e-18
2. Use .Invert(&det)
Maximum off-diagonal = 1.74623e-10
Determinant = 5.3673e-18
3. Use TDecompLU
Maximum off-diagonal = 1.74623e-10
Determinant = 5.3673e-18
4. Use TDecompSVD on non-square matrix
Maximum off-diagonal = 5.45697e-12
Determinant = 1.34646e-11
#include <iostream>
void invertMatrix(
Int_t msize=6)
{
if (msize < 2 || msize > 10) {
std::cout << "2 <= msize <= 10" <<std::endl;
return;
}
std::cout << "--------------------------------------------------------" <<std::endl;
std::cout << "Inversion results for a ("<<msize<<","<<msize<<") matrix" <<std::endl;
std::cout << "For each inversion procedure we check the maximum size " <<std::endl;
std::cout << "of the off-diagonal elements of Inv(A) * A " <<std::endl;
std::cout << "--------------------------------------------------------" <<std::endl;
std::cout << "1. Use .InvertFast(&det)" <<std::endl;
if (msize > 6)
std::cout << " for ("<<msize<<","<<msize<<") this is identical to .Invert(&det)" <<std::endl;
std::cout << " Maximum off-diagonal = " << U1_max_offdiag << std::endl;
std::cout << " Determinant = " << det1 << std::endl;
std::cout << "2. Use .Invert(&det)" << std::endl;
std::cout << " Maximum off-diagonal = " << U2_max_offdiag << std::endl;
std::cout << " Determinant = " << det2 << std::endl;
std::cout << "3. Use TDecompLU" << std::endl;
lu.Invert(H3);
lu.Det(d1_lu,d2_lu);
std::cout << " Maximum off-diagonal = " << U3_max_offdiag << std::endl;
std::cout << " Determinant = " << det3 << std::endl;
std::cout << "4. Use TDecompSVD on non-square matrix" << std::endl;
svd.Det(d1_svd,d2_svd);
std::cout << " Maximum off-diagonal = " << U4_max_offdiag << std::endl;
std::cout << " Determinant = " << det4 << std::endl;
}
- Author
- Eddy Offermann
Definition in file invertMatrix.C.