dense_matrix.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library 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
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



#ifndef LIBMESH_DENSE_MATRIX_H
#define LIBMESH_DENSE_MATRIX_H

// Local Includes
#include "libmesh/libmesh_common.h"
#include "libmesh/dense_matrix_base.h"

// C++ includes
#include <vector>
#include <algorithm>

namespace libMesh
{

// Forward Declarations
template <typename T> class DenseVector;



// ------------------------------------------------------------
// Dense Matrix class definition
template<typename T>
class DenseMatrix : public DenseMatrixBase<T>
{
public:

  DenseMatrix(const unsigned int new_m=0,
              const unsigned int new_n=0);

  //DenseMatrix (const DenseMatrix<T>& other_matrix);

  virtual ~DenseMatrix() {}


  virtual void zero();

  T operator() (const unsigned int i,
                const unsigned int j) const;

  T & operator() (const unsigned int i,
                  const unsigned int j);

  virtual T el(const unsigned int i,
               const unsigned int j) const { return (*this)(i,j); }

  virtual T & el(const unsigned int i,
                 const unsigned int j)     { return (*this)(i,j); }

  virtual void left_multiply (const DenseMatrixBase<T>& M2);

  template <typename T2>
  void left_multiply (const DenseMatrixBase<T2>& M2);

  virtual void right_multiply (const DenseMatrixBase<T>& M2);

  template <typename T2>
  void right_multiply (const DenseMatrixBase<T2>& M2);

  void vector_mult (DenseVector<T>& dest,
                    const DenseVector<T>& arg) const;

  template <typename T2>
  void vector_mult (DenseVector<typename CompareTypes<T,T2>::supertype>& dest,
                    const DenseVector<T2>& arg) const;

  void vector_mult_transpose (DenseVector<T>& dest,
                              const DenseVector<T>& arg) const;

  template <typename T2>
  void vector_mult_transpose (DenseVector<typename CompareTypes<T,T2>::supertype>& dest,
                              const DenseVector<T2>& arg) const;

  void vector_mult_add (DenseVector<T>& dest,
                        const T factor,
                        const DenseVector<T>& arg) const;

  template <typename T2, typename T3>
  void vector_mult_add (DenseVector<typename CompareTypes<T, typename CompareTypes<T2,T3>::supertype>::supertype>& dest,
                        const T2 factor,
                        const DenseVector<T3>& arg) const;

  void get_principal_submatrix (unsigned int sub_m, unsigned int sub_n, DenseMatrix<T>& dest) const;

  void get_principal_submatrix (unsigned int sub_m, DenseMatrix<T>& dest) const;

  DenseMatrix<T>& operator = (const DenseMatrix<T>& other_matrix);

  template <typename T2>
  DenseMatrix<T>& operator = (const DenseMatrix<T2>& other_matrix);

  void swap(DenseMatrix<T>& other_matrix);

  void resize(const unsigned int new_m,
              const unsigned int new_n);

  void scale (const T factor);


  void scale_column (const unsigned int col, const T factor);

  DenseMatrix<T>& operator *= (const T factor);

  template<typename T2, typename T3>
  typename boostcopy::enable_if_c<
    ScalarTraits<T2>::value, void >::type add (const T2 factor,
                                               const DenseMatrix<T3>& mat);

  bool operator== (const DenseMatrix<T> &mat) const;

  bool operator!= (const DenseMatrix<T> &mat) const;

  DenseMatrix<T>& operator+= (const DenseMatrix<T> &mat);

  DenseMatrix<T>& operator-= (const DenseMatrix<T> &mat);

  Real min () const;

  Real max () const;

  Real l1_norm () const;

  Real linfty_norm () const;

  void left_multiply_transpose (const DenseMatrix<T>& A);

  template <typename T2>
  void left_multiply_transpose (const DenseMatrix<T2>& A);


  void right_multiply_transpose (const DenseMatrix<T>& A);

  template <typename T2>
  void right_multiply_transpose (const DenseMatrix<T2>& A);

  T transpose (const unsigned int i,
               const unsigned int j) const;

  void get_transpose(DenseMatrix<T>& dest) const;

  std::vector<T>& get_values() { return _val; }

  const std::vector<T>& get_values() const { return _val; }

  void condense(const unsigned int i,
                const unsigned int j,
                const T val,
                DenseVector<T>& rhs)
  { DenseMatrixBase<T>::condense (i, j, val, rhs); }

  void lu_solve (const DenseVector<T>& b,
                 DenseVector<T>& x);



  template <typename T2>
  void cholesky_solve(const DenseVector<T2>& b,
                      DenseVector<T2>& x);


  void svd(DenseVector<T>& sigma);


  void svd(DenseVector<T>& sigma, DenseMatrix<T>& U, DenseMatrix<T>& VT);


  void evd(DenseVector<T>& lambda_real, DenseVector<T>& lambda_imag);


  T det();

  // void inverse();

  bool use_blas_lapack;

private:

  std::vector<T> _val;

  void _lu_decompose ();

  void _lu_back_substitute (const DenseVector<T>& b,
                            DenseVector<T>& x) const;

  void _cholesky_decompose();

  template <typename T2>
  void _cholesky_back_substitute(const DenseVector<T2>& b,
                                 DenseVector<T2>& x) const;

  enum DecompositionType {LU=0, CHOLESKY=1, LU_BLAS_LAPACK, NONE};

  DecompositionType _decomposition_type;

  enum _BLAS_Multiply_Flag {
    LEFT_MULTIPLY = 0,
    RIGHT_MULTIPLY,
    LEFT_MULTIPLY_TRANSPOSE,
    RIGHT_MULTIPLY_TRANSPOSE
  };

  void _multiply_blas(const DenseMatrixBase<T>& other,
                      _BLAS_Multiply_Flag flag);

  void _lu_decompose_lapack();

  void _svd_lapack(DenseVector<T>& sigma);

  void _svd_lapack(DenseVector<T>& sigma,
                   DenseMatrix<T>& U,
                   DenseMatrix<T>& VT);

  void _svd_helper (char JOBU,
                    char JOBVT,
                    std::vector<T>& sigma_val,
                    std::vector<T>& U_val,
                    std::vector<T>& VT_val);

  void _evd_lapack(DenseVector<T>& lambda_real, DenseVector<T>& lambda_imag);

  std::vector<int> _pivots;

  void _lu_back_substitute_lapack (const DenseVector<T>& b,
                                   DenseVector<T>& x);

  void _matvec_blas(T alpha, T beta,
                    DenseVector<T>& dest,
                    const DenseVector<T>& arg,
                    bool trans=false) const;
};





// ------------------------------------------------------------
namespace DenseMatrices
{

typedef DenseMatrix<Real> RealDenseMatrix;

typedef DenseMatrix<Complex> ComplexDenseMatrix;

}



using namespace DenseMatrices;





// ------------------------------------------------------------
// Dense Matrix member functions
template<typename T>
inline
DenseMatrix<T>::DenseMatrix(const unsigned int new_m,
                            const unsigned int new_n)
  : DenseMatrixBase<T>(new_m,new_n),
#if defined(LIBMESH_HAVE_PETSC) && defined(LIBMESH_USE_REAL_NUMBERS) && defined(LIBMESH_DEFAULT_DOUBLE_PRECISION)
    use_blas_lapack(true),
#else
    use_blas_lapack(false),
#endif
    _val(),
    _decomposition_type(NONE),
    _pivots()
{
  this->resize(new_m,new_n);
}



// FIXME[JWP]: This copy ctor has not been maintained along with
// the rest of the class...
// Can we just use the compiler-generated copy ctor here?
// template<typename T>
// inline
// DenseMatrix<T>::DenseMatrix (const DenseMatrix<T>& other_matrix)
//   : DenseMatrixBase<T>(other_matrix._m, other_matrix._n)
// {
//   _val = other_matrix._val;
// }



template<typename T>
inline
void DenseMatrix<T>::swap(DenseMatrix<T>& other_matrix)
{
  std::swap(this->_m, other_matrix._m);
  std::swap(this->_n, other_matrix._n);
  _val.swap(other_matrix._val);
  DecompositionType _temp = _decomposition_type;
  _decomposition_type = other_matrix._decomposition_type;
  other_matrix._decomposition_type = _temp;
}


template <typename T>
template <typename T2>
inline
DenseMatrix<T>&
DenseMatrix<T>::operator=(const DenseMatrix<T2>& mat)
{
  unsigned int mat_m = mat.m(), mat_n = mat.n();
  this->resize(mat_m, mat_n);
  for (unsigned int i=0; i<mat_m; i++)
    for (unsigned int j=0; j<mat_n; j++)
      (*this)(i,j) = mat(i,j);

  return *this;
}



template<typename T>
inline
void DenseMatrix<T>::resize(const unsigned int new_m,
                            const unsigned int new_n)
{
  _val.resize(new_m*new_n);

  this->_m = new_m;
  this->_n = new_n;

  // zero and set decomposition_type to NONE
  this->zero();
}



template<typename T>
inline
void DenseMatrix<T>::zero()
{
  _decomposition_type = NONE;

  std::fill (_val.begin(), _val.end(), static_cast<T>(0));
}



template<typename T>
inline
DenseMatrix<T>& DenseMatrix<T>::operator = (const DenseMatrix<T>& other_matrix)
{
  this->_m = other_matrix._m;
  this->_n = other_matrix._n;

  _val                = other_matrix._val;
  _decomposition_type = other_matrix._decomposition_type;

  return *this;
}



template<typename T>
inline
T DenseMatrix<T>::operator () (const unsigned int i,
                               const unsigned int j) const
{
  libmesh_assert_less (i*j, _val.size());
  libmesh_assert_less (i, this->_m);
  libmesh_assert_less (j, this->_n);


  //  return _val[(i) + (this->_m)*(j)]; // col-major
  return _val[(i)*(this->_n) + (j)]; // row-major
}



template<typename T>
inline
T & DenseMatrix<T>::operator () (const unsigned int i,
                                 const unsigned int j)
{
  libmesh_assert_less (i*j, _val.size());
  libmesh_assert_less (i, this->_m);
  libmesh_assert_less (j, this->_n);

  //return _val[(i) + (this->_m)*(j)]; // col-major
  return _val[(i)*(this->_n) + (j)]; // row-major
}





template<typename T>
inline
void DenseMatrix<T>::scale (const T factor)
{
  for (unsigned int i=0; i<_val.size(); i++)
    _val[i] *= factor;
}


template<typename T>
inline
void DenseMatrix<T>::scale_column (const unsigned int col, const T factor)
{
  for (unsigned int i=0; i<this->m(); i++)
    (*this)(i, col) *= factor;
}



template<typename T>
inline
DenseMatrix<T>& DenseMatrix<T>::operator *= (const T factor)
{
  this->scale(factor);
  return *this;
}



template<typename T>
template<typename T2, typename T3>
inline
typename boostcopy::enable_if_c<
  ScalarTraits<T2>::value, void >::type
DenseMatrix<T>::add (const T2 factor,
                     const DenseMatrix<T3>& mat)
{
  libmesh_assert_equal_to (this->m(), mat.m());
  libmesh_assert_equal_to (this->n(), mat.n());

  for (unsigned int i=0; i<this->m(); i++)
    for (unsigned int j=0; j<this->n(); j++)
      (*this)(i,j) += factor * mat(i,j);
}



template<typename T>
inline
bool DenseMatrix<T>::operator == (const DenseMatrix<T> &mat) const
{
  for (unsigned int i=0; i<_val.size(); i++)
    if (_val[i] != mat._val[i])
      return false;

  return true;
}



template<typename T>
inline
bool DenseMatrix<T>::operator != (const DenseMatrix<T> &mat) const
{
  for (unsigned int i=0; i<_val.size(); i++)
    if (_val[i] != mat._val[i])
      return true;

  return false;
}



template<typename T>
inline
DenseMatrix<T>& DenseMatrix<T>::operator += (const DenseMatrix<T> &mat)
{
  for (unsigned int i=0; i<_val.size(); i++)
    _val[i] += mat._val[i];

  return *this;
}



template<typename T>
inline
DenseMatrix<T>& DenseMatrix<T>::operator -= (const DenseMatrix<T> &mat)
{
  for (unsigned int i=0; i<_val.size(); i++)
    _val[i] -= mat._val[i];

  return *this;
}



template<typename T>
inline
Real DenseMatrix<T>::min () const
{
  libmesh_assert (this->_m);
  libmesh_assert (this->_n);
  Real my_min = libmesh_real((*this)(0,0));

  for (unsigned int i=0; i!=this->_m; i++)
    {
      for (unsigned int j=0; j!=this->_n; j++)
        {
          Real current = libmesh_real((*this)(i,j));
          my_min = (my_min < current? my_min : current);
        }
    }
  return my_min;
}



template<typename T>
inline
Real DenseMatrix<T>::max () const
{
  libmesh_assert (this->_m);
  libmesh_assert (this->_n);
  Real my_max = libmesh_real((*this)(0,0));

  for (unsigned int i=0; i!=this->_m; i++)
    {
      for (unsigned int j=0; j!=this->_n; j++)
        {
          Real current = libmesh_real((*this)(i,j));
          my_max = (my_max > current? my_max : current);
        }
    }
  return my_max;
}



template<typename T>
inline
Real DenseMatrix<T>::l1_norm () const
{
  libmesh_assert (this->_m);
  libmesh_assert (this->_n);

  Real columnsum = 0.;
  for (unsigned int i=0; i!=this->_m; i++)
    {
      columnsum += std::abs((*this)(i,0));
    }
  Real my_max = columnsum;
  for (unsigned int j=1; j!=this->_n; j++)
    {
      columnsum = 0.;
      for (unsigned int i=0; i!=this->_m; i++)
        {
          columnsum += std::abs((*this)(i,j));
        }
      my_max = (my_max > columnsum? my_max : columnsum);
    }
  return my_max;
}



template<typename T>
inline
Real DenseMatrix<T>::linfty_norm () const
{
  libmesh_assert (this->_m);
  libmesh_assert (this->_n);

  Real rowsum = 0.;
  for (unsigned int j=0; j!=this->_n; j++)
    {
      rowsum += std::abs((*this)(0,j));
    }
  Real my_max = rowsum;
  for (unsigned int i=1; i!=this->_m; i++)
    {
      rowsum = 0.;
      for (unsigned int j=0; j!=this->_n; j++)
        {
          rowsum += std::abs((*this)(i,j));
        }
      my_max = (my_max > rowsum? my_max : rowsum);
    }
  return my_max;
}



template<typename T>
inline
T DenseMatrix<T>::transpose (const unsigned int i,
                             const unsigned int j) const
{
  // Implement in terms of operator()
  return (*this)(j,i);
}





// template<typename T>
// inline
// void DenseMatrix<T>::condense(const unsigned int iv,
//       const unsigned int jv,
//       const T val,
//       DenseVector<T>& rhs)
// {
//   libmesh_assert_equal_to (this->_m, rhs.size());
//   libmesh_assert_equal_to (iv, jv);


//   // move the known value into the RHS
//   // and zero the column
//   for (unsigned int i=0; i<this->m(); i++)
//     {
//       rhs(i) -= ((*this)(i,jv))*val;
//       (*this)(i,jv) = 0.;
//     }

//   // zero the row
//   for (unsigned int j=0; j<this->n(); j++)
//     (*this)(iv,j) = 0.;

//   (*this)(iv,jv) = 1.;
//   rhs(iv) = val;

// }




} // namespace libMesh




#endif // LIBMESH_DENSE_MATRIX_H