matrix.hpp

dysii是一款非常出色的滤波函数库

#ifndef INDII_ML_AUX_MATRIX_HPP
#define INDII_ML_AUX_MATRIX_HPP

#include "vector.hpp"

#include "boost/numeric/ublas/matrix.hpp"
#include "boost/numeric/ublas/matrix_sparse.hpp"
#include "boost/numeric/ublas/matrix_proxy.hpp"
#include "boost/numeric/ublas/symmetric.hpp"
#include "boost/numeric/ublas/triangular.hpp"
#include "boost/numeric/ublas/banded.hpp"
#include "boost/numeric/ublas/io.hpp"
#include "boost/numeric/ublas/storage.hpp"

namespace indii {
  namespace ml {
    namespace aux {

      /**
       * General matrix.
       */
      typedef boost::numeric::ublas::matrix<double,
          boost::numeric::ublas::column_major,
	  boost::numeric::ublas::unbounded_array<double> > matrix;

      /**
       * Symmetric matrix.
       */
      typedef boost::numeric::ublas::symmetric_matrix<double,
          boost::numeric::ublas::lower, boost::numeric::ublas::column_major,
	  boost::numeric::ublas::unbounded_array<double> > symmetric_matrix;

      /**
       * Lower triangular matrix.
       */
      typedef boost::numeric::ublas::triangular_matrix<double,
          boost::numeric::ublas::lower, boost::numeric::ublas::column_major,
	  boost::numeric::ublas::unbounded_array<double> >
          lower_triangular_matrix;

      /**
       * Upper triangular matrix.
       */
      typedef boost::numeric::ublas::triangular_matrix<double,
          boost::numeric::ublas::upper, boost::numeric::ublas::column_major,
          boost::numeric::ublas::unbounded_array<double> >
          upper_triangular_matrix;

      /**
       * Identity matrix.
       */
      typedef boost::numeric::ublas::identity_matrix<double> identity_matrix;
      
      /**
       * Zero matrix.
       */
      typedef boost::numeric::ublas::zero_matrix<double> zero_matrix;

      /**
       * Scalar matrix.
       */
      typedef boost::numeric::ublas::scalar_matrix<double> scalar_matrix;

      /**
       * Banded matrix.
       */
      typedef boost::numeric::ublas::banded_matrix<double> banded_matrix;

      /**
       * Sparse matrix.
       */
      typedef boost::numeric::ublas::mapped_matrix<double,
          boost::numeric::ublas::column_major> sparse_matrix;
      //typedef boost::numeric::ublas::compressed_matrix<double,
      //    boost::numeric::ublas::column_major> sparse_matrix;
      // ^ seems to cause segfaults in some situations, perhaps when
      //   assertions are disabled?

      /**
       * Projection matrix.
       */
      typedef boost::numeric::ublas::mapped_matrix<short int,
          boost::numeric::ublas::column_major> projection_matrix;

      /**
       * Inverse of a square matrix.
       *
       * @param A matrix to invert.
       * @param AI matrix into which to write the inverse of A.
       */
      void inv(matrix& A, matrix& AI);

      /**
       * Diagonal of a square matrix.
       *
       * @param A square matrix.
       */
      template <class T>
      boost::numeric::ublas::matrix_vector_range<T> diag(T& A);

    }
  }
}

template <class T>
boost::numeric::ublas::matrix_vector_range<T> indii::ml::aux::diag(T& A) {
  /* pre-condition */
  assert (A.size1() == A.size2());

  const boost::numeric::ublas::range step(0, A.size1());

  return boost::numeric::ublas::matrix_vector_range<T>(A, step, step);
}

#endif