matpackI.h

Go to the documentation of this file.

/* Copyright (C) 2001-2008 Stefan Buehler <sbuehler@ltu.se>

   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, 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 matpackI_h
#define matpackI_h

#include "matpack.h"
#include <cassert>
#include "array.h"

// Declare existance of some classes
class bofstream;

class Joker {
  // Nothing here.
};

// Declare existence of the global joker object:
extern Joker joker;

// Declare the existence of class VectorView:
class VectorView;

// Declare the existence of class ConstVectorView:
class ConstVectorView;

// Declare the existence of class ConstMatrixView:
class ConstMatrixView;

class Range{
public:
  // Constructors:
  Range(Index start, Index extent, Index stride=1);
  Range(Index start, Joker      j, Index stride=1);
  Range(Joker     j, Index stride=1);
  Range(Index max_size, const Range& r);
  Range(const Range& p, const Range& n);

  // Friends:
  friend class ConstVectorView;
  friend class VectorView;
  friend class Vector;
  friend class ConstMatrixView;
  friend class MatrixView;
  friend class Matrix;
  friend class Iterator2D;
  friend class Iterator3D;
  friend class Iterator4D;
  friend class Iterator5D;
  friend class Iterator6D;
  friend class Iterator7D;
  friend class ConstIterator2D;
  friend class ConstIterator3D;
  friend class ConstIterator4D;
  friend class ConstIterator5D;
  friend class ConstIterator6D;
  friend class ConstIterator7D;
  friend class ConstTensor3View;
  friend class Tensor3View;
  friend class Tensor3;
  friend class ConstTensor4View;
  friend class Tensor4View;
  friend class Tensor4;
  friend class ConstTensor5View;
  friend class Tensor5View;
  friend class Tensor5;
  friend class ConstTensor6View;
  friend class Tensor6View;
  friend class Tensor6;
  friend class ConstTensor7View;
  friend class Tensor7View;
  friend class Tensor7;
  friend class Sparse;
  friend void mult (VectorView, const ConstMatrixView&,
                    const ConstVectorView&);

  Index get_start () const { return mstart; }
  Index get_extent () const { return mextent; }
  Index get_stride () const { return mstride; }

private:
  Index mstart;
  Index mextent;
  Index mstride;
};

class Iterator1D {
public:
  Iterator1D() : mx(NULL), mstride(0) { /* Nothing to do here. */ }

  Iterator1D(const Iterator1D& o) : mx(o.mx), mstride(o.mstride)
    { /* Nothing to do here. */ }

  Iterator1D(Numeric *x, Index stride) : mx(x), mstride(stride)
      { /* Nothing to do here. */ }

  // Operators:

  Iterator1D& operator++()
    { mx += mstride; return *this; }

  Numeric& operator*() const { return *mx; }

  bool operator!=(const Iterator1D& other) const
    { if (mx != other.mx) return true; else return false; }


private:
  Numeric *mx;
  Index mstride;
};

class ConstIterator1D {
public:
  ConstIterator1D() : mx(NULL), mstride(0)
    { /* Nothing to do here. */ }

  ConstIterator1D(const ConstIterator1D& o) : mx(o.mx), mstride(o.mstride)
      { /* Nothing to do here. */ }

  ConstIterator1D(Numeric *x, Index stride) : mx(x), mstride(stride)
      { /* Nothing to do here. */ }

  // Operators:
  ConstIterator1D& operator++()
    { mx += mstride; return *this; }

  const Numeric& operator*() const { return *mx; }

  bool operator!=(const ConstIterator1D& other) const
    { if (mx != other.mx) return true; else return false; }

private:
  const Numeric *mx;
  Index mstride;
};

// Declare the vector class:
class Vector;

// Declare the MatrixView class
class MatrixView;

class ConstVectorView {
public:
  // Typedef for compatibility with STL
  typedef ConstIterator1D const_iterator;

  // Member functions:
  Index nelem() const;
  Numeric sum() const;

  // Const index operators:
  Numeric operator[](Index n) const
    { // Check if index is valid:
      assert( 0<=n );
      assert( n<mrange.mextent );
      return *( mdata +
                mrange.mstart +
                n*mrange.mstride );
    }

