Num2DArray.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.
 *  _________________________________________________________________________
 */

//
// Num2DArray.h
//
#ifndef __Num2DArray_h
#define __Num2DArray_h

#ifdef __GNUC__
#pragma interface
#endif

#ifndef ANSI_HDRS
#include <math.h>
#else
#include <cmath>
#endif

#include "_math.h"
#include "NumArray.h"
#include "Simple2DArray.h"

template <class T>
class UTILIB_API Num2DArray;


template <class T>
Num2DArray<T> operator%    (const Num2DArray<T> & m1, const Num2DArray<T> & m2)
{
static Num2DArray<T> temp_mat(m1.nrows(), m2.ncols());
matmult(temp_mat, m1, m2);
return temp_mat;
}
 

template <class T>
NumArray<T> operator%    (const Num2DArray<T> & m, const NumArray<T> & v)
{
static NumArray<T> temp_vec(m.nrows());
matmult(temp_vec, m, v);
return temp_vec;
}
 

template <class T>
NumArray<T> operator%    (const NumArray<T> & v, const Num2DArray<T> & m)
{
static NumArray<T> temp_vec(m.ncols());
matmult(temp_vec,v,m);
return temp_vec;
}
 
 
template <class T>
void matmult(Num2DArray<T>& res, const Num2DArray<T> & m1, 
                                        const Num2DArray<T> & m2)
{
if (m1.ncols() != m2.nrows())
   ErrAbort(errmsg("matmult: ncols(M1) != nrows(M2) %d %d",m1.ncols(),
                                                m2.nrows()) );
if ((res.nrows() != m1.nrows()) || (res.ncols() != m2.ncols()))
   ErrAbort("matmult: bad dimensions for result.");
 
T temp;
size_type M1Ncols = m1.ncols();
size_type M1Nrows = m1.nrows();
size_type M2Ncols = m2.ncols();
for (size_type i=0; i<M1Nrows; i++){
  for (size_type j=0; j<M2Ncols; j++){
    temp = m1[i][0] * m2[0][j];
    for (size_type k=1; k<M1Ncols; k++){
      temp += m1[i][k] * m2[k][j];
      }
    res[i][j] = temp;
    }
  }
}
 

template <class T>
void matmult(NumArray<T>& res, const Num2DArray<T>& m, const NumArray<T>& v)
{
if (m.ncols() != v.size())
   ErrAbort( "matmult: ncols(M) != size(V)" );
if (res.size() != m.nrows())
   ErrAbort( "matmult: size(RES) != nrows(M)" );
 
T temp;
size_type NRows = m.nrows();
size_type NCols = v.size();
for (size_type i=0; i<NRows; i++) {
  temp = m[i][0] * v[0];
  for (size_type j=1; j<NCols; j++)
    temp += m[i][j] * v[j];
  res[i] = temp;
  }
}
 
 
template <class T>
void matmult(NumArray<T>& res, const NumArray<T> & v, const Num2DArray<T> & m)
{
if (m.nrows() != v.size())
   ErrAbort( "matmult: size(V) != nrows(M)" );
if (res.size() != m.ncols())
   ErrAbort( "matmult: size(RES) != ncols(M)" );
 
T temp;
size_type NRows = v.size();
size_type NCols = m.ncols();
for (size_type i=0; i<NCols; i++) {
  temp = m[0][i] * v[0];
  for (size_type j=1; j<NRows; j++)
    temp += m[j][i] * v[j];
  res[i] = temp;
  }
}



