[Home]LU Matrix Inversion

BOOST WIKI | RecentChanges | Preferences | Page List | Links List

Showing revision 11
The following code inverts the matrix input using LU-decomposition with backsubstitution of unit vectors. Reference: Numerical Recipies in C, 2nd ed., by Press, Teukolsky, Vetterling & Flannery.
 #ifndef INVERT_MATRIX_HPP
 #define INVERT_MATRIX_HPP

 // REMEMBER to update "lu.hpp" header includes from boost-CVS
 #include <boost/numeric/ublas/vector.hpp>
 #include <boost/numeric/ublas/vector_proxy.hpp>
 #include <boost/numeric/ublas/matrix.hpp>
 #include <boost/numeric/ublas/triangular.hpp>
 #include <boost/numeric/ublas/lu.hpp>
 #include <boost/numeric/ublas/io.hpp>

 namespace ublas = boost::numeric::ublas;

 /* Matrix inversion routine.
    Uses lu_factorize and lu_substitute in uBLAS to invert a matrix */
 template<class T>
 void InvertMatrix? (const ublas::matrix<T>& input, ublas::matrix<T>& inverse) {
 	using namespace boost::numeric::ublas;
 	typedef permutation_matrix<std::size_t> pmatrix;
 	// create a working copy of the input
 	matrix<T> A(input);
 	// create a permutation matrix for the LU-factorization
 	pmatrix pm(A.size1());

 	// perform LU-factorization
 	lu_factorize(A,pm);

 	// create identity matrix of "inverse"
 	inverse.assign(ublas::identity_matrix<T>(A.size1()));

 	// backsubstitute to get the inverse
 	lu_substitute(A, pm, inverse);
 }

 #endif //INVERT_MATRIX_HPP
Hope someone finds this useful. Regards, Fredrik Orderud.

Matrix inversion with Singularity Detection

This code snipet is concise and effective, but when input matrix is singular, the function simply crashes. A more complex version is at the second section of Effective UBLAS/Matrix Inversion with the capability to detect and report singularity of the input matrix.

edit: lu_factorize() returns 0 if it was successful. It returns (k+1) if it detects singularity after processing row k. So one should always check its return value.


BOOST WIKI | RecentChanges | Preferences | Page List | Links List
Edit revision 11 of this page | View other revisions | View current revision
Edited June 13, 2006 2:41 pm (diff)
Search:
Disclaimer: This site not officially maintained by Boost Developers