Re: [ublas] determinant and norm

From: Thomas Lemaire (thomas.lemaire_at_[hidden])
Date: 2005-12-20 06:12:57

Hello Eddie,

I copy/paste my functions to compute determinant and inverse of a general 
dense matrix using lu factorisation. I am not using any bindings, but the lu 
decomposition from ublas itself, the advantage is that it works with any kind 
of dense matrix (float, double, column or row major,...), whereas Lapack 
needs column_major matrix.
My functions could be better using expression_template, but I did not spend 
much time on this.

=============== code ================

#if BOOST_VERSION < 103301 // missing includes in lu.hpp
#include "boost/numeric/ublas/vector.hpp"
#include "boost/numeric/ublas/vector_proxy.hpp"
#include "boost/numeric/ublas/matrix.hpp"
#include "boost/numeric/ublas/matrix_proxy.hpp"
#include "boost/numeric/ublas/triangular.hpp"
#include "boost/numeric/ublas/operation.hpp"
#endif

#include <boost/numeric/ublas/lu.hpp>

// forget about JFR_PRECOND macros...
  typedef boost::numeric::ublas::matrix<double> mat;
  typedef boost::numeric::ublas::identity_matrix<double> identity_mat;

      /** General matrix inversion routine.
       *  It uses lu_factorize and lu_substitute in uBLAS to invert a matrix 
       */
      template<class M1, class M2>
      void lu_inv(M1 const& m, M2& inv) {
        JFR_PRECOND(m.size1() == m.size2(),
                    "ublasExtra::lu_inv(): input matrix must be squared");
        JFR_PRECOND(inv.size1() == m.size1() && inv.size1() == m.size2(),
                    "ublasExtra::lu_inv(): invalid size for inverse matrix");

        using namespace boost::numeric::ublas;
        // create a working copy of the input
        mat mLu(m);

        // perform LU-factorization
        lu_factorize(mLu);

        // create identity matrix of "inverse"
        inv.assign(identity_mat(m.size1()));

        // backsubstitute to get the inverse
        lu_substitute<mat const, M2 >(mLu, inv);
      }

      /** General matrix determinant. 
       *  It uses lu_factorize in uBLAS.
       */
      template<class M>
      double lu_det(M const& m) {
        JFR_PRECOND(m.size1() == m.size2(),
                    "ublasExtra::lu_det: matrix must be square");

        // create a working copy of the input
        mat mLu(m);
        ublas::permutation_matrix<std::size_t> pivots(m.size1());

        lu_factorize(mLu, pivots);

        double det = 1.0;

        for (std::size_t i=0; i < pivots.size(); ++i) {
          if (pivots(i) != i)
            det *= -1.0;
          det *= mLu(i,i);
        }
        return det;
      }

=============== code ================

thomas

On Monday 19 December 2005 15:10, Eddie.Vedder_at_[hidden] wrote:
> Hello list,
>
> I have two questions:
>
> Is there a function to determine the determinant of a matrix in ublas or
> one of the bindings (lapack, atlas)? I could not find one.
>
> Is there a function to determin the "norm" which is computed by Matlabs
> norm-function which "returns the largest singular value of matrix A"? Or do
> I have to compute the singular value decomposition explicitly?
>
> Thanks for any help.
> Ed

-- 
thomas