DVector.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 DVector_H
#define DVector_H

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

using namespace std;


class DVector
{
 public:
          friend class DMatrix;

 private:
  int n_1;


  double* pArray;


  int lda;

public:
  int dim1(){return n_1;}

  int size(){return n_1;}

  double* array(){return pArray;}

  void allocate(int i)
    {
      if( n_1 != i )
        {
          deallocate();
          
          n_1 = i;
          
          if(n_1 >= 1) 
            {
              pArray = new 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(dnrm2,DNRM2)(&n_1,pArray,&lda);
    }
 
    double AbsoluteSum()
    {
      return F77_FUNC(dasum,DASUM)(&n_1,pArray,&lda);
    }
   
    double Sum()
    {
      double sum=0;
      for(int i=0;i<n_1;i++) sum+=pArray[i];
      return sum;
    }


    int iMax()
    {
      return F77_FUNC(idamax,IDAMAX)(&n_1,pArray,&lda);
    }


  void operator=(const DVector & rhs)
  {
   if(n_1 != rhs.n_1) allocate(rhs.n_1);
   F77_FUNC(dcopy,DCOPY)(&n_1,rhs.pArray,&lda,pArray,&lda);
  }
  
  void operator=(const double C)
  {
    for(int i=0; i<n_1; i++)
      {
        pArray[i] = C;
      }
  }
  

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


  DVector operator*(double rhs)
    {
      DVector  A(n_1);
      double zero=0.0;
      F77_FUNC(dscal,DSCAL)(&n_1,&zero,A.pArray,&lda);
      F77_FUNC(daxpy,DAXPY)(&n_1,&rhs,pArray,&lda,A.pArray,&lda);
      return A;
    }


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


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


  void operator*=(double C)
    {
    F77_FUNC(dscal,DSCAL)(&n_1,&C,pArray,&lda);
    }


  void operator/=(double C)
    {
    double inv=1.0/C;
    F77_FUNC(dscal,DSCAL)(&n_1,&inv,pArray,&lda);
    }


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


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


  DVector(const DVector & rhs)
    {
      n_1 = 0;
      pArray = 0;
      lda = 1;
      allocate(rhs.n_1);
      F77_FUNC(dcopy,DCOPY)(&n_1,rhs.pArray,&lda,pArray,&lda);
    }

  DVector(int dim,double *rhs)
    {
      n_1 = 0;
      pArray = 0;
      lda = 1;
      allocate(dim);
      n_1=dim;
      F77_FUNC(dcopy,DCOPY)(&n_1,rhs,&lda,pArray,&lda);
    }

   
  void set(int ndim,double *rhs)
    {
#ifdef DEBUG
       if( n_1 != ndim )
        {
          cerr << "ERROR: incorrect sizes for two Vector's in "
               << " DVector.set(n_1,c_array) " << endl;
          exit(1);
        }
#endif
      F77_FUNC(dcopy,DCOPY)(&n_1,rhs,&lda,pArray,&lda);
    }
  
  void get(int dim,double *rhs)
    {
#ifdef DEBUG
       if( n_1 != dim )
        {
          cerr << "ERROR: incorrect sizes for two Vector's in "
               << " Vector.get(n_1,c_array) " << endl;
          exit(1);
        }
#endif
      F77_FUNC(dcopy,DCOPY)(&n_1,pArray,&lda,rhs,&lda);
    }
  
  ~DVector(){deallocate();}

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

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

#endif

---