ZMatrix.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 ZMatrix_H
#define ZMatrix_H
#if HAVE_CONFIG
#include <config.h>
#endif

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


using namespace std;


class ZMatrix
{
  private:
  int n_1;


  int n_2;

  complex<double>* pArray;

  complex<double>* tmpArray;

  complex<double> temp;

  int one;

public:
  int dim1(){return n_1;}

  int dim2(){return n_2;}

  int size(){return n_1*n_2;}

  complex<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 complex<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 ZMatrix & 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(zcopy,ZCOPY)(&tot,(double*)rhs.pArray,&one,(double*)pArray,&one);
    }


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

  void operator*=(complex<double> C)
    {
        int tot=n_1*n_2;
        F77_FUNC(zscal,ZSCAL)(&tot,(double*)&C,(double*)pArray,&one);
    }

  void operator/=(const complex<double> C)
    {
        int tot=n_1*n_2;
        complex<double> inv=dcmplx(1.0,0.0)/C;
        F77_FUNC(zscal,ZSCAL)(&tot,(double*)&inv,(double*)pArray,&one);
    }


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


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


  ZMatrix( const ZMatrix & 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(zcopy,ZCOPY)(&tot,(double*)rhs.pArray,&one,(double*)pArray,&one);
    }


  ~ZMatrix(){deallocate();}

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

---