matrix_template.h

Go to the documentation of this file.

/*
 *  Copyright (C) 2002 Stefan Weinzierl
 *
 *  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 __PDF_MATRIX_TEMPLATE_H__
#define __PDF_MATRIX_TEMPLATE_H__

#include <cstddef>
#include <valarray>
#include <iostream>

namespace pdf {

  // forward declaration
  template<class T> class matrix_template;

  template<class T> class vector_template {

        // ctors
    public :  
      vector_template(void);
      explicit vector_template(size_t n);
      vector_template(const T & val, size_t n);
      vector_template(const T* p, size_t n);
      explicit vector_template(const std::valarray<T> & vv);

        // functions 
    public :
      size_t size(void) const;


        // subscripting
    public :  
      T operator[](size_t i) const;
      T & operator[](size_t i);

      T operator()(size_t i) const;
      T & operator()(size_t i);

        // friends
      template<class TT> friend vector_template<TT> operator+ (const vector_template<TT> & v1, const vector_template<TT> & v2);
      template<class TT> friend vector_template<TT> operator- (const vector_template<TT> & v1, const vector_template<TT> & v2);
      template<class TT> friend vector_template<TT> operator+ (const vector_template<TT> & v1);
      template<class TT> friend vector_template<TT> operator- (const vector_template<TT> & v1);
      template<class TT, class SS> friend vector_template<TT> operator* (const SS & c, const vector_template<TT> & v2);
      template<class TT, class SS> friend vector_template<TT> operator* (const vector_template<TT> & v1, const SS & c);
      template<class TT, class SS> friend vector_template<TT> operator/ (const vector_template<TT> & v1, const SS & c);

      template<class TT> friend TT operator* (const vector_template<TT> & v1, const vector_template<TT> & v2);
      template<class TT> friend vector_template<TT> operator* (const matrix_template<TT> & v1, const vector_template<TT> & v2);
      template<class TT> friend vector_template<TT> operator* (const vector_template<TT> & v1, const matrix_template<TT> & v2);
 
        // I/O operators
      template<class TT> friend std::ostream & operator<< (std::ostream & os, const vector_template<TT> & v);

        // member variables :
    protected : 
      std::valarray<T> v;

  };

  template<class T> class matrix_template {

        // ctors
    public :  
      matrix_template(void);
      matrix_template(size_t n, size_t m);
      matrix_template(const T & val, size_t n, size_t m);
      matrix_template(const T* p, size_t n, size_t m);
      matrix_template(const std::valarray<T> & vv, size_t n, size_t m);

        // functions 
    public :
      size_t rows(void) const;
      size_t columns(void) const;


        // subscripting
    public :  
      T operator()(size_t i, size_t j) const;
      T & operator()(size_t i, size_t j);

        // friends
      template<class TT> friend matrix_template<TT> operator+ (const matrix_template<TT> & v1, const matrix_template<TT> & v2);
      template<class TT> friend matrix_template<TT> operator- (const matrix_template<TT> & v1, const matrix_template<TT> & v2);
      template<class TT> friend matrix_template<TT> operator+ (const matrix_template<TT> & v1);
      template<class TT> friend matrix_template<TT> operator- (const matrix_template<TT> & v1);
      template<class TT, class SS> friend matrix_template<TT> operator* (const SS & c, const matrix_template<TT> & v2);
      template<class TT, class SS> friend matrix_template<TT> operator* (const matrix_template<TT> & v1, const SS & c);
      template<class TT, class SS> friend matrix_template<TT> operator/ (const matrix_template<TT> & v1, const SS & c);

      template<class TT> friend matrix_template<TT> operator* (const matrix_template<TT> & v1, const matrix_template<TT> & v2);
      template<class TT> friend vector_template<TT> operator* (const matrix_template<TT> & v1, const vector_template<TT> & v2);
      template<class TT> friend vector_template<TT> operator* (const vector_template<TT> & v1, const matrix_template<TT> & v2);
 
        // I/O operators
      template<class TT> friend std::ostream & operator<< (std::ostream & os, const matrix_template<TT> & v);

        // member variables :
    protected : 
      size_t nrows; 
      size_t ncols; 
      std::valarray<T> v;

  };

  // ----------------------------------------------------------------------------

  // performance-critical inlined method implementations

  // ctors

  template<class T> inline vector_template<T>::vector_template(void) : v()
    {}

  template<class T> inline vector_template<T>::vector_template(size_t n) : v(n)
    {}

  template<class T> inline vector_template<T>::vector_template(const T & val, size_t n) : v(val, n)
    {}

  template<class T> inline vector_template<T>::vector_template(const T* p, size_t n) : v(p,n)
    {}

  template<class T> inline vector_template<T>::vector_template(const std::valarray<T> & vv) : v(vv)
    {}

  // functions 

  template<class T> inline size_t vector_template<T>::size(void) const
    {
      return v.size();
    }

  // subscripting

  template<class T> inline T vector_template<T>::operator[](size_t i) const
    {
      return v[i];
    }

  template<class T> inline T & vector_template<T>::operator[](size_t i)
    {
      return v[i];
    }

  template<class T> inline T vector_template<T>::operator()(size_t i) const
    {
      return v[i];
    }

  template<class T> inline T & vector_template<T>::operator()(size_t i)
    {
      return v[i];
    }

  // friends

  // operators

  template<class T> inline vector_template<T> operator+ (const vector_template<T> & v1, const vector_template<T> & v2)
    {
      return vector_template<T>(v1.v+v2.v);
    }

  template<class T> inline vector_template<T> operator- (const vector_template<T> & v1, const vector_template<T> & v2)
    {
      return vector_template<T>(v1.v-v2.v);
    }

  template<class T> inline vector_template<T> operator+ (const vector_template<T> & v1)
    {
      return vector_template<T>(+v1.v);
    }

  template<class T> inline vector_template<T> operator- (const vector_template<T> & v1)
    {
      return vector_template<T>(-v1.v);
    }

  template<class T, class S> inline vector_template<T> operator* (const S & c, const vector_template<T> & v2)
    {
      T cc = c;
      return vector_template<T>( cc*v2.v);
    }

  template<class T, class S> inline vector_template<T> operator* (const vector_template<T> & v1, const S & c)
    {
      T cc = c;
      return vector_template<T>(v1.v*cc);
    }

  template<class T, class S> inline vector_template<T> operator/ (const vector_template<T> & v1, const S & c)
    {
      T cc = c;
      return vector_template<T>(v1.v/cc);
    }

  template<class T> inline T operator* (const vector_template<T> & v1, const vector_template<T> & v2)
    {
      T res = 0;
      for ( size_t i=0; i<v1.size(); i++) res += v1(i)*v2(i);

      return res;
    }

  // I/O operators

  template<class T> inline std::ostream & operator<< (std::ostream & os, const vector_template<T> & v)
    {
      os << "[";
      for ( size_t i=0; i<v.size(); i++)
        {
          if (i>0) os << ",";
          os << v(i);
        }

      return os << "]";
    }


  // --------------------------------------------------------------

  // ctors

  template<class T> inline matrix_template<T>::matrix_template(void) : nrows(), ncols(), v()
    {}

  template<class T> inline matrix_template<T>::matrix_template(size_t n, size_t m) : nrows(n), ncols(m), v(n*m)
    {}

  template<class T> inline matrix_template<T>::matrix_template(const T & val, size_t n, size_t m) : 
    nrows(n), ncols(m), v(val, n*m)
    {}

  template<class T> inline matrix_template<T>::matrix_template(const T* p, size_t n, size_t m) : 
    nrows(n), ncols(m), v(p,n*m)
    {}

  template<class T> inline matrix_template<T>::matrix_template(const std::valarray<T> & vv, size_t n, size_t m) : 
    nrows(n), ncols(m), v(vv)
    {}

  // functions 

  template<class T> inline size_t matrix_template<T>::rows(void) const
    {
      return nrows;
    }

  template<class T> inline size_t matrix_template<T>::columns(void) const
    {
      return ncols;
    }

  // subscripting

  template<class T> inline T matrix_template<T>::operator()(size_t i, size_t j) const
    {
      // this fixes the physical layout
      return v[i*ncols+j];
    }

  template<class T> inline T & matrix_template<T>::operator()(size_t i, size_t j)
    {
      // this fixes the physical layout
      return v[i*ncols+j];
    }

  // friends

  // operators

  template<class T> inline matrix_template<T> operator+ (const matrix_template<T> & v1, const matrix_template<T> & v2)
    {
      return matrix_template<T>(v1.v+v2.v, v1.nrows, v1.ncols);
    }

  template<class T> inline matrix_template<T> operator- (const matrix_template<T> & v1, const matrix_template<T> & v2)
    {
      return matrix_template<T>(v1.v-v2.v, v1.nrows, v1.ncols);
    }

  template<class T> inline matrix_template<T> operator+ (const matrix_template<T> & v1)
    {
      return matrix_template<T>(+v1.v, v1.nrows, v1.ncols);
    }

  template<class T> inline matrix_template<T> operator- (const matrix_template<T> & v1)
    {
      return matrix_template<T>(-v1.v, v1.nrows, v1.ncols);
    }

  template<class T, class S> inline matrix_template<T> operator* (const S & c, const matrix_template<T> & v2)
    {
      T cc = c;
      return matrix_template<T>(cc*v2.v, v2.nrows, v2.ncols);
    }

  template<class T, class S> inline matrix_template<T> operator* (const matrix_template<T> & v1, const S & c)
    {
      T cc = c;
      return matrix_template<T>(v1.v*cc, v1.nrows, v1.ncols);
    }

  template<class T, class S> inline matrix_template<T> operator/ (const matrix_template<T> & v1, const S & c)
    {
      T cc = c;
      return matrix_template<T>(v1.v/cc, v1.nrows, v1.ncols);
    }

  template<class T> inline matrix_template<T> operator* (const matrix_template<T> & v1, const matrix_template<T> & v2)
    {
      std::valarray<T> res(v1.nrows*v2.ncols);

      for ( size_t i=0; i<v1.nrows; i++) 
        {
          for ( size_t j=0; j<v2.ncols; j++)
            {
              T res_scal = 0;
              for ( size_t k=0; k<v1.ncols; k++) res_scal += v1.v[i*v1.ncols+k]*v2.v[k*v2.ncols+j];
              res[i*v2.ncols+j] = res_scal;
            }
        }

      return matrix_template<T>(res, v1.nrows, v2.ncols);
    }

  template<class T> inline vector_template<T> operator* (const matrix_template<T> & v1, const vector_template<T> & v2)
    {
      std::valarray<T> res(v1.nrows);

      for ( size_t i=0; i<v1.nrows; i++) 
        {
          T res_scal = 0;
          for ( size_t k=0; k<v1.ncols; k++) res_scal += v1.v[i*v1.ncols+k]*v2.v[k];
          res[i] = res_scal;
        }

      return vector_template<T>(res);
    }

  template<class T> inline vector_template<T> operator* (const vector_template<T> & v1, const matrix_template<T> & v2)
    {
      std::valarray<T> res(v2.ncols);

      for ( size_t j=0; j<v2.ncols; j++)
        {
          T res_scal = 0;
          for ( size_t k=0; k<v2.nrows; k++) res_scal += v1.v[k]*v2.v[k*v2.ncols+j];
          res[j] = res_scal;
        }

      return vector_template<T>(res);
    }

  // I/O operators

  template<class T> inline std::ostream & operator<< (std::ostream & os, const matrix_template<T> & v)
    {
      os << "[";
      for ( size_t i=0; i<v.rows(); i++)
        {
          if (i>0) os << ",";
          os << "[";
          for ( size_t j=0; j<v.columns(); j++)
            {
              if (j>0) os << ",";
              os << v(i,j);
            }
          os << "]";
        }

      return os << "]";
    }



} // namespace pdf

#endif // ndef __PDF_MATRIX_TEMPLATE_H__