linalg.h

Go to the documentation of this file.

#ifndef LINALG_H
#define LINALG_H

#include <boost/numeric/ublas/banded.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <cassert>
#include "mimasconfig.h"
#include "angle.h"
#include "mimasexception.h"

namespace mimas {

template< typename T >
boost::numeric::ublas::vector< T > unit
   ( const boost::numeric::ublas::vector< T > &x )
{
  assert( boost::numeric::ublas::norm_2( x ) > 0 );
  return x / boost::numeric::ublas::norm_2( x );
}

template< typename T >
T scalar_cross_product( const boost::numeric::ublas::vector< T > &a, 
                       const boost::numeric::ublas::vector< T > &b )
{
  assert( a.size() == 2 );
  assert( b.size() == 2 );
  return a(0) * b(1) - a(1) * b(0);
}

template< typename T >
angle getAngle( const boost::numeric::ublas::vector< T > &a, 
                const boost::numeric::ublas::vector< T > &b )
{
  return angle( scalar_cross_product< T >( a, b ), inner_prod( a, b ) );
}

template< typename T >
boost::numeric::ublas::vector< T > rotate
   ( const boost::numeric::ublas::vector< T > &v, const angle &a );

template< typename T >
boost::numeric::ublas::vector< T > crossProduct
   ( boost::numeric::ublas::vector< T > &a,
     boost::numeric::ublas::vector< T > &b );

template< typename T >
boost::numeric::ublas::vector< T > rodrigues
   ( boost::numeric::ublas::vector< T > const &u,
     boost::numeric::ublas::vector< T > const &v,
     double angle);

template< typename T >
double determinant (boost::numeric::ublas::matrix< T > const &M);


#ifdef HAVE_LIBLAPACK

template< typename T >
struct gglse_t
{
  typedef boost::numeric::ublas::vector< T > Vector;
  typedef boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > Matrix;
  Vector operator()( const Matrix &_A, const Matrix &_B,
                     const Vector &_c, const Vector &_d ) const;
};

template< typename T >
boost::numeric::ublas::vector< T > gglse
  ( const boost::numeric::ublas::matrix
      < T, boost::numeric::ublas::column_major > &A,
    const boost::numeric::ublas::matrix
      < T, boost::numeric::ublas::column_major > &B,
    const boost::numeric::ublas::vector< T > &c,
    const boost::numeric::ublas::vector< T > &d )
{ return gglse_t< T >()( A, B, c, d ); }

template< typename T >
struct inv_t
{
  typedef boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > Matrix;
  Matrix operator()( const Matrix &A ) const;
};

template< typename T >
boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > inv
  ( const boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > &A )
{ return inv_t< T >()( A ); }

template< typename T >
struct gels_t
{
  typedef boost::numeric::ublas::vector< T > Vector;
  typedef boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > Matrix;
  Vector operator()( const Matrix &_A, const Vector &_b ) const;
};

template< typename T >
boost::numeric::ublas::vector< T > gels
  ( const boost::numeric::ublas::matrix
      < T, boost::numeric::ublas::column_major > &A,
    const boost::numeric::ublas::vector< T > &b )
{ return gels_t< T >()( A, b ); }

template< typename T >
struct gesvd_t
{
  typedef boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > Matrix;
  typedef boost::numeric::ublas::banded_matrix< T > BandedMatrix;
  BandedMatrix operator()( const Matrix &_A, Matrix *_U, Matrix *_Vt )
       throw(mimasexception);
};

template< typename T >
boost::numeric::ublas::banded_matrix< T > gesvd
  ( const boost::numeric::ublas::matrix
      < T, boost::numeric::ublas::column_major > &A,
    boost::numeric::ublas::matrix
      < T, boost::numeric::ublas::column_major > *U,
    boost::numeric::ublas::matrix
      < T, boost::numeric::ublas::column_major > *Vt )
{ return gesvd_t< T >()( A, U, Vt ); }

template< typename T >
struct gesdd_t
{
  typedef boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > Matrix;
  typedef boost::numeric::ublas::banded_matrix< T > BandedMatrix;
  BandedMatrix operator()( const Matrix &_A, Matrix *_U, Matrix *_Vt )
       throw (mimasexception);
};

template< typename T >
boost::numeric::ublas::banded_matrix< T > gesdd
  ( const boost::numeric::ublas::matrix
      < T, boost::numeric::ublas::column_major > &A,
    boost::numeric::ublas::matrix
      < T, boost::numeric::ublas::column_major > *U,
    boost::numeric::ublas::matrix
      < T, boost::numeric::ublas::column_major > *Vt )
{ return gesdd_t< T >()( A, U, Vt ); }

template< typename T >
struct syev_t
{
  typedef boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > Matrix;
  typedef boost::numeric::ublas::symmetric_matrix< T, boost::numeric::ublas::upper > SymmetricMatrix;
  typedef boost::numeric::ublas::diagonal_matrix< T > EigenValues;
  EigenValues operator()( const SymmetricMatrix &_A, Matrix *_E ) const
    throw (mimasexception);
};

template< typename T >
boost::numeric::ublas::diagonal_matrix< T > syev
  ( const boost::numeric::ublas::symmetric_matrix
      < T, boost::numeric::ublas::upper > &A,
    boost::numeric::ublas::matrix
      < T, boost::numeric::ublas::column_major > *E )
{ return syev_t< T >()( A, E ); }


template< typename T >
struct pptrf_t
{
  typedef boost::numeric::ublas::triangular_matrix 
    < T, boost::numeric::ublas::upper, boost::numeric::ublas::column_major > TriangularMatrix;
  typedef boost::numeric::ublas::symmetric_matrix 
    < T, boost::numeric::ublas::upper, boost::numeric::ublas::column_major > SymmetricMatrix;

  TriangularMatrix operator()( const SymmetricMatrix &_A )
       throw(mimasexception);
};

template< typename T >
boost::numeric::ublas::triangular_matrix 
 < T, boost::numeric::ublas::upper, boost::numeric::ublas::column_major  > pptrf
  ( const boost::numeric::ublas::symmetric_matrix
      < T, boost::numeric::ublas::upper, boost::numeric::ublas::column_major  > &A)
{ return pptrf_t< T >()( A ); }


#endif


};

#include "linalg.tcc"

#endif

---