MTK++ Latest version: 0.2.0

inverseMatrix.h
Go to the documentation of this file.
00001 #ifndef INVERT_MATRIX_HPP
00002 #define INVERT_MATRIX_HPP
00003 
00004 // REMEMBER to update "lu.hpp" header includes from boost-CVS
00005 #include <boost/numeric/ublas/vector.hpp>
00006 #include <boost/numeric/ublas/vector_proxy.hpp>
00007 #include <boost/numeric/ublas/matrix.hpp>
00008 #include <boost/numeric/ublas/triangular.hpp>
00009 #include <boost/numeric/ublas/lu.hpp>
00010 #include <boost/numeric/ublas/io.hpp>
00011 
00012 namespace ublas = boost::numeric::ublas;
00013 
00014 /* Matrix inversion routine.
00015    Uses lu_factorize and lu_substitute in uBLAS to invert a matrix
00016 */
00017 template<class T>
00018 void InvertMatrix(const ublas::matrix<T>& input,
00019                         ublas::matrix<T>& inverse)
00020 {
00021     using namespace boost::numeric::ublas;
00022     typedef permutation_matrix<std::size_t> pmatrix;
00023     // create a working copy of the input
00024     matrix<T> A(input);
00025     // create a permutation matrix for the LU-factorization
00026     pmatrix pm(A.size1());
00027 
00028     // perform LU-factorization
00029     lu_factorize(A,pm);
00030 
00031     // create identity matrix of "inverse"
00032     inverse.assign(ublas::identity_matrix<T>(A.size1()));
00033 
00034     // backsubstitute to get the inverse
00035     lu_substitute(A, pm, inverse);
00036 }
00037 
00038 #endif //INVERT_MATRIX_HPP
00039 
00040 /*
00041 #include <boost/numeric/ublas/matrix.hpp>
00042 #include <boost/numeric/ublas/matrix_proxy.hpp>
00043 #include <boost/numeric/ublas/io.hpp>
00044 #include <boost/numeric/ublas/matrix_expression.hpp>
00045 
00046 #include <iostream>
00047 #include <fstream>
00048 #include <vector>
00049 
00050  //
00051  // Invert a matrix via gauss-jordan algorithm (PARTIAL PIVOT)
00052  //
00053  // @param m The matrix to invert. Must be square.
00054  // @param singular If the matrix was found to be singular, then this
00055  //        is set to true, else set to false.
00056  // @return If singular is false, then the inverted matrix is returned.
00057  //         Otherwise it contains random values.
00058  //
00059 template<class T>
00060 boost::numeric::ublas::matrix<T>
00061 gjinverse(const boost::numeric::ublas::matrix<T> &m, bool &singular)
00062 {
00063     using namespace boost::numeric::ublas;
00064     const int size = m.size1();
00065 
00066     // Cannot invert if non-square matrix or 0x0 matrix.
00067     // Report it as singular in these cases, and return 
00068     // a 0x0 matrix.
00069     if (size != m.size2() || size == 0) {
00070       singular = true;
00071       matrix<T> A(0,0);
00072       return A;
00073     }
00074 
00075     // Handle 1x1 matrix edge case as general purpose 
00076     // inverter below requires 2x2 to function properly.
00077     if (size == 1) {
00078       matrix<T> A(1, 1);
00079       if (m(0,0) == 0.0) {
00080         singular = true;
00081         return A;
00082       }
00083       singular = false;
00084       A(0,0) = 1/m(0,0);
00085          return A;
00086     }
00087 
00088     // Create an augmented matrix A to invert. Assign the
00089     // matrix to be inverted to the left hand side and an
00090     // identity matrix to the right hand side.
00091     matrix<T> A(size, 2*size);
00092     matrix_range<matrix<T> > Aleft(A, range(0, size), range(0, size));
00093     Aleft = m;
00094     matrix_range<matrix<T> > Aright(A, range(0, size), range(size, 2*size));
00095     Aright = identity_matrix<T>(size);
00096 
00097     // Doing partial pivot
00098     for (int k = 0; k < size; k++) {
00099       // Swap rows to eliminate zero diagonal elements.
00100       for (int k = 0; k < size; k++) {
00101         if ( A(k,k) == 0 ) { // XXX: test for "small" instead
00102           // Find a row(l) to swap with row(k)
00103           int l = -1;
00104           for (int i = k+1; i < size; i++) {
00105             if ( A(i,k) != 0 ) {
00106               l = i; 
00107               break;
00108             }
00109           }
00110 
00111           // Swap the rows if found
00112           if ( l < 0 ) {
00113             std::cerr << "Error:" <<  __FUNCTION__ << ":"
00114                       << "Input matrix is singular, because cannot find"
00115                       << " a row to swap while eliminating zero-diagonal.";
00116             singular = true;
00117             return Aleft;
00118           }
00119           else {
00120             matrix_row<matrix<T> > rowk(A, k);
00121             matrix_row<matrix<T> > rowl(A, l);
00122             rowk.swap(rowl);
00123 
00124 #if defined(DEBUG) || !defined(NDEBUG)
00125                  std::cerr << __FUNCTION__ << ":"
00126                            << "Swapped row " << k << " with row " << l 
00127                            << ":" << A << "\n";
00128  #endif
00129           }
00130         }
00131       }
00132 
00134       // normalize the current row
00135       for (int j = k+1; j < 2*size; j++)
00136         A(k,j) /= A(k,k);
00137         A(k,k) = 1;
00138 
00139         // normalize other rows
00140         for (int i = 0; i < size; i++) {
00141           if ( i != k ) { // other rows  // FIX: PROBLEM HERE
00142             if ( A(i,k) != 0 ) {
00143               for (int j = k+1; j < 2*size; j++)
00144                 A(i,j) -= A(k,j) * A(i,k);
00145                 A(i,k) = 0;
00146               }
00147             }
00148           }
00149 
00150  #if defined(DEBUG) || !defined(NDEBUG)
00151          std::cerr << __FUNCTION__ << ":"
00152                    << "GJ row " << k << " : " << A << "\n";
00153  #endif
00154      }
00155 
00156      singular = false;
00157      return Aright;
00158 }
00159 */

Generated on Fri Dec 23 2011 09:28:50 for MTK++ by Doxygen 1.7.5.1