ublasTools.h

Go to the documentation of this file.

// ----------------------------------------------------------------------------
//
// $Id$
//
// Copyright 2008, 2009, 2010, 2011, 2012  Antonio Franchi and Paolo Stegagno
//
// This file is part of MIP.
//
// MIP is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// MIP is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with MIP. If not, see <http://www.gnu.org/licenses/>.
//
// Contact info: antonio.franchi@tuebingen.mpg.de stegagno@diag.uniroma1.it
//
// ----------------------------------------------------------------------------
 
#ifndef __UBLAS_TOOLS_H__
#define __UBLAS_TOOLS_H__
 
 
 // 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.
 //
 // http://www.crystalclearsoftware.com/cgi-bin/boost_wiki/wiki.pl?action=browse&diff=1&id=LU_Matrix_Inversion
 // Hope someone finds this useful. Regards, Fredrik Orderud.

 #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>
 bool invert (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
   int res = lu_factorize(A,pm);
   if( res != 0 ) return false;
 
   // create identity matrix of "inverse"
   inverse.assign(ublas::identity_matrix<T>(A.size1()));
 
   // backsubstitute to get the inverse
   lu_substitute(A, pm, inverse);
 
   return true;
 }
 
 template<class T>
 boost::numeric::ublas::matrix<T>
 invert(const boost::numeric::ublas::matrix<T> &m, bool &is_singular)
 { 
    ublas::matrix<T> inv_m (m.size1(), m.size2());
    is_singular = ! invert (m, inv_m);
    return inv_m;
 }
 // http://archives.free.net.ph/message/20080909.064313.59c122c4.fr.html 
 template<class matrix_T>
 double determinant(ublas::matrix_expression<matrix_T> const& mat_r) {
   double det = 1.0;
   matrix_T mLu(mat_r());
   ublas::permutation_matrix<std::size_t> pivots(mat_r().size1() );
   bool is_singular = lu_factorize(mLu, pivots);
   if (!is_singular) {
     for (std::size_t i = 0; i < pivots.size(); ++i) {
       if (pivots(i) != i) {
         det *= -1.0;
       }
       det *= mLu(i,i);
     }
   } else {
     det = 0.0;
   }
   return det;
 } 
 
#endif