DMatrix.h

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2004 by Claudio Attaccalite                             *
 *   claudio@freescience.info                                              *
 *                                                                         *
 *   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; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   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.             *
 ***************************************************************************/


#ifndef DMatrix_H
#define DMatrix_H
#if HAVE_CONFIG
#include <config.h>
#endif

#include <iostream>
#include "DVector.h"
#include "lapack_wrapper.h"
#include "blas_wrapper.h"


using namespace std;


class DMatrix
{
  private:
  int n_1;


  int n_2;

  double* pArray;

  double* tmpArray;

  double temp;

  int one;

public:
  int dim1(){return n_1;}

  int dim2(){return n_2;}

  int size(){return n_1*n_2;}

  double* array(){return pArray;}

 void transpose()
 {
 if(n_1 !=n_2)
 {
         cerr << " Transpose only for square matrix ! " << endl;
         exit(1);
 }
 for(int i=0;i<n_1;i++)
        for(int j=0;j<i;j++)
        {
                temp=pArray[j+i*n_1];
                pArray[j+i*n_1]=pArray[i+j*n_1];
                pArray[i+j*n_1]=temp;
        }
 }


  void allocate(int i, int j)
    {
      if( n_1 != i || n_2 != j )
        {
          deallocate();
          
          n_1 = i;
          n_2 = j;
          
          if(n_1 >= 1 && n_2 >= 1)
              pArray = new double[n_1*n_2];
          else 
              pArray = 0;
        }
    }


  void deallocate()
    {
      if(n_2>0 && n_1>0)
      {
       delete [] pArray;
      }
       pArray = 0;

      n_1 = 0;
      n_2 = 0;
    }


  void operator=(const DMatrix & rhs)
    {
     int tot=n_1*n_2;
#ifdef DEBUG
     if(n_1 != rhs.n_1 || n_2 != rhs.n_2) allocate(rhs.n_1,rhs.n_2);
#endif
     
         F77_FUNC(dcopy,DCOPY)(&tot,rhs.pArray,&one,pArray,&one);
    }


  void operator=(const double C)
    {
      for(int i=0; i<n_1*n_2;i++) pArray[i] = C;
    }

  void operator*=(double C)
    {
        int tot=n_1*n_2;
        F77_FUNC(dscal,DSCAL)(&tot,&C,pArray,&one);
    }

  void operator/=(const double C)
    {
        int tot=n_1*n_2;
        double inv=1.0/C;
        F77_FUNC(dscal,DSCAL)(&tot,&inv,pArray,&one);
    }


  DMatrix(){pArray = 0; n_1 = 0; n_2 = 0; one = 1;}


  DMatrix(int i, int j) {pArray = 0; n_1 = 0; n_2 = 0; one=1; allocate(i,j);}


  DMatrix( const DMatrix & rhs)
    {
      n_1 = 0;
      n_2 = 0;
      pArray = 0;
      allocate(rhs.n_1,rhs.n_2);
      int tot=n_1*n_2;
      one=1;
      F77_FUNC(dcopy,DCOPY)(&tot,rhs.pArray,&one,pArray,&one);
    }


  ~DMatrix(){deallocate();}

  double SumAll()
  {
   double sum=0;
   for(int i=0;i<n_1*n_2;i++) sum+=pArray[i];
   return sum;
  }
 
  double Trace()
  {
#ifdef DEBUG
   if(n_1 !=n_2)
   {
         cerr << " Trace only for Squared DMatrix ! " << endl;
         exit(1);
   } 
#endif
   double sum=0;
   for(int i=0;i<n_1;i++) sum+=pArray[i*n_1+i];
   return sum;
  }
  
  void setrow(int irow,DVector &vector)
  {
#ifdef DEBUG
          if(vector.n_1!=n_2 || irow>=n_1)
          { 
                  cerr << " ERROR! dimension array copying row! " << endl;
                  exit(1);
          }
#endif
   F77_FUNC(dcopy,DCOPY)(&n_2,vector.pArray,&one,pArray+irow,&n_1);
  }

  void setrow(int irow,double *vector)
  {
   F77_FUNC(dcopy,DCOPY)(&n_2,vector,&one,pArray+irow,&n_1);
  }

  void setrow(int irow,double *vector,int lda)
  {
   F77_FUNC(dcopy,DCOPY)(&n_2,vector,&lda,pArray+irow,&n_1);
  }
  
  void setcol(int icol, const DVector & vector)
  {
#ifdef DEBUG
          if(vector.n_1!=n_1 || icol>=n_2)
          { 
                  cerr << " ERROR! dimension array copying coloumn! " << endl;
                  exit(1);
          }
#endif
   F77_FUNC(dcopy,DCOPY)(&n_1,vector.pArray,&one,pArray+icol*n_1,&one);
  }

  void setcol(int icol,double *vector)
  {
#ifdef DEBUG
          if(icol>=n_2)
          { 
                  cerr << " ERROR! dimension array copying coloumn! " << endl;
                  exit(1);
          }
#endif
   F77_FUNC(dcopy,DCOPY)(&n_1,vector,&one,pArray+icol*n_1,&one);
  }
 
  void row_scale(int irow, double C)
  {
#ifdef DEBUG
          if(irow>=n_1)
          { 
                  cerr << " ERROR! dimension array in row_scale! " << endl;
                  exit(1);
          }
#endif
   F77_FUNC(dscal,DSCAL)(&n_2,&C,pArray+irow,&n_1);
  }



  
  double rowdot(int irow, const DVector & vector)
  {
#ifdef DEBUG
          if(vector.n_1!=n_2 || irow>=n_1)
          { 
                  cerr << " ERROR! dimension array in row dot! " << endl;
                  exit(1);
          }
#endif
   return F77_FUNC(ddot,DDOT)(&n_2,pArray+irow,&n_1,vector.pArray,&one);
  }

