SVector.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 SVector_H
#define SVector_H

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

using namespace std;


class SVector
{
 public:
          friend class SMatrix;

 private:
  int n_1;


  float* pArray;


  int lda;

public:
  int dim1(){return n_1;}

  int size(){return n_1;}

  float* array(){return pArray;}

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


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

          n_1 = 0;
        }
    }

   
    float Norm()
    {
      return F77_FUNC(snrm2,SNRM2)(&n_1,pArray,&lda);
    }
 
    float AbsoluteSum()
    {
      return F77_FUNC(sasum,SASUM)(&n_1,pArray,&lda);
    }
   
    float Sum()
    {
      float sum=0;
      for(int i=0;i<n_1;i++) sum+=pArray[i];
      return sum;
    }


    int iMax()
    {
      return F77_FUNC(isamax,ISAMAX)(&n_1,pArray,&lda);
    }


  void operator=(const SVector & rhs)
  {
   if(n_1 != rhs.n_1) allocate(rhs.n_1);
   F77_FUNC(scopy,SCOPY)(&n_1,rhs.pArray,&lda,pArray,&lda);
  }
  
  void operator=(const float C)
  {
    for(int i=0; i<n_1; i++)
      {
        pArray[i] = C;
      }
  }
  

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


  SVector operator*(float rhs)
    {
      SVector  A(n_1);
      float zero=0.0;
      F77_FUNC(sscal,SSCAL)(&n_1,&zero,A.pArray,&lda);
      F77_FUNC(saxpy,SAXPY)(&n_1,&rhs,pArray,&lda,A.pArray,&lda);
      return A;
    }


  SVector operator+( const SVector & rhs)
    {
#ifdef DEBUG
      if( n_1 != rhs.n_1 )
        {
          cerr << "ERROR: incorrect sizes for two SVector's in SVector+"
               << "SVector" << endl;
          exit(1);
        }
      if( n_1 <= 0 )
        {
          cerr << "ERROR: SVector of size 0 used in SVector+SVector" << endl;
          exit(1);
        }
#endif
      SVector  A(n_1);
      float d_one=1.0;
      F77_FUNC(scopy,SCOPY)(&n_1,pArray,&lda,A.pArray,&lda);
      F77_FUNC(saxpy,SAXPY)(&n_1,&d_one,rhs.pArray,&lda,A.pArray,&lda);
      return A;
    }


  SVector operator-( const SVector & rhs)
    {
#ifdef DEBUG
      if( n_1 != rhs.n_1 )
        {
          cerr << "ERROR: incorrect sizes for two SVector's in SVector-"
               << "SVector" << endl;
          exit(1);
        }
      if( n_1 <= 0 )
        {
          cerr << "ERROR: SVector of size 0 used in SVector-SVector" << endl;
          exit(1);
        }
#endif
      SVector  A(n_1);
      float d_m_one=-1.0;
      F77_FUNC(scopy,SCOPY)(&n_1,pArray,&lda,A.pArray,&lda);
      F77_FUNC(saxpy,SAXPY)(&n_1,&d_m_one,rhs.pArray,&lda,A.pArray,&lda);
      return A;
    }


  void operator*=(float C)
    {
    F77_FUNC(sscal,SSCAL)(&n_1,&C,pArray,&lda);
    }


  void operator/=(float C)
    {
    float inv=1.0/C;
    F77_FUNC(sscal,SSCAL)(&n_1,&inv,pArray,&lda);
    }


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


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


  SVector(const SVector & rhs)
    {
      n_1 = 0;
      pArray = 0;
      lda = 1;
      allocate(rhs.n_1);
      F77_FUNC(scopy,SCOPY)(&n_1,rhs.pArray,&lda,pArray,&lda);
    }

  SVector(int dim,float *rhs)
    {
      n_1 = 0;
      pArray = 0;
      lda = 1;
      allocate(dim);
      n_1=dim;
      F77_FUNC(scopy,SCOPY)(&n_1,rhs,&lda,pArray,&lda);
    }

   
  void set(int ndim,float *rhs)
    {
#ifdef DEBUG
       if( n_1 != ndim )
        {
          cerr << "ERROR: incorrect sizes for two Vector's in "
               << " SVector.set(n_1,c_array) " << endl;
          exit(1);
        }
#endif
      F77_FUNC(scopy,SCOPY)(&n_1,rhs,&lda,pArray,&lda);
    }
  
  void get(int dim,float *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(scopy,SCOPY)(&n_1,pArray,&lda,rhs,&lda);
    }
  
  ~SVector(){deallocate();}

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

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

#endif

---