Array1D.h

Go to the documentation of this file.

//            QMcBeaver
//
//         Constructed by
//
//     Michael Todd Feldmann
//              and
//   David Randall "Chip" Kent IV
//
// Copyright 2000.  All rights reserved.
//
// drkent@users.sourceforge.net mtfeldmann@users.sourceforge.net

#ifndef Array1D_H
#define Array1D_H

#include <iostream>
#include <iomanip>
#include <math.h>
#include "cppblas.h"
#include <assert.h>

using namespace std;

template <class T> class Array1D
  {
  private:

    int n_1;

    T* pArray;

  public:

    int dim1() const
      {
        return n_1;
      }

    int size() const
      {
        return n_1;
      }

    T* array()
    {
      return pArray;
    }

    void allocate(int i)
    {
#ifdef QMC_DEBUG
      /*
        //There are a number of places in the code that allocate
        //negative numbers, e.g. Polynomial::initialize,
        //i don't want to change them, so i'm
        //commenting out this warning.
      if ( i < 0 )
        {
          cerr << "Warning: invalid dimension in Array1D::allocate\n"
               << " i = " << i << endl;
        }
      */
#endif
      if ( n_1 != i )
        {
          deallocate();
          n_1 = i;
          if (n_1 >= 1)
            {
              pArray = new T[n_1];
            }
          else
            {
              n_1 = 0;
              pArray = 0;
            }
        }
    }

    void resize(int i)
    {
      allocate(i);
    }

    void clear()
    {
      deallocate();
    }

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

#ifdef USEBLAS
    T operator*( const Array1D<double> & rhs)
    {
#ifdef QMC_DEBUG
      assert(n_1 > 0);
      assert(pArray);
      assert(rhs.pArray);
#endif
#ifdef USE_CBLAS
      return cblas_ddot(n_1, pArray, 1,rhs.pArray, 1);
#else
      int inc = 1;
      return FORTRAN_FUNC(ddot)(&n_1, pArray, &inc, rhs.pArray, &inc);
#endif
    }

    T operator*( const Array1D<float> & rhs)
    {
#ifdef QMC_DEBUG
      assert(n_1 > 0);
      assert(pArray);
      assert(rhs.pArray);
#endif
#ifdef USE_CBLAS
      return cblas_sdot(n_1, pArray, 1,rhs.pArray, 1);
#else
      int inc = 1;
      return FORTRAN_FUNC(sdot)(&n_1, pArray, &inc, rhs.pArray, &inc);
#endif
    }

    //this method has not been tested with ATLAS
    Array1D operator+( const Array1D & rhs)
    {
      Array1D <T> A(rhs);
#ifdef USE_CBLAS
      cblas_daxpy(n_1, 1.0, pArray,1, A.pArray, 1);
#else
      int inc = 1;
      double a = 1.0;
      FORTRAN_FUNC(daxpy)(&n_1, &a, pArray, &inc, A.pArray, &inc);
#endif
      return A;
    }

    //this method has not been tested with ATLAS
    void operator+=( const Array1D & rhs)
    {
#ifdef USE_CBLAS
      cblas_daxpy(n_1, 1.0, rhs.pArray,1, pArray, 1);
#else
      int inc = 1;
      double a = 1.0;
      FORTRAN_FUNC(daxpy)(&n_1, &a, rhs.pArray, &inc, pArray, &inc);
#endif
    }

#else

    T operator*( const Array1D & rhs)
    {
#ifdef QMC_DEBUG
      if ( n_1 != rhs.n_1 )
        {
          cerr << "ERROR: incorrect sizes for two Array1D's in Array1D*"
          << "Array1D" << endl;
          assert(0);
        }
      if ( n_1 <= 0 )
        {
          cerr << "ERROR: Array1D of size 0 used in Array1D*Array1D" << endl;
          assert(0);
        }
      assert(pArray);
      assert(rhs.pArray);
#endif
      T temp = 0.0;
      for ( int i=0; i<n_1; i++ )
        {
          temp += pArray[i] * rhs.pArray[i];
        }
      return temp;
    }

    Array1D operator+( const Array1D & rhs)
    {
#ifdef QMC_DEBUG
      if ( n_1 != rhs.n_1 )
        {
          cerr << "ERROR: incorrect sizes for two Array1D's in Array1D+"
          << "Array1D" << endl;
          exit(1);
        }
      if ( n_1 <= 0 )
        {
          cerr << "ERROR: Array1D of size 0 used in Array1D+Array1D" << endl;
          exit(1);
        }
      assert(pArray);
      assert(rhs.pArray);
#endif
      Array1D <T> A(n_1);

      for ( int i=0; i<n_1; i++ )
        {
          A.pArray[i] = pArray[i] + rhs.pArray[i];
        }
      return A;
    }

    void operator+=( const Array1D<T> & rhs)
      {
        for(int i=0; i<n_1; i++)
          pArray[i] += rhs.pArray[i];
      }
    
#endif

    void operator=(const Array1D & rhs)
    {
      if (n_1 != rhs.n_1) allocate(rhs.n_1);

      for (int i=0; i<n_1;i++) pArray[i] = rhs.pArray[i];
      //memcpy(pArray, rhs.pArray, sizeof(T)*n_1);
    }

    void operator=(const T C)
    {
      if (C == 0)
        {
          memset(pArray,0,sizeof(T)*n_1);
          return;
        }
      for (int i=0; i<n_1; i++)
        pArray[i] = C;
    }

    Array1D operator*( const double rhs)
    {
      Array1D <T> A(n_1);
      for ( int i=0; i<n_1; i++ )
        {
          A.pArray[i] = rhs*pArray[i];
        }
      return A;
    }

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

    void operator/=(const T C)
    {
      T inv = 1.0/C;
      operator*=(inv);
    }

    Array1D operator-( const Array1D & rhs)
    {
      if ( n_1 != rhs.n_1 )
        {
          cerr << "ERROR: incorrect sizes for two Array1D's in Array1D-"
          << "Array1D" << endl;
          exit(1);
        }
      if ( n_1 <= 0 )
        {
          cerr << "ERROR: Array1D of size 0 used in Array1D-Array1D" << endl;
          exit(1);
        }
      Array1D <T> A(n_1);

      for ( int i=0; i<n_1; i++ )
        {
          A.pArray[i] = pArray[i] - rhs.pArray[i];
        }
      return A;
    }

    void quaternion_conjugate()
    {
      if ( n_1 != 4 )
        {
          cerr << "ERROR: incorrect size for quaternion conjugate." << endl;
          exit(1);
        }
      for (int i=1; i<4; i++)
        {
          pArray[i] *= -1;
        }
    }

    Array1D quaternion_product( const Array1D & rhs )
    {
      if ( n_1 != 4 || rhs.n_1 != 4 )
        {
          cerr << "ERROR: incorrect size for quaternion product." << endl;
          exit(1);
        }
      Array1D<T> product(4);
      product.pArray[0] = pArray[0]*rhs.pArray[0] - pArray[1]*rhs.pArray[1] -
                          pArray[2]*rhs.pArray[2] - pArray[3]*rhs.pArray[3];
      product.pArray[1] = pArray[0]*rhs.pArray[1] + pArray[1]*rhs.pArray[0] +
                          pArray[2]*rhs.pArray[3] - pArray[3]*rhs.pArray[2];
      product.pArray[2] = pArray[0]*rhs.pArray[2] - pArray[1]*rhs.pArray[3] +
                          pArray[2]*rhs.pArray[0] + pArray[3]*rhs.pArray[1];
      product.pArray[3] = pArray[0]*rhs.pArray[3] + pArray[1]*rhs.pArray[2] -
                          pArray[2]*rhs.pArray[1] + pArray[3]*rhs.pArray[0];
      return product;
    }

    void rotate(Array1D axis, double angle)
    {
#ifdef QMC_DEBUG
      if (n_1 != 3)
        {
          cerr << "ERROR: incorrect size for rotating point." << endl;
          exit(1);
        }
#endif
      double x = axis(0);
      double y = axis(1);
      double z = axis(2);
      double s  = sin(angle);
      double c  = cos(angle);
      double v  = 1.0 - c;
      double xs = x*s;
      double ys = y*s;
      double zs = z*s;
      double xv = x*v;
      double yv = y*v;
      double zv = z*v;
      double xyv = x*yv;
      double yzv = y*zv;
      double zxv = z*xv;

      double i = (x*xv+c)*pArray[0] + (xyv-zs)*pArray[1] + (zxv+ys)*pArray[2];
      double j = (xyv+zs)*pArray[0] + (y*yv+c)*pArray[1] + (yzv-xs)*pArray[2];
      double k = (zxv-ys)*pArray[0] + (yzv+xs)*pArray[1] + (z*zv+c)*pArray[2];
      
      pArray[0] = i;
      pArray[1] = j;
      pArray[2] = k;
    }

    void rotateQ(Array1D axis, double angle)
    {
      if (n_1 != 3)
        {
          cerr << "ERROR: incorrect size for rotating point." << endl;
          exit(1);
        }
      Array1D<T> q_point(4), q_axis(4), q1(4), q2(4);
      q_point.pArray[0] = 0.0;
      q_axis.pArray[0] = cos(angle/2);
      for (int i=1; i<4; i++)
        {
          q_point.pArray[i] = pArray[i-1];
          q_axis.pArray[i] = axis.pArray[i-1]*sin(angle/2);
        }
      q1 = q_axis.quaternion_product(q_point);
      q_axis.quaternion_conjugate();
      q2 = q1.quaternion_product(q_axis);
      for (int i=1; i<4; i++)
        {
          pArray[i-1] = q2.pArray[i];
        }
    }

    Array1D<T> cross(Array1D<T> & rhs)
      {
#ifdef QMC_DEBUG
        if(n_1 != 3 || rhs.dim1() != 3)
          {
            cout << "Error: cross product vectors are " << n_1 << "x" << rhs.dim1() << endl;
          }
#endif  
        Array1D<T> crossed(3);
        crossed[0] = pArray[1]*rhs[2] - pArray[2]*rhs[1];
        crossed[1] = pArray[2]*rhs[0] - pArray[0]*rhs[2];
        crossed[2] = pArray[0]*rhs[1] - pArray[1]*rhs[0];
        return crossed;
      }
    
    /*
      Returns the length of the vector.
    */
    T mag(){
      double sum = 0.0;
      for(int i=0; i<n_1; i++)
        sum += pArray[i]*pArray[i];
      return sqrt(sum);
    }

    Array1D()
    {
      pArray = 0; n_1 = 0;
    }

    Array1D(int i)
    {
      pArray = 0; n_1 = 0; allocate(i);
    }

  Array1D(int i, double shift, double range, long & iseed)
    {
      pArray = 0;
      n_1 = 0;
      allocate(i);
      
      if(iseed > 0) srand(iseed);
      for(int idx=0; idx<n_1; idx++)
        {
          // between 0 and 1.0
          double val = (double)(rand())/RAND_MAX;
          val = range * (val + shift);
          pArray[idx] = (T)val;
        }
    }
  
    Array1D(const Array1D & rhs)
    {
      n_1 = 0;
      pArray = 0;
      allocate(rhs.n_1);
      for (int i=0; i<n_1; i++) pArray[i] = rhs.pArray[i];
    }

    ~Array1D()
    {
      deallocate();
    }

    T& operator()(int i)
    {
#ifdef QMC_DEBUG
      if( i < 0 || i >= n_1 )
        {
          cerr << "Error: Array1D::operator() has i = " << i
               << " where n_1 = " << n_1 << endl;
          assert( i >= 0 && i < n_1 );
        }
      assert(pArray);
#endif
      return pArray[i];
    }

    T get(int i) const
      {
#ifdef QMC_DEBUG
        assert( i >= 0 && i < n_1 );
        assert(pArray);
#endif
        return pArray[i];
      }

    T& operator[](int i)
    {
#ifdef QMC_DEBUG
      assert( i >= 0 && i < n_1 );
      assert(pArray);
#endif
      return pArray[i];
    }

    friend ostream& operator<<(ostream & strm, const Array1D<T> & rhs)
    {
      rhs.printArray1D(strm, 15, 5, 8, ',', false);
      return strm;
    }

  void printArray1D(ostream & strm, int numSigFig, int columnSpace, int maxI, char sep, bool sci) const
    {
      strm.precision(4);
      if(maxI < 0) maxI = n_1;
      strm << "{";
      for(int i=0; i<n_1;i++)
        {
          if(i % maxI == 0 && i != 0)
            strm << endl << " ";

          if(sci) strm.setf(ios::scientific);
          else strm.unsetf(ios::scientific);

          strm << setprecision(numSigFig)
               << setw(numSigFig+columnSpace)
               << pArray[i];
          if(i<n_1-1) strm << sep;
        }
      strm << " }";
    }

    int read(istream & strm, T initialization, string label)
      {
          int pi = 0;
          strm >> ws;
          string temp;
          while(isdigit(strm.peek()) ||
                strm.peek() == '+' ||
                strm.peek() == '-' ||
                strm.peek() == '.')
            {
              strm >> temp;
              if(pi < n_1)
                pArray[pi] = (T)(atof(temp.c_str()));
              pi++;
              strm >> ws;
            }

          if(pi > n_1)
            clog << "Warning: you entered " << pi << " parameters, when only "
                 << n_1 << " were expected for " << label << endl;  

          for (int i=pi; i < n_1; i++)
            pArray[i] = initialization;

          return pi;
      }
  };

#endif

---