ZVector.h

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2005 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 ZVector_H
#define ZVector_H

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

using namespace std;


class ZVector
{
 public:
          friend class ZMatrix;

 private:
  int n_1;


  complex<double>* pArray;


  int lda;

public:
  int dim1(){return n_1;}

  int size(){return n_1;}

  complex<double>* array(){return pArray;}

  void allocate(int i)
    {
      if( n_1 != i )
        {
          deallocate();
          
          n_1 = i;
          
          if(n_1 >= 1) 
            {
              pArray = new complex<double>[n_1];
            }
          else 
            {
              n_1 = 0;
              pArray = 0;
            }
        }

      lda = 1;
    }


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

          n_1 = 0;
        }
    }

   
    double Norm()
    {
      return F77_FUNC(dznrm2,DZNRM2)(&n_1,(double*)pArray,&lda);
    }
 
    double AbsoluteSum()
    {
      return F77_FUNC(dzasum,DZASUM)(&n_1,(double*)pArray,&lda);
    }
   
    complex<double> Sum()
    {
      complex<double> sum=dcmplx(0.0,0.0);
      for(int i=0;i<n_1;i++) sum=sum+pArray[i];
      return sum;
    }


    int iMax()
    {
      return F77_FUNC(izamax,IZAMAX)(&n_1,(double*)pArray,&lda);
    }


  void operator=(const ZVector & rhs)
  {
   if(n_1 != rhs.n_1) allocate(rhs.n_1);
   F77_FUNC(zcopy,ZCOPY)(&n_1,(double*)rhs.pArray,&lda,(double*)pArray,&lda);
  }
  
  void operator=(const double C)
  {
    for(int i=0; i<n_1; i++)
      {
        pArray[i] = dcmplx(C,0);
      }
  }
  
  void operator=(const complex<double> C)
  {
    for(int i=0; i<n_1; i++)
      {
        pArray[i] = C;
      }
  }
  
  complex<double> zdotu(const ZVector & rhs)
    {
#ifdef DEBUG
      if( n_1 != rhs.n_1 )
        {
          cerr << "ERROR: incorrect sizes for two ZVector's in ZVector*"
               << "ZVector" << endl;
          exit(1);
        }
      if( n_1 <= 0 )
        {
          cerr << "ERROR: ZVector of size 0 used in ZVector*ZVector" << endl;
          exit(1);
        }
#endif
      return  F77_FUNC(zdotu,ZDOTU)(&n_1,(double*)pArray,&lda,(double*)rhs.pArray,&lda);
    }

    complex<double> zdotc(const ZVector & rhs)
    {
#ifdef DEBUG
      if( n_1 != rhs.n_1 )
        {
          cerr << "ERROR: incorrect sizes for two ZVector's in ZVector*"
               << "ZVector" << endl;
          exit(1);
        }
      if( n_1 <= 0 )
        {
          cerr << "ERROR: ZVector of size 0 used in ZVector*ZVector" << endl;
          exit(1);
        }
#endif
      return  F77_FUNC(zdotc,ZDOTC)(&n_1,(double*)pArray,&lda,(double*)rhs.pArray,&lda);
    }

  ZVector operator+( const ZVector & rhs)
    {
#ifdef DEBUG
      if( n_1 != rhs.n_1 )
        {
          cerr << "ERROR: incorrect sizes for two ZVector's in ZVector+"
               << "ZVector" << endl;
          exit(1);
        }
      if( n_1 <= 0 )
        {
          cerr << "ERROR: ZVector of size 0 used in ZVector+ZVector" << endl;
          exit(1);
        }
#endif
      ZVector  A(n_1);
      double d_one=1.0;
      F77_FUNC(zcopy,ZCOPY)(&n_1,(double*)pArray,&lda,(double*)A.pArray,&lda);
      F77_FUNC(zaxpy,ZAXPY)(&n_1,&d_one,(double*)rhs.pArray,&lda,(double*)A.pArray,&lda); 
      return A;
    }


  ZVector operator-( const ZVector & rhs)
    {
#ifdef DEBUG
      if( n_1 != rhs.n_1 )
        {
          cerr << "ERROR: incorrect sizes for two ZVector's in ZVector-"
               << "ZVector" << endl;
          exit(1);
        }
      if( n_1 <= 0 )
        {
          cerr << "ERROR: ZVector of size 0 used in ZVector-ZVector" << endl;
          exit(1);
        }
#endif
      ZVector  A(n_1);
      complex<double> d_m_one(-1.0,0.0);
      F77_FUNC(zcopy,ZCOPY)(&n_1,(double*)pArray,&lda,(double*)A.pArray,&lda);
      F77_FUNC(zaxpy,ZAXPY)(&n_1,(double*)&d_m_one,(double*)rhs.pArray,&lda,(double*)A.pArray,&lda); 
      return A;
    }


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


  void operator/=(complex<double> C)
    {
    complex<double> inv=1.0/C;
    F77_FUNC(zscal,ZSCAL)(&n_1,(double*)&inv,(double*)pArray,&lda);
    }


  ZVector(){pArray = 0; n_1 = 0; lda=1;}


  ZVector(int i){pArray = 0; n_1 = 0; allocate(i); lda=1; }


  ZVector(const ZVector & rhs)
    {
      n_1 = 0;
      pArray = 0;
      lda = 1;
      allocate(rhs.n_1);
      F77_FUNC(zcopy,ZCOPY)(&n_1,(double*)rhs.pArray,&lda,(double*)pArray,&lda);
    }

  ~ZVector(){deallocate();}

  complex<double>& operator()(int i){return pArray[i];}

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

#endif

---