#ifdef SWIG
template <class T>
class Num2DArray: public Simple2DArray<T>
#else
template <class T>
class UTILIB_API Num2DArray: public Simple2DArray<T>
#endif
{
public:

  Num2DArray()  : Simple2DArray<T>() {}
  Num2DArray(const BasicArray<T>& array, const int nrows=1,
                const EnumDataOwned own=DataNotOwned)
                : Simple2DArray<T>(array,nrows,own) {}
  Num2DArray(const int nrows, const int ncols, T *d=((T*)0),
                const EnumDataOwned own=DataNotOwned)
                : Simple2DArray<T>(nrows,ncols,d,own) {}
  Num2DArray(const int nrows, const int ncols, const BasicArray<T>& array,
                const EnumDataOwned own=DataNotOwned)
                : Simple2DArray<T>(nrows,ncols,array,own) {}
  Num2DArray(const Num2DArray& array)
                : Simple2DArray<T>(array) {}

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

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

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

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

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

  Num2DArray<T>& operator+=  (const NumArray<T>& y);
  Num2DArray<T>& operator-=  (const NumArray<T>& y);
  Num2DArray<T>& operator*=  (const NumArray<T>& y);
  Num2DArray<T>& operator/=  (const NumArray<T>& y);

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

  friend void matmult<>(Num2DArray<T>& res, const Num2DArray<T>&, 
                                const Num2DArray<T>&);
  friend void matmult<>(NumArray<T>& res, const Num2DArray<T>&, 
                                const NumArray<T>&);
  friend void matmult<>(NumArray<T>& res, const NumArray<T>&, 
                                const Num2DArray<T>&);
  friend Num2DArray<T> operator%<>(const Num2DArray<T> &, const Num2DArray<T> &);
  friend NumArray<T> operator%<>(const Num2DArray<T> &, const NumArray<T> &);
  friend NumArray<T> operator%<>(const NumArray<T> &, const Num2DArray<T> &);

};
 



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

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

#undef BINARYOP


#define BINARYOP(opname,op1,pseudonym, op)\
template <class T>\
void Num2DArray<T>::pseudonym(const Num2DArray<T>& a1, const Num2DArray<T>& a2)\
{\
if ((a1.nrows() != a2.nrows()) || (a1.ncols() != a2.ncols()))\
   ErrAbort(errmsg("Num2DArray<T>::pseudonym : Unequal 2D array sizes %dx%d and %dx%d",a1.nrows(),a1.ncols(),a2.nrows(),a2.ncols()));\
resize(a1.nrows(),a2.ncols());\
for (size_type i=0; i<a->Nrows; i++)\
  for (size_type j=0; j<a->Ncols; j++)\
    a->Data[i][j] = a1[i][j] op a2[i][j];\
}\
\
template <class T>\
void Num2DArray<T>::pseudonym(const Num2DArray<T>& a1, const NumArray<T>& array)\
{\
if (a1.ncols() != array.size())\
   ErrAbort(errmsg("Num2DArray<T>::pseudonym : Number of columns %d does not equal array length %d", a1.ncols(), array.size()));\
resize(a1.nrows(),a1.ncols());\
for (size_type i=0; i<a->Nrows; i++)\
  for (size_type j=0; j<a->Ncols; j++)\
    a->Data[i][j] = a1[i][j] op array[j];\
}\
\
template <class T>\
void Num2DArray<T>::pseudonym(const Num2DArray<T>& a1, const T& val)\
{\
resize(a1.nrows(),a1.ncols());\
for (size_type i=0; i<a->Nrows; i++)\
  for (size_type j=0; j<a->Ncols; j++)\
    a->Data[i][j] = a1[i][j] op val;\
}\
template <class T>\
Num2DArray<T>& Num2DArray<T>::opname(const Num2DArray<T>& a1)\
{\
if ((nrows() != a1.nrows()) || (ncols() != a1.ncols()))\
   ErrAbort(errmsg("Num2DArray<T>::opname : Unequal 2D array sizes %dx%d and %dx%d",nrows(),ncols(),a1.nrows(),a1.ncols()));\
for (size_type i=0; i<nrows(); i++)\
  for (size_type j=0; j<ncols(); j++)\
    a->Data[i][j] op1 a1.a->Data[i][j];\
return *this;\
}\
\
template <class T>\
Num2DArray<T>& Num2DArray<T>::opname(const NumArray<T>& a1)\
{\
if (ncols() != a1.size())\
   ErrAbort(errmsg("Num2DArray<T>::opname : Number of columns %d does not equal array length %d", ncols(), a1.size()));\
for (size_type i=0; i<nrows(); i++)\
  for (size_type j=0; j<ncols(); j++)\
    a->Data[i][j] op1 a1.data()[j];\
return *this;\
}\
\
template <class T>\
Num2DArray<T>& Num2DArray<T>::opname(const T& val)\
{\
for (size_type i=0; i<nrows(); i++)\
  for (size_type j=0; j<ncols(); j++)\
    a->Data[i][j] op1 val;\
return *this;\
}

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

#undef BINARYOP


template <class T>
Num2DArray<T> Num2DArray<T>::operator-    ()
{
Num2DArray<T> res;
for (size_type i=0; i<a->Nrows; i++)
  for (size_type j=0; j<a->Ncols; j++)
    res.a->Data[i][j] = - a->Data[i][j];
return res;
}


#endif