  ConstVectorView operator[](const Range& r) const;

  // Functions returning iterators:
  ConstIterator1D begin() const;
  ConstIterator1D end() const;
  
  // Conversion to 1 column matrix:
  operator ConstMatrixView() const;

  virtual ~ConstVectorView() {};

  // Friends:
  friend class VectorView;
  friend class ConstIterator2D;
  friend class ConstMatrixView;
  friend class ConstTensor3View;
  friend class ConstTensor4View;
  friend class ConstTensor5View;
  friend class ConstTensor6View;
  friend class ConstTensor7View;
  friend int poly_root_solve (Matrix& roots, Vector& coeffs);
  friend void mult (VectorView, const ConstMatrixView&,
                    const ConstVectorView&);
  
  // A special constructor, that allows to make a ConstVectorView of a scalar.
  ConstVectorView(const Numeric& a);

protected:
  // Constructors:
  ConstVectorView();
  ConstVectorView(Numeric *data, const Range& range);
  ConstVectorView(Numeric *data, const Range& p, const Range& n);

  // Data members:
  // -------------
  Range mrange;
  Numeric *mdata;
};

class VectorView : public ConstVectorView {
public:
  VectorView (const Vector&);
  VectorView (Vector& v);

  // Typedef for compatibility with STL
  typedef Iterator1D iterator;

  // Const index operators:
  Numeric operator[](Index n) const
    { return ConstVectorView::operator[](n); }

  ConstVectorView operator[](const Range& r) const;

  Numeric& operator[](Index n)
    { // Check if index is valid:
      assert( 0<=n );
      assert( n<mrange.mextent );
      return *( mdata +
                mrange.mstart +
                n*mrange.mstride );
    }

  VectorView operator[](const Range& r);

  // Constant iterators:
  ConstIterator1D begin() const;
  ConstIterator1D end() const;
  // Iterators:
  Iterator1D begin();
  Iterator1D end();

  // Assignment operators:
  VectorView& operator=(const ConstVectorView& v);
  VectorView& operator=(const VectorView& v);
  VectorView& operator=(const Vector& v);
  VectorView& operator=(const Array<Numeric>& v); 
  VectorView& operator=(Numeric x);

  // Other operators:
  VectorView operator*=(Numeric x);
  VectorView operator/=(Numeric x);
  VectorView operator+=(Numeric x);
  VectorView operator-=(Numeric x);

  VectorView operator*=(const ConstVectorView& x);
  VectorView operator/=(const ConstVectorView& x);
  VectorView operator+=(const ConstVectorView& x);
  VectorView operator-=(const ConstVectorView& x);

  // Conversion to 1 column matrix:
  operator MatrixView();
  // Conversion to a plain C-array
  const Numeric *get_c_array() const;
  Numeric *get_c_array();

  virtual ~VectorView() {};

  // Friends:
  friend class ConstIterator2D;
  friend class Iterator2D;
  friend class MatrixView;
  friend class Tensor3View;
  friend class Tensor4View;
  friend class Tensor5View;
  friend class Tensor6View;
  friend class Tensor7View;

  // A special constructor, that allows to make a VectorView of a scalar.
  VectorView(Numeric& a);


protected:
  // Constructors:
  VectorView();
  VectorView(Numeric *data, const Range& range);
  VectorView(Numeric *data, const Range& p, const Range& n);
};

class Iterator2D {
public:
  // Constructors:
  Iterator2D() : msv(), mstride(0)  { /* Nothing to do here. */ }

  Iterator2D(const Iterator2D& o) : msv(o.msv), mstride(o.mstride)
      { /* Nothing to do here. */ }

  Iterator2D(const VectorView& x, Index stride) : msv(x), mstride(stride)
      { /* Nothing to do here. */ }

  // Operators:
  Iterator2D& operator++() { msv.mdata += mstride; return *this; }

  bool operator!=(const Iterator2D& other) const
    { if ( msv.mdata + msv.mrange.mstart !=
           other.msv.mdata + other.msv.mrange.mstart )
        return true;
      else
        return false;
    }

  VectorView * operator->() { return &msv; }

  VectorView& operator*() { return msv; }

private:
  VectorView msv;
  Index mstride;
};

