NumArray.h

/*  _________________________________________________________________________
 *
 *  UTILIB: A utility library for developing portable C++ codes.
 *  Copyright (c) 2001, Sandia National Laboratories.
 *  This software is distributed under the GNU Lesser General Public License.
 *  For more information, see the README file in the top UTILIB directory.
 *  _________________________________________________________________________
 */

//
// NumArray.h
//
#ifndef __NumArray_h
#define __NumArray_h

#if !defined (ANSI_HDRS) || defined(SGI)
#include <math.h>
#else
#include <cmath>
#endif

#include "SimpleArray.h"

#ifndef SWIG
template <class T>
class UTILIB_API NumArray;

template <class T>
class UTILIB_API Num2DArray;
#endif // SWIG


#ifdef SWIG
template <class T>
class NumArray: public SimpleArray<T> {
#else
template <class T>
class UTILIB_API NumArray: public SimpleArray<T> {
#endif

  friend Num2DArray<T>;

public:

  NumArray() : SimpleArray<T>() {}
#ifndef NON_ANSI
  explicit 
#endif
  NumArray(const size_type len, T *d=((T*)0), const EnumDataOwned o=DataNotOwned)
                : SimpleArray<T>(len,d,o)
                {if (d == NULL) initialize(Data,0,len);}
  NumArray(const NumArray& array)
                : SimpleArray<T>(array) {}
  ~NumArray() {}

  NumArray<T>& operator=(const NumArray<T>& array)
                {SimpleArray<T>::operator=(array); return *this;}
  NumArray<T>& operator=(const T& array)
                {SimpleArray<T>::operator=(array); return *this;}

  void plus  (const NumArray<T>& x, const NumArray<T>& y);
  void minus (const NumArray<T>& x, const NumArray<T>& y);
  void times (const NumArray<T>& x, const NumArray<T>& y);
  void divide(const NumArray<T>& x, const NumArray<T>& y);

  void plus  (const NumArray<T>& x, const T& z);
  void minus (const NumArray<T>& x, const T& z);
  void times (const NumArray<T>& x, const T& z);
  void divide(const NumArray<T>& x, const T& z);

  NumArray<T> operator-    ();
 
  NumArray<T>& operator+= (const NumArray<T>& x);
  NumArray<T>& operator-= (const NumArray<T>& x);
  NumArray<T>& operator*= (const NumArray<T>& x);
  NumArray<T>& operator/= (const NumArray<T>& x);
 
  NumArray<T>& operator+=  (const T& z);
  NumArray<T>& operator-=  (const T& z);
  NumArray<T>& operator*=  (const T& z);
  NumArray<T>& operator/=  (const T& z);

protected:

  void initialize(T* data, const size_type start, const size_type stop);
};
 


#ifndef SWIG
#define BINARYOP(opname,pseudonym)\
template <class T>\
inline UTILIB_API NumArray<T> opname (const NumArray<T>& a1, const NumArray<T>& a2)\
{\
NumArray<T> res;\
res.pseudonym(a1,a2);\
return res;\
}\
\
template <class T>\
inline UTILIB_API NumArray<T> opname (const NumArray<T>& a1, const T& val)\
{\
NumArray<T> res;\
res.pseudonym(a1,val);\
return res;\
}\
\
template <class T>\
inline UTILIB_API NumArray<T> opname (const T& val, const NumArray<T>& a1)\
{\
NumArray<T> res;\
res.pseudonym(a1,val);\
return res;\
}

BINARYOP( operator+ , plus )
BINARYOP( operator- , minus )
BINARYOP( operator* , times )
BINARYOP( operator/ , divide )

#undef BINARYOP
#endif // SWIG


/*
template <class T>
T NumArray<T>::operator%(const NumArray<T>& a1, const NumArray<T>& a2)
{
T cum;
T *ptr1, *ptr2;
ptr1 = a1.a->Data;
ptr2 = a2.a->Data;
cum = *ptr1 * *ptr2;
ptr1++; ptr2++;
for (size_type i=1; i<a1.a->len; i++, ptr1++, ptr2++)
  cum += *ptr1 * *ptr2;
return cum;
}
*/


#define BINARYOP(opname,op1,pseudonym, op)\
template <class T>\
void NumArray<T>::pseudonym(const NumArray<T>& a1, const NumArray<T>& a2)\
{\
if (a1.size() != a2.size())\
   ErrAbort(errmsg("NumArray<T>::pseudonym : Unequal array lengths %d and %d", a1.size(),a2.size()));\
resize(a1.size());\
T* tmp = Data;\
size_type mylen = size();\
for (size_type i=0; i<mylen; i++)\
  tmp[i] = a1[i] op a2[i];\
}\
\
template <class T>\
void NumArray<T>::pseudonym(const NumArray<T>& a1, const T& val)\
{\
resize(a1.size());\
T* tmp = Data;\
size_type mylen = size();\
for (size_type i=0; i<mylen; i++)\
  tmp[i] = a1[i] op val;\
}\
\
template <class T>\
NumArray<T>& NumArray<T>::opname(const T& val)\
{\
T* tmp = Data;\
size_type mylen = size();\
for (size_type i=0; i<mylen; i++) \
  tmp[i] op1 val;\
return(*this); \
}\
\
template <class T>\
NumArray<T>& NumArray<T>::opname(const NumArray<T>& array)\
{\
T* tmp = Data;\
T* arraytmp = array.Data;\
size_type mylen = size();\
for (size_type i=0; i<mylen; i++) \
  tmp[i] op1 arraytmp[i]; \
return(*this); \
}


BINARYOP(operator+=,+=,plus, + )
BINARYOP(operator-=,-=,minus, - )
BINARYOP(operator*=,*=,times, * )
BINARYOP(operator/=,/=,divide, / )

#undef BINARYOP


template <class T>
NumArray<T> NumArray<T>::operator-    ()
{
size_type mylen = size();
NumArray<T> res(mylen);
T* tmp = Data;
T* restmp = res.Data;
for (size_type i=0; i<mylen; i++)
  restmp[i] = - tmp[i];
return res;
}


template <class T>
void NumArray<T>::initialize(T* data, const size_type start, 
                                                        const size_type stop)
{
for (size_type i=start; i < stop; i++)
  data[i] = ( T ) 0;
}
#endif