  double coldot(int icol, const DVector &vector)
  {
#ifdef DEBUG
          if(vector.n_1!=n_1 || icol>=n_2)
          { 
                  cerr << " ERROR! dimension array in coloumn dot! " << endl;
                  exit(1);
          }
#endif
   return F77_FUNC(ddot,DDOT)(&n_1,pArray+icol*n_1,&one,vector.pArray,&one);
  }
  
  void col_scale(int icol, double C)
  {
#ifdef DEBUG
          if(icol>=n_2)
          { 
                  cerr << " ERROR! dimension array in col_rescale! " << endl;
                  exit(1);
          }
#endif
   F77_FUNC(dscal,DSCAL)(&n_1,&C,pArray+icol*n_1,&one);
  }


  double rowdot(int irow, double *vector)
  {
   return F77_FUNC(ddot,DDOT)(&n_2,pArray+irow,&n_1,vector,&one);
  }

  double coldot(int icol, double *vector)
  {
   return F77_FUNC(ddot,DDOT)(&n_1,pArray+icol*n_1,&one,vector,&one);
  }



  double *row_pointer(int irow)
  {
#ifdef DEBUG
          if(irow>=n_1)
          { 
                  cerr << " ERROR! dimension array getting a row! " << endl;
                  exit(1);
          }
#endif
   return(pArray+irow);
  }
  
  DVector row(int irow)
  {
#ifdef DEBUG
          if(irow>=n_1)
          { 
                  cerr << " ERROR! dimension array getting a row! " << endl;
                  exit(1);
          }
#endif
   DVector TEMP(n_2);
   F77_FUNC(dcopy,DCOPY)(&n_2,pArray+irow,&n_1,TEMP.pArray,&one);
   return TEMP;
  }




  double *col_pointer(int icol)
  {
#ifdef DEBUG
          if(icol>=n_2)
          { 
                  cerr << " ERROR! dimension array getting a coloumn! " << endl;
                  exit(1);
          }
#endif
   return(pArray+icol*n_1);
  }


  
  DVector col(int icol)
  {
#ifdef DEBUG
          if(icol>=n_2)
          { 
                  cerr << " ERROR! dimension array getting a coloumn! " << endl;
                  exit(1);
          }
#endif
   DVector TEMP(n_1);
   F77_FUNC(dcopy,DCOPY)(&n_1,pArray+icol*n_1,&one,TEMP.pArray,&one);
   return TEMP;
  }

  double SumRow(const int irow)
  {
#ifdef DEBUG
          if(irow>=n_1)
          { 
                  cerr << " ERROR! dimension array getting a row! " << endl;
                  exit(1);
          }
#endif
   double sum=0;
   for(int i=0;i<n_2;i++) sum+=pArray[irow+i*n_1];
   return sum;
  }


  double SumCol(const int icol)
  {
#ifdef DEBUG
          if(icol>=n_2)
          { 
                  cerr << " ERROR! dimension array getting a row! " << endl;
                  exit(1);
          }
#endif
   double sum=0;
   for(int i=0;i<n_1;i++) sum+=pArray[i+icol*n_1];
   return sum;
  }





  DMatrix operator+( const DMatrix & rhs)
    {
#ifdef DEBUG
      if( n_1 != rhs.n_1 || n_2 !=rhs.n_2 )
        {
          cerr << "ERROR: incorrect sizes for two DMatrix's in DMatrix+"
               << "DMatrix" << endl;
          exit(1);
        }
      if( n_1 <= 0 || n_2 <=0 )
        {
          cerr << "ERROR: DMatrix of size 0 used in DMatrix+DMatrix" << endl;
          exit(1);
        }
#endif
      int tot=n_1*n_2;
      double d_one=1.0;
      DMatrix  A(n_1,n_2);
      F77_FUNC(dcopy,DCOPY)(&tot,pArray,&one,A.pArray,&one);
      F77_FUNC(daxpy,DAXPY)(&tot,&d_one,rhs.pArray,&one,A.pArray,&one); 
      return A;
    }


  DMatrix operator-( const DMatrix & rhs)
    {
#ifdef DEBUG
      if( n_1 != rhs.n_1 || n_2 !=rhs.n_2 )
        {
          cerr << "ERROR: incorrect sizes for two DMatrix's in DMatrix-"
               << "DMatrix" << endl;
          exit(1);
        }
      if( n_1 <= 0 || n_2 <=0)
        {
          cerr << "ERROR: DMatrix of size 0 used in DMatrix-DMatrix" << endl;
          exit(1);
        }
#endif
      DMatrix A(n_1,n_2);
      int tot=n_1*n_2;
      double d_m_one=-1.0;
      F77_FUNC(dcopy,DCOPY)(&tot,pArray,&one,A.pArray,&one);
      F77_FUNC(daxpy,DAXPY)(&tot,&d_m_one,rhs.pArray,&one,A.pArray,&one); 
      return A;
    }

  
  double & operator()(int i,int j)
  {
#ifdef DEBUG
          if(i>=n_1 || j>=n_2)
          { 
                  cerr << " ERROR! index out of range! " 
                       << "("<<i<<","<<j<<")" <<  endl;
                  exit(1);
          }
#endif
          return pArray[i+j*n_1];
  }


  friend ostream& operator<<(ostream & strm, const DMatrix & rhs)
    {
      for(int i=0; i<rhs.n_1;i++)
        {
          for(int j=0; j<rhs.n_2; j++)
            {
              strm << rhs.pArray[i+j*rhs.n_1] << "\t";
            }
          strm << endl;
        }
      return strm;
    }

};

#endif

---