class ConstIterator2D {
public:
  // Constructors:
  ConstIterator2D() : msv(), mstride(0) { /* Nothing to do here. */ }

  ConstIterator2D(const ConstIterator2D& o) : msv(o.msv), mstride(o.mstride)
      { /* Nothing to do here. */ }

  ConstIterator2D(const ConstVectorView& x, Index stride)
    : msv(x), mstride(stride)
      { /* Nothing to do here. */ }

  // Operators:
  ConstIterator2D& operator++() { msv.mdata += mstride; return *this; }

  bool operator!=(const ConstIterator2D& other) const
    { if ( msv.mdata + msv.mrange.mstart !=
           other.msv.mdata + other.msv.mrange.mstart )
        return true;
      else
        return false;
    }

  const ConstVectorView* operator->() const { return &msv; }

  const ConstVectorView& operator*() const { return msv; }

private:
  ConstVectorView msv;
  Index mstride;
};

class Vector : public VectorView {
public:
  // Constructors:
  Vector();
  explicit Vector(Index n);
  Vector(Index n, Numeric fill);
  Vector(Numeric start, Index extent, Numeric stride);
  Vector(const ConstVectorView& v);
  Vector(const Vector& v);

  // Assignment operators:
  Vector& operator=(const Vector& v);
  Vector& operator=(const Array<Numeric>& v); 
  Vector& operator=(Numeric x);

  // Resize function:
  void resize(Index n);

  // Destructor:
  virtual ~Vector();
};

// Declare class Matrix:
class Matrix;


class ConstMatrixView {
public:
  // Typedef for compatibility with STL
  typedef ConstIterator2D const_iterator;

  // Member functions:
  Index nrows() const;
  Index ncols() const;

  // Const index operators:
  Numeric operator()(Index r, Index c) const
    { // Check if indices are valid:
      assert( 0<=r );
      assert( 0<=c );
      assert( r<mrr.mextent );
      assert( c<mcr.mextent );

      return *( mdata +
                mrr.mstart +
                r*mrr.mstride +
                mcr.mstart +
                c*mcr.mstride );
    }

  ConstMatrixView operator()(const Range& r, const Range& c) const;
  ConstVectorView operator()(const Range& r, Index c) const;
  ConstVectorView operator()(Index r, const Range& c) const;

  // Functions returning iterators:
  ConstIterator2D begin() const;
  ConstIterator2D end() const;

  virtual ~ConstMatrixView() {};

  // Friends:
  friend class MatrixView;
  friend class ConstIterator3D;
  friend class ConstVectorView;
  friend class ConstTensor3View;
  friend class ConstTensor4View;
  friend class ConstTensor5View;
  friend class ConstTensor6View;
  friend class ConstTensor7View;
  friend ConstMatrixView transpose(ConstMatrixView m);
  friend int poly_root_solve (Matrix& roots, Vector& coeffs);
  friend void mult (VectorView, const ConstMatrixView&,
                    const ConstVectorView&);

protected:
  // Constructors:
  ConstMatrixView();
  ConstMatrixView(Numeric *data, const Range& r, const Range& c);
  ConstMatrixView(Numeric *data,
                  const Range& pr, const Range& pc,
                  const Range& nr, const Range& nc);

  // Data members:
  // -------------
  Range mrr;
  Range mcr;
  Numeric *mdata;
};

class MatrixView : public ConstMatrixView {
public:
  // Typedef for compatibility with STL
  typedef Iterator2D iterator;

  // Const index operators:
  Numeric operator()(Index r, Index c) const
    { return ConstMatrixView::operator()(r,c); }

  ConstMatrixView operator()(const Range& r, const Range& c) const;
  ConstVectorView operator()(const Range& r, Index c) const;
  ConstVectorView operator()(Index r, const Range& c) const;
  // Index Operators:
  Numeric& operator()(Index r, Index c)
    { // Check if indices are valid:
      assert( 0<=r );
      assert( 0<=c );
      assert( r<mrr.mextent );
      assert( c<mcr.mextent );

      return *( mdata +
                mrr.mstart +
                r*mrr.mstride +
                mcr.mstart +
                c*mcr.mstride );
    }

  MatrixView operator()(const Range& r, const Range& c);
  VectorView operator()(const Range& r, Index c);
  VectorView operator()(Index r, const Range& c);

  // Functions returning const iterators:
  ConstIterator2D begin() const;
  ConstIterator2D end() const;
  // Functions returning iterators:
  Iterator2D begin();
  Iterator2D end();
  
  // Assignment operators:
  MatrixView& operator=(const ConstMatrixView& v);
  MatrixView& operator=(const MatrixView& v);
  MatrixView& operator=(const Matrix& v);
  MatrixView& operator=(const ConstVectorView& v);
  MatrixView& operator=(Numeric x);

  // Other operators:
  MatrixView& operator*=(Numeric x);
  MatrixView& operator/=(Numeric x);
  MatrixView& operator+=(Numeric x);
  MatrixView& operator-=(Numeric x);

  MatrixView& operator*=(const ConstMatrixView& x);
  MatrixView& operator/=(const ConstMatrixView& x);
  MatrixView& operator+=(const ConstMatrixView& x);
  MatrixView& operator-=(const ConstMatrixView& x);

  MatrixView& operator*=(const ConstVectorView& x);
  MatrixView& operator/=(const ConstVectorView& x);
  MatrixView& operator+=(const ConstVectorView& x);
  MatrixView& operator-=(const ConstVectorView& x);

  // Conversion to a plain C-array
  const Numeric *get_c_array() const;
  Numeric *get_c_array();

  virtual ~MatrixView() {};

  // Friends:
  friend class VectorView;
  friend class Iterator3D;
  friend class Tensor3View;
  friend class Tensor4View;
  friend class Tensor5View;
  friend class Tensor6View;
  friend class Tensor7View;
  friend ConstMatrixView transpose(ConstMatrixView m);
  friend MatrixView transpose(MatrixView m);

protected:
  // Constructors:
  MatrixView();
  MatrixView(Numeric *data, const Range& r, const Range& c);
  MatrixView(Numeric *data,
             const Range& pr, const Range& pc,
             const Range& nr, const Range& nc);
};

class Matrix : public MatrixView {
public:
  // Constructors:
  Matrix();
  Matrix(Index r, Index c);
  Matrix(Index r, Index c, Numeric fill);
  Matrix(const ConstMatrixView& v);
  Matrix(const Matrix& v);

  // Assignment operators:
  Matrix& operator=(const Matrix& x);
  Matrix& operator=(Numeric x);
  Matrix& operator=(const ConstVectorView& v);

  // Resize function:
  void resize(Index r, Index c);

  // Destructor:
  virtual ~Matrix();

  Numeric *get_raw_data() { return mdata; };
};

// Function declarations:
// ----------------------

void copy(ConstIterator1D origin,
          const ConstIterator1D& end,
          Iterator1D target);

void copy(Numeric x,
          Iterator1D target,
          const Iterator1D& end);

void copy(ConstIterator2D origin,
          const ConstIterator2D& end,
          Iterator2D target);

void copy(Numeric x,
          Iterator2D target,
          const Iterator2D& end);

void mult( VectorView y,
           const ConstMatrixView& M,
           const ConstVectorView& x );

void mult( MatrixView A,
           const ConstMatrixView& B,
           const ConstMatrixView& C );

ConstMatrixView transpose(ConstMatrixView m);

MatrixView transpose(MatrixView m);

void transform( VectorView y,
                double (&my_func)(double),
                ConstVectorView x );

void transform( MatrixView y,
                double (&my_func)(double),
                ConstMatrixView x );

Numeric max(const ConstVectorView& x);

Numeric max(const ConstMatrixView& x);

Numeric min(const ConstVectorView& x);

Numeric min(const ConstMatrixView& x);

Numeric mean(const ConstVectorView& x);

Numeric mean(const ConstMatrixView& x);

Numeric operator*(const ConstVectorView& a, const ConstVectorView& b);

ostream& operator<<(ostream& os, const ConstVectorView& v);

ostream& operator<<(ostream& os, const ConstMatrixView& v);

// Helper function for debugging
#ifndef NDEBUG

Numeric debug_matrixview_get_elem (MatrixView& mv, Index r, Index c);

#endif

#endif    // matpackI_h

---