An Object Finite Element Library

DMatrix.h

Go to the documentation of this file.

/*==============================================================================

                                    O  F  E  L  I

                           Object  Finite  Element  Library

  ==============================================================================

   Copyright (C) 1998 - 2008 Rachid Touzani

   This program 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; Version 2 of the License.

   This program 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 this program; if not, write to the :

   Free Software Foundation
   Inc., 59 Temple Place - Suite 330
   Boston, MA  02111-1307, USA

  ==============================================================================

               Definition of Template Class DMatrix for dense matrices

  ==============================================================================*/


#ifndef __DMATRIX_H
#define __DMATRIX_H

#include "Matrix.h"


namespace OFELI {

template<class T_> class DMatrix;

#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_>
struct ErrorInDMatrix {
   void Message(const char *file, int line, int code, int p1=0, int p2=0, T_ p3=T_(0)) {
     cerr << "\n*** Fatal Error in OFELI ***\n";
     cerr << "File : " << file << ", line : " << line << "\nIn DMatrix<>::";
     switch (code) {

        case 51: 
           cerr << "Set(i,j,x) : Illegal pair of indices : " << p1 << "," << p2 << endl;
           break;

        case 61: 
           cerr << "Factor() : The " << p1 << "-th pivot = " << p3 << " is too small\n";
           break;

        case 62: 
           cerr << "Factor() : Can't factorize a rectangle matrix\n";
           break;

        case 71: 
           cerr << "Solve(b) : The " << p1 << "-th diagonal entry = " << p3 << " is too small\n";
           break;

        case 72: 
           cerr << "Solve(b) : Warning : Matrix has not been factorized\n";
           return;

        case 73: 
           cerr << "Solve(b) : Can't solve with a rectangle matrix\n";
           return;
     }
     cerr << "Program stops" << endl;
     exit(1);
   }
};
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

template<class T_>
class DMatrix : public Matrix<T_>
{

   using Matrix<T_>::_rmin;
   using Matrix<T_>::_cmin;
   using Matrix<T_>::_rmax;
   using Matrix<T_>::_cmax;
   using Matrix<T_>::_nb_rows;
   using Matrix<T_>::_nb_cols;
   using Matrix<T_>::_size;
   using Matrix<T_>::_length;
   using Matrix<T_>::_a;
   using Matrix<T_>::_diag;
   using Matrix<T_>::_ch;
   using Matrix<T_>::_zero;


 public:

/*----------------------------  BASIC OPERATORS  ----------------------------*/

    DMatrix();

    DMatrix(size_t nr);

    DMatrix(size_t nr, size_t nc);

    DMatrix(Vect<T_> &v);

    DMatrix(const DMatrix<T_> &m);

    ~DMatrix();

#ifndef DOXYGEN_SHOULD_SKIP_THIS
    virtual void setMesh(const class Mesh &mesh, size_t dof=0, size_t type=0);
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

#ifndef DOXYGEN_SHOULD_SKIP_THIS
    void setMesh(size_t dof, const class Mesh &mesh, int code=0);
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

#ifndef DOXYGEN_SHOULD_SKIP_THIS
    void setMesh(size_t dof, size_t nb_eq, const class Mesh &mesh);
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

    void setDiag();

    void Resize(size_t size);

    void Resize(size_t nr, size_t nc);

    Vect<T_> getColumn(size_t j) const;

    Vect<T_> getRow(size_t i) const;

    void Set(size_t i, size_t j, const T_ &x);

    void MultAdd(T_ a, const Vect<T_> &x, Vect<T_> &y) const;

    void MultAdd(const Vect<T_> &x, Vect<T_> &y) const;

    void Mult(const Vect<T_> &x, Vect<T_> &y) const;

    void TMult(const Vect<T_> &x, Vect<T_> &y) const;

    void Add(size_t i, size_t j, const T_ &x);

    T_ operator()(size_t i, size_t j) const;

    T_ & operator()(size_t i, size_t j);

    const T_ operator()(size_t i) const;

    int Factor();

    int Solve(Vect<T_> &b);

    int Solve(const Vect<T_> &b, Vect<T_> &x);

    DMatrix & operator=(DMatrix<T_> &m);

    DMatrix & operator+=(const DMatrix<T_> &m);

    DMatrix & operator-=(const DMatrix<T_> &m);

    DMatrix & operator=(const T_ &x);

    DMatrix & operator*=(const T_ &x);

    DMatrix & operator+=(const T_ &x);

    DMatrix & operator-=(const T_ &x);

    size_t getColInd(size_t i) const;

#ifndef DOXYGEN_SHOULD_SKIP_THIS
    size_t getRowPtr(size_t i) const;
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

    T_ *getArray() const;

    T_ getEntry(size_t i, size_t j) const;

    void setPrintView(size_t rmin, size_t rmax, size_t cmin, size_t cmax);

private:

   mutable ErrorInDMatrix<T_> _e;
   bool _fact;
   void set_(size_t i, size_t j, const T_ &x) { _a[_nb_cols*(i-1)+j-1] = x; }
};


//                         I M P L E M E N T A T I O N                       //


template<class T_>
DMatrix<T_>::DMatrix() : _fact(false)
{
   _length = _size = _nb_rows = _nb_cols = 0;
#ifdef _OFELI_DEBUG
   std::clog << "An instance of class DMatrix is constructed.\n";
   std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
#endif
}


template<class T_>
DMatrix<T_>::DMatrix(size_t nr)
{
   Resize(nr);
   _rmin = 1;
   _rmax = _nb_rows;
   _cmin = 1;
   _cmax = _nb_cols;
#ifdef _OFELI_DEBUG
   std::clog << "An instance of class DMatrix is constructed.\n";
   std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
#endif
}


template<class T_>
DMatrix<T_>::DMatrix(size_t nr, size_t nc)
{
   Resize(nr,nc);
   _rmin = 1;
   _rmax = _nb_rows;
   _cmin = 1;
   _cmax = _nb_cols;
#ifdef _OFELI_DEBUG
   std::clog << "An instance of class DMatrix is constructed.\n";
   std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
#endif
}


template<class T_>
DMatrix<T_>::DMatrix(Vect<T_> &v)
{
   _length = v.getSize();
   _nb_rows = _nb_cols = _size = sqrt(double(_length));
   _rmin = 1;
   _rmax = _nb_rows;
   _cmin = 1;
   _cmax = _nb_cols;
   _a = v;
}


template<class T_>
DMatrix<T_>::DMatrix(const DMatrix<T_> &m)
{
   Resize(m._nb_rows,m._nb_cols);
   _a = m._a;
   _diag = m._diag;
   _fact = m._fact;
   _rmin = m._rmin;
   _rmax = m._rmax;
   _cmin = m._cmin;
   _cmax = m._cmax;
#ifdef _OFELI_DEBUG
   std::clog << "An instance of class DMatrix is constructed.\n";
   std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
#endif
}


template<class T_>
DMatrix<T_>::~DMatrix()
{
#ifdef _OFELI_DEBUG
   std::clog << "An instance of class DMatrix is destructed.\n";
   std::clog << "File : " << __FILE__ << ", Line : " << __LINE__ << endl;
#endif
}


#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_>
void DMatrix<T_>::setMesh(const class Mesh &mesh, size_t dof, size_t type)
{
   Matrix<T_>::_init_set_mesh(mesh,dof,type);
   _length = lsize_t(_size*_size);
   _diag.resize(_size);
   _a.resize(_length);
   Clear(_a);
   _zero = T_(0);
   _fact = false;
}
#endif /* DOXYGEN_SHOULD_SKIP_THIS */


#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_>
void DMatrix<T_>::setMesh(size_t dof, const class Mesh &mesh, int code)
{
// This is just to avoid warning on unused variable
   dof  = 0;
   code = 0;
   if (mesh.getDim()==0) { }
}
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_>
void DMatrix<T_>::setMesh(size_t dof, size_t nb_eq, const class Mesh &mesh)
{
// This is just to avoid warning on unused variable
   dof = 0;
   nb_eq = 0;
   if (mesh.getDim()==0) { }
}
#endif /* DOXYGEN_SHOULD_SKIP_THIS */


template<class T_>
void DMatrix<T_>::setDiag()
{
   for (size_t i=0; i<_size; i++)
      _diag[i] = _a[_nb_cols*i+i];
}


template<class T_>
void DMatrix<T_>::Resize(size_t size)
{
   _nb_rows = _nb_cols = _size = size;
   _length = lsize_t(_nb_rows*_nb_cols);
   _a.resize(_length);
   Clear(_a);
   _rmin = 1;
   _rmax = _nb_rows;
   _cmin = 1;
   _cmax = _nb_cols;
}


template<class T_>
void DMatrix<T_>::Resize(size_t nr, size_t nc)
{
   _nb_rows = nr;
   _nb_cols = nc;
   _size = 0;
   if (_nb_rows==_nb_cols)
      _size = _nb_rows;
   _length = lsize_t(_nb_rows*_nb_cols);
   _a.resize(_length);
   Clear(_a);
   _rmin = 1;
   _rmax = _nb_rows;
   _cmin = 1;
   _cmax = _nb_cols;
}


template<class T_>
Vect<T_>DMatrix<T_>::getColumn(size_t j) const
{
   Vect<T_> v(_nb_rows);
   for (size_t i=0; i<_nb_rows; i++)
      v[i] = _a[_nb_cols*i+j-1];
   return v;
}


template<class T_>
Vect<T_> DMatrix<T_>::getRow(size_t i) const
{
   Vect<T_> v(_nb_cols);
   for (size_t j=0; j<_nb_cols; j++)
      v[j] = _a[_nb_cols*(i-1)+j];
   return v;
}


template<class T_>
void DMatrix<T_>::Set(size_t i, size_t j, const T_ &x)
{
   _a[_nb_cols*(i-1)+j-1] = x;
}


template<class T_>
void DMatrix<T_>::MultAdd(T_ a, const Vect<T_> &x, Vect<T_> &y) const
{
  for (size_t i=0; i<_nb_rows; i++)
     for (size_t j=0; j<_nb_cols; j++)
        y[i] += a * _a[_nb_cols*i+j] * x[j];
}


template<class T_>
void DMatrix<T_>::MultAdd(const Vect<T_> &x, Vect<T_> &y) const
{
   for (size_t i=0; i<_nb_rows; i++)
      for (size_t j=0; j<_nb_cols; j++)
         y[i] += _a[_nb_cols*i+j] * x[j];
}


template<class T_>
void DMatrix<T_>::Mult(const Vect<T_> &x, Vect<T_> &y) const
{
   Clear(y);
   MultAdd(x,y);
}


template<class T_>
void DMatrix<T_>::TMult(const Vect<T_> &x, Vect<T_> &y) const
{
   for (size_t i=0; i<_nb_rows; i++)
      for (size_t j=0; j<_nb_cols; j++)
         y[i] += _a[_nb_cols*(j-1)+i-1] * x[j];
}


template<class T_>
void DMatrix<T_>::Add(size_t i, size_t j, const T_ &x)
{
   _a[_nb_cols*(i-1)+j-1] += x; 
}


template<class T_>
T_ DMatrix<T_>::operator()(size_t i, size_t j) const
{
#ifdef _BOUNDS
   assert(i>0);
   assert(i<=_nb_rows);
   assert(j>0);
   assert(j<=_nb_cols);
#endif
   return _a[_nb_cols*(i-1)+j-1];
}


template<class T_>
T_ & DMatrix<T_>::operator()(size_t i, size_t j)
{
#ifdef _BOUNDS
   assert(i>0);
   assert(i<=_nb_rows);
   assert(j>0);
   assert(j<=_nb_cols);
#endif
   return _a[_nb_cols*(i-1)+j-1];
}


template<class T_>
const T_ DMatrix<T_>::operator()(size_t i) const
{
   return _a[i-1];
}


template<class T_>
int DMatrix<T_>::Factor()
{
   register size_t j, k;
   if (_size==0)
      _e.Message(__FILE__,__LINE__,62);
   for (size_t i=1; i<_size; i++) {
      for (j=1; j<=i; j++) {
         if (Abs(_a[_nb_rows*(j-1)+j-1]) < OFELI_TOLERANCE)
            _e.Message(__FILE__,__LINE__,61,int(j),0,_a[_nb_rows*(j-1)+j-1]);
         _a[_nb_rows*i+j-1] /= _a[_nb_rows*(j-1)+j-1];
         for (k=0; k<j; k++)
            _a[_nb_rows*i+j] -= _a[_nb_rows*i+k]*_a[_nb_rows*k+j];
      }
      for (j=i+1; j<_size; j++)
         for (k=0; k<i; k++)
            _a[_nb_rows*i+j] -= _a[_nb_rows*i+k]*_a[_nb_rows*k+j];
   }
   _fact = true;
   return 0;
}


template<class T_>
int DMatrix<T_>::Solve(Vect<T_> &b)
{
   if (!_fact)
      _e.Message(__FILE__,__LINE__,72);
   if (_size==0)
      _e.Message(__FILE__,__LINE__,73);
   for (size_t i=0; i<_size; i++) {
      T_ s = 0;
      for (size_t j=0; j<i; j++)
         s += _a[_nb_cols*i+j] * b[j];
      b[i] -= s;
   }
   for (int ii=int(_size)-1; ii>-1; ii--) {
      if (Abs(_a[_nb_cols*ii+ii]) < OFELI_TOLERANCE)
         _e.Message(__FILE__,__LINE__,71,ii+1,0,_a[_nb_cols*ii+ii]);
      b[ii] /= _a[_nb_cols*ii+ii];
      for (size_t j=0; j<size_t(ii); j++)
         b[j] -= b[ii] * _a[_nb_cols*j+ii];
   }
   return 0;
}


template<class T_>
int DMatrix<T_>::Solve(const Vect<T_> &b, Vect<T_> &x)
{
   x = b;
   return Solve(x);
}


template<class T_>
DMatrix<T_> & DMatrix<T_>::operator=(DMatrix<T_> &m)
{
   _a = m._a;
   return *this;
}


template<class T_>
DMatrix<T_> & DMatrix<T_>::operator+=(const DMatrix<T_> &m)
{
   for (size_t i=0; i<_length; i++)
      _a[i] += m._a[i];
   return *this;
}


template<class T_>
DMatrix<T_> & DMatrix<T_>::operator-=(const DMatrix<T_> &m)
{
   for (size_t i=0; i<_length; i++)
      _a[i] -= m._a[i];
   return *this;
}


template<class T_>
DMatrix<T_> & DMatrix<T_>::operator=(const T_ &x)
{
   for (size_t i=0; i<_length; i++)
      _a[i] = x;
   return *this;
}


template<class T_>
DMatrix<T_> & DMatrix<T_>::operator*=(const T_ &x)
{
   for (size_t i=0; i<_length; i++)
      _a[i] *= x;
   return *this;
}


template<class T_>
DMatrix<T_> & DMatrix<T_>::operator+=(const T_ &x)
{
   for (size_t i=0; i<_length; i++)
      _a[i] += x;
   return *this;
}


template<class T_>
DMatrix<T_> & DMatrix<T_>::operator-=(const T_ &x)
{
   for (size_t i=0; i<_length; i++)
      _a[i] -= x;
   return *this;
}


template<class T_>
size_t DMatrix<T_>::getColInd(size_t i) const
{
#ifdef _BOUNDS_DEBUG
   assert(i<=_length && i>0);
#endif
   return _ch[i-1];
}


#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_>
size_t DMatrix<T_>::getRowPtr(size_t i) const
{
#ifdef _BOUNDS_DEBUG
   assert(i<=_nb_rows+1 && i>0);
#endif
   return _ch[i-1];
}
#endif /* DOXYGEN_SHOULD_SKIP_THIS */


template<class T_>
T_ * DMatrix<T_>::getArray() const
{
   return _a;
}


template<class T_>
T_ DMatrix<T_>::getEntry(size_t i, size_t j) const 
{
   return _a[_nb_cols*(i-1)+j-1]; 
}


template<class T_>
void DMatrix<T_>::setPrintView(size_t rmin, size_t rmax, size_t cmin, size_t cmax)
{
   _rmin = rmin;
   _rmax = rmax;
   if (_rmin==0)
      _rmin = 1;
   if (_rmax>_nb_rows)
      _rmax = _nb_rows;
   _cmin = cmin;
   _cmax = cmax;
   if (_cmin==0)
      _cmin = 1;
   if (_cmax>_nb_cols)
      _cmax = _nb_cols;
}


//                           ASSOCIATED  FUNCTIONS                           //


template<class T_>
ostream& operator<<(ostream& s, const DMatrix<T_> &a)
{
   size_t rmin, rmax, cmin, cmax;
   a.getPrintView(rmin,rmax,cmin,cmax);
   s.setf(ios::right|ios::scientific);
   s << endl;
   for (size_t i=rmin-1; i<rmax; i++) {
      s << "\nRow  " << setw(6) << i+1 << endl;
      for (size_t j=cmin-1; j<cmax; j++)
         s << "  " << setprecision(8) << std::setfill(' ') << setw(18) << a(i+1,j+1);
      s << endl;
   }
   return s;
}

} /* namespace OFELI */

#endif

Copyright © 1998-2007 Rachid Touzani