matpackI.h

Go to the documentation of this file.

/* Copyright (C) 2001 Stefan Buehler <sbuehler@uni-bremen.de>

   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 <iostream>
#include <iomanip>
#include "arts.h"

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 SparseMatrixView:
class SparseMatrixView;

class Range{
public:
  // Constructors:
  Range(Index start, Index extent, Index stride=1);
  Range(Index start, Joker joker, Index stride=1);
  Range(Joker joker, 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 ConstIterator2D;
  friend class ConstIterator3D;
  friend class SparseMatrixView;
  friend class ConstTensor3View;
  friend class Tensor3View;
  friend class Tensor3;
  friend std::ostream& operator<<(std::ostream& os, const SparseMatrixView& v);


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

class Iterator1D {
public:
  // Constructors:
  Iterator1D();
  Iterator1D(const Iterator1D& o);
  Iterator1D(Numeric *x, Index stride);

  // Operators:
  Iterator1D& operator++();
  Numeric& operator*() const;
  bool operator!=(const Iterator1D& other) const;

private:
  Numeric *mx;
  Index mstride;
};

class ConstIterator1D {
public:
  // Constructors:
  ConstIterator1D();
  ConstIterator1D(const ConstIterator1D& o);
  ConstIterator1D(Numeric *x, Index stride);

  // Operators:
  ConstIterator1D& operator++();
  const Numeric& operator*() const;
  bool operator!=(const ConstIterator1D& other) const;

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

// Declare the vector class:
class Vector;

// Declare the MatrixView class
class MatrixView;

// Declare the ConstMatrixView class
class ConstMatrixView;

class ConstVectorView {
public:

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

  // Const index operators:
  Numeric  operator[](Index n) const;
  ConstVectorView operator[](const Range& r) const;

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

  // Friends:
  friend class VectorView;
  friend class ConstIterator2D;
  friend class ConstMatrixView;
  friend class ConstTensor3View;

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:

  // Const index operators:
  Numeric  operator[](Index n) const;
  ConstVectorView operator[](const Range& r) const;
  // Index Operators:
  Numeric& operator[](Index n);
  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();

  // Friends:
  friend class ConstIterator2D;
  friend class Iterator2D;
  friend class MatrixView;
  friend class Tensor3View;

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

class Iterator2D {
public:
  // Constructors:
  Iterator2D();
  Iterator2D(const Iterator2D& o);
  Iterator2D(const VectorView& x, Index stride);

  // Operators:
  Iterator2D& operator++();
  bool operator!=(const Iterator2D& other) const;
  VectorView* const operator->();
  VectorView& operator*();
  
private:
  VectorView msv;
  Index mstride;
};

class ConstIterator2D {
public:
  // Constructors:
  ConstIterator2D();
  ConstIterator2D(const ConstIterator2D& o);
  ConstIterator2D(const ConstVectorView& x, Index stride);

  // Operators:
  ConstIterator2D& operator++();
  bool operator!=(const ConstIterator2D& other) const;
  const ConstVectorView* operator->() const;
  const ConstVectorView& operator*() const;

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=(VectorView x);
  Vector& operator=(const Vector& v);
  Vector& operator=(const Array<Numeric>& v);
  Vector& operator=(Numeric x);

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

  // Destructor:
  ~Vector();
};

// Declare class Matrix:
class Matrix;


class ConstMatrixView {
public:
  // Member functions:
  Index nrows() const;
  Index ncols() const;

  // Const index operators:
  Numeric  operator()(Index r, Index c) const;
  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;
  
  // Friends:
  friend class ConstVectorView;
  friend class MatrixView;
  friend class ConstIterator3D;
  friend class ConstTensor3View;
  friend ConstMatrixView transpose(ConstMatrixView m);

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:

  // Const index operators:
  Numeric  operator()(Index r, Index c) const;
  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);
  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);

  // Friends:
  friend class VectorView;
  friend class Iterator3D;
  friend class Tensor3View;
  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:
  ~Matrix();
};


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

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

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

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

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



// Declare the existance of class Array:
template<class base>
class Array;

typedef Array<Vector> ArrayOfVector;

typedef Array<Matrix> ArrayOfMatrix;


// Functions for Range:
// --------------------

/* Explicit constructor. 

  \param Start must be >= 0.

  \param Extent also. Although internally negative extent means "to the end",
  this can not be created this way, only with the joker. Zero
  extent is allowed, though, which corresponds to an empty range.

  \param Stride can be anything. It can be omitted, in which case the
  default value is 1. */
inline Range::Range(Index start, Index extent, Index stride) :
  mstart(start), mextent(extent), mstride(stride)
{
  // Start must be >= 0:
  assert( 0<=mstart );
  // Extent also. Although internally negative extent means "to the end",
  // this can not be created this way, only with the joker. Zero
  // extent is allowed, though, which corresponds to an empty range.
  assert( 0<=mextent );
  // Stride can be anything except 0.
  // SAB 2001-09-21: Allow 0 stride.
  //  assert( 0!=mstride);
}

inline Range::Range(Index start, Joker /* joker */, Index stride) :
  mstart(start), mextent(-1), mstride(stride)
{
  // Start must be >= 0:
  assert( 0<=mstart );
}

inline Range::Range(Joker /* joker */, Index stride) :
  mstart(0), mextent(-1), mstride(stride)
{
  // Nothing to do here.
}

inline Range::Range(Index max_size, const Range& r) :
  mstart(r.mstart),
  mextent(r.mextent),
  mstride(r.mstride)
{
  // Start must be >= 0:
  assert( 0<=mstart );
  // ... and < max_size:
  assert( mstart<max_size );

  // Stride must be != 0:
  assert( 0!=mstride);
  
  // Convert negative extent (joker) to explicit extent
  if ( mextent<0 )
    {
      if ( 0<mstride )
        mextent = 1 + (max_size-1-mstart)/mstride;
      else
        mextent = 1 + (0-mstart)/mstride;
    }
  else
    {
#ifndef NDEBUG
      // Check that extent is ok:
      Index fin = mstart+(mextent-1)*mstride;
      assert( 0   <= fin );
      assert( fin <  max_size );
#endif
    }
}

inline Range::Range(const Range& p, const Range& n) :
  mstart(p.mstart + n.mstart*p.mstride),
  mextent(n.mextent),
  mstride(p.mstride*n.mstride)
{
  // We have to juggle here a bit with previous, new, and resulting
  // quantities. I.e.;
  // p.mstride: Previous stride
  // n.mstride: New stride (as specified)
  // mstride:   Resulting stride (old*new)

  // Get the previous final element:
  Index prev_fin = p.mstart+(p.mextent-1)*p.mstride;

  // Resulting start must be >= previous start:
  assert( p.mstart<=mstart );
  // and <= prev_fin:
  assert( mstart<=prev_fin );

  // Resulting stride must be != 0:
  assert( 0!=mstride);

  // Convert negative extent (joker) to explicit extent
  if ( mextent<0 )
    {
      if ( 0<mstride )
        mextent = 1 + (prev_fin-mstart)/mstride;
      else
        mextent = 1 + (p.mstart-mstart)/mstride;
    }
  else
    {
#ifndef NDEBUG
      // Check that extent is ok:
      Index fin = mstart+(mextent-1)*mstride;
      assert( p.mstart <= fin      );
      assert( fin      <= prev_fin );
#endif
    }

}


// Functions for Iterator1D
// ------------------------

inline Iterator1D::Iterator1D()
{
  // Nothing to do here.
}

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

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

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

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

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

// Functions for ConstIterator1D
// -----------------------------

inline ConstIterator1D::ConstIterator1D()
{
  // Nothing to do here.
}

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

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

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

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

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

// Functions for Iterator2D
// ------------------------

inline Iterator2D::Iterator2D()
{
  // Nothing to do here.
}

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

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

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

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

inline VectorView* const Iterator2D::operator->()
{
  return &msv;
}

inline VectorView& Iterator2D::operator*()
{
  return msv;
}

// Functions for ConstIterator2D
// -----------------------------

inline ConstIterator2D::ConstIterator2D()
{
  // Nothing to do here.
}

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

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

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

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

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

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


// Functions for ConstVectorView:
// ------------------------------

inline Index ConstVectorView::nelem() const
{
  return mrange.mextent;
}

inline Numeric ConstVectorView::sum() const
{
  Numeric s=0;
  ConstIterator1D i = begin();
  const ConstIterator1D e = end();

  for ( ; i!=e; ++i )
    s += *i;

  return s;
}

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

inline ConstVectorView ConstVectorView::operator[](const Range& r) const
{
  return ConstVectorView(mdata, mrange, r);
}

inline ConstIterator1D ConstVectorView::begin() const
{
  return ConstIterator1D(mdata+mrange.mstart, mrange.mstride);
}

inline ConstIterator1D ConstVectorView::end() const
{
  return ConstIterator1D( mdata +
                          mrange.mstart +
                          (mrange.mextent)*mrange.mstride,
                          mrange.mstride );
}

inline ConstVectorView::operator ConstMatrixView() const
{
  return ConstMatrixView(mdata,mrange,Range(mrange.mstart,1));
}

inline ConstVectorView::ConstVectorView() :
  mrange(0,0), mdata(NULL)
{
  // Nothing to do here.
}

inline ConstVectorView::ConstVectorView(Numeric *data,
                                        const Range& range) :
  mrange(range),
  mdata(data)
{
  // Nothing to do here.
}

inline ConstVectorView::ConstVectorView(Numeric *data,
                                        const Range& p,
                                        const Range& n) :
  mrange(p,n),
  mdata(data)
{
  // Nothing to do here.
}

inline std::ostream& operator<<(std::ostream& os, const ConstVectorView& v)
{
  ConstIterator1D i=v.begin();
  const ConstIterator1D end=v.end();

  if ( i!=end )
    {
      os << setw(3) << *i;
      ++i;
    }
  for ( ; i!=end ; ++i )
    {
      os << "\n" << setw(3) << *i;
    }

  return os;
}


// Functions for VectorView:
// ------------------------

inline Numeric VectorView::operator[](Index n) const
{
  return ConstVectorView::operator[](n);
}

inline ConstVectorView VectorView::operator[](const Range& r) const
{
  return ConstVectorView::operator[](r);
}

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

inline VectorView VectorView::operator[](const Range& r)
{
  return VectorView(mdata, mrange, r);
}

inline ConstIterator1D VectorView::begin() const
{
  return ConstVectorView::begin();
}

inline ConstIterator1D VectorView::end() const
{
  return ConstVectorView::end();
}

inline Iterator1D VectorView::begin()
{
  return Iterator1D(mdata+mrange.mstart, mrange.mstride);
}

inline Iterator1D VectorView::end()
{
  return Iterator1D( mdata +
                     mrange.mstart +
                     (mrange.mextent)*mrange.mstride,
                     mrange.mstride );
}

inline VectorView VectorView::operator=(const ConstVectorView& v)
{
  //  cout << "Assigning VectorView from ConstVectorView.\n";

  // Check that sizes are compatible:
  assert(mrange.mextent==v.mrange.mextent);
  
  copy( v.begin(), v.end(), begin() );

  return *this;
}

inline VectorView VectorView::operator=(const VectorView& v)
{
  //  cout << "Assigning VectorView from VectorView.\n";

  // Check that sizes are compatible:
  assert(mrange.mextent==v.mrange.mextent);
  
  copy( v.begin(), v.end(), begin() );

  return *this;
}

inline VectorView VectorView::operator=(const Vector& v)
{
  //  cout << "Assigning VectorView from Vector.\n";

  // Check that sizes are compatible:
  assert(mrange.mextent==v.mrange.mextent);
  
  copy( v.begin(), v.end(), begin() );

  return *this;
}

inline VectorView VectorView::operator=(Numeric x)
{
  copy( x, begin(), end() );
  return *this;
}

inline VectorView VectorView::operator*=(Numeric x)
{
  const Iterator1D e=end();
  for ( Iterator1D i=begin(); i!=e ; ++i )
    *i *= x;
  return *this;
}

inline VectorView VectorView::operator/=(Numeric x)
{
  const Iterator1D e=end();
  for ( Iterator1D i=begin(); i!=e ; ++i )
    *i /= x;
  return *this;
}

inline VectorView VectorView::operator+=(Numeric x)
{
  const Iterator1D e=end();
  for ( Iterator1D i=begin(); i!=e ; ++i )
    *i += x;
  return *this;
}

inline VectorView VectorView::operator-=(Numeric x)
{
  const Iterator1D e=end();
  for ( Iterator1D i=begin(); i!=e ; ++i )
    *i -= x;
  return *this;
}

inline VectorView VectorView::operator*=(const ConstVectorView& x)
{
  assert( nelem()==x.nelem() );

  ConstIterator1D s=x.begin();

  Iterator1D i=begin();
  const Iterator1D e=end();

  for ( ; i!=e ; ++i,++s )
    *i *= *s;
  return *this;
}

inline VectorView VectorView::operator/=(const ConstVectorView& x)
{
  assert( nelem()==x.nelem() );

  ConstIterator1D s=x.begin();

  Iterator1D i=begin();
  const Iterator1D e=end();

  for ( ; i!=e ; ++i,++s )
    *i /= *s;
  return *this;
}

inline VectorView VectorView::operator+=(const ConstVectorView& x)
{
  assert( nelem()==x.nelem() );

  ConstIterator1D s=x.begin();

  Iterator1D i=begin();
  const Iterator1D e=end();

  for ( ; i!=e ; ++i,++s )
    *i += *s;
  return *this;
}

inline VectorView VectorView::operator-=(const ConstVectorView& x)
{
  assert( nelem()==x.nelem() );

  ConstIterator1D s=x.begin();

  Iterator1D i=begin();
  const Iterator1D e=end();

  for ( ; i!=e ; ++i,++s )
    *i -= *s;
  return *this;
}

inline VectorView::operator MatrixView()
{
  return MatrixView(mdata,mrange,Range(mrange.mstart,1));
}

inline VectorView::VectorView() :
  ConstVectorView()
{
  // Nothing to do here.
}

inline VectorView::VectorView(Numeric *data,
                            const Range& range) :
  ConstVectorView(data,range)
{
  // Nothing to do here.
}

inline VectorView::VectorView(Numeric *data,
                              const Range& p,
                              const Range& n) :
  ConstVectorView(data,p,n)
{
  // Nothing to do here.
}

inline void copy(ConstIterator1D origin,
                 const ConstIterator1D& end,
                 Iterator1D target)
{
  for ( ; origin!=end ; ++origin,++target )
    *target = *origin;
}

inline void copy(Numeric x,
                 Iterator1D target,
                 const Iterator1D& end)
{
  for ( ; target!=end ; ++target )
    *target = x;
}


// Functions for Vector:
// ---------------------

inline Vector::Vector() 
{
  // Nothing to do here
}

inline Vector::Vector(Index n) :
  VectorView( new Numeric[n],
             Range(0,n))
{
  // Nothing to do here.
}

inline Vector::Vector(Index n, Numeric fill) :
  VectorView( new Numeric[n],
             Range(0,n))
{
  // Here we can access the raw memory directly, for slightly
  // increased efficiency:
  for ( Numeric *x=mdata; x<mdata+n; ++x )
    *x = fill;
}

inline Vector::Vector(Numeric start, Index extent, Numeric stride) :
  VectorView( new Numeric[extent],
             Range(0,extent))
{
  // Fill with values:
  Numeric x = start;
  Iterator1D       i=begin();
  const Iterator1D e=end();
  for ( ; i!=e; ++i )
    {
      *i = x;
      x += stride;
    }
}

inline Vector::Vector(const ConstVectorView& v) :
  VectorView( new Numeric[v.nelem()],
              Range(0,v.nelem()))
{
  copy(v.begin(),v.end(),begin());
}

inline Vector::Vector(const Vector& v) :
  VectorView( new Numeric[v.nelem()],
              Range(0,v.nelem()))
{
  copy(v.begin(),v.end(),begin());
}

inline Vector& Vector::operator=(const Vector& v)
{
  //  cout << "Assigning VectorView from Vector View.\n";

  // Check that sizes are compatible:
  assert(mrange.mextent==v.mrange.mextent);
  copy( v.begin(), v.end(), begin() );
  return *this;
}

inline Vector& Vector::operator=(const Array<Numeric>& x)
{
  VectorView::operator=(x);
  return *this;
}

inline Vector& Vector::operator=(Numeric x)
{
  VectorView::operator=(x);
  return *this;
}

// /** Assignment operator from VectorView. Assignment operators are not
//     inherited. */  
// inline Vector& Vector::operator=(const VectorView v)
// {
//   cout << "Assigning Vector from Vector View.\n";
//  // Check that sizes are compatible:
//   assert(mrange.mextent==v.mrange.mextent);
//   VectorView::operator=(v);
//   return *this;
// }

inline void Vector::resize(Index n)
{
  if ( mrange.mextent != n )
    {
      delete mdata;
      mdata = new Numeric[n];
      mrange.mstart = 0;
      mrange.mextent = n;
      mrange.mstride = 1;
    }
}

inline Vector::~Vector()
{
  delete [] mdata;
}


// Functions for ConstMatrixView:
// ------------------------------

inline Index ConstMatrixView::nrows() const
{
  return mrr.mextent;
}

inline Index ConstMatrixView::ncols() const
{
  return mcr.mextent;
}

inline Numeric ConstMatrixView::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 );
}

inline ConstMatrixView ConstMatrixView::operator()(const Range& r,
                                                   const Range& c) const
{
  return ConstMatrixView(mdata, mrr, mcr, r, c);
}

inline ConstVectorView ConstMatrixView::operator()(const Range& r, Index c) const
{
  // Check that c is valid:
  assert( 0 <= c );
  assert( c < mcr.mextent );

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

inline ConstVectorView ConstMatrixView::operator()(Index r, const Range& c) const
{
  // Check that r is valid:
  assert( 0 <= r );
  assert( r < mrr.mextent );

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

inline ConstIterator2D ConstMatrixView::begin() const
{
  return ConstIterator2D(ConstVectorView(mdata+mrr.mstart,
                                         mcr),
                         mrr.mstride);
}

inline ConstIterator2D ConstMatrixView::end() const
{
  return ConstIterator2D( ConstVectorView(mdata + mrr.mstart + (mrr.mextent)*mrr.mstride,
                                          mcr),
                          mrr.mstride );
}

inline ConstMatrixView::ConstMatrixView() :
  mrr(0,0,1), mcr(0,0,1), mdata(NULL)
{
  // Nothing to do here.
}

inline ConstMatrixView::ConstMatrixView(Numeric *data,
                                        const Range& rr,
                                        const Range& cr) :
  mrr(rr),
  mcr(cr),
  mdata(data)
{
  // Nothing to do here.
}

inline ConstMatrixView::ConstMatrixView(Numeric *data,
                                        const Range& pr, const Range& pc,
                                        const Range& nr, const Range& nc) :
  mrr(pr,nr),
  mcr(pc,nc),
  mdata(data)
{
  // Nothing to do here.
}

inline std::ostream& operator<<(std::ostream& os, const ConstMatrixView& v)
{
  // Row iterators:
  ConstIterator2D ir=v.begin();
  const ConstIterator2D end_row=v.end();

  if ( ir!=end_row )
    {
      ConstIterator1D ic =  ir->begin();
      const ConstIterator1D end_col = ir->end();

      if ( ic!=end_col )
        {
          os << setw(3) << *ic;
          ++ic;
        }
      for ( ; ic!=end_col ; ++ic )
        {
          os << " " << setw(3) << *ic;
        }
      ++ir;
    }
  for ( ; ir!=end_row ; ++ir )
    {
      ConstIterator1D ic =  ir->begin();
      const ConstIterator1D end_col = ir->end();

      os << "\n";
      if ( ic!=end_col )
        {
          os << setw(3) << *ic;
          ++ic;
        }
      for ( ; ic!=end_col ; ++ic )
        {
          os << " " << setw(3) << *ic;
        }
    }

  return os;
}


// Functions for MatrixView:
// -------------------------

inline Numeric MatrixView::operator()(Index r, Index c) const
{
  return ConstMatrixView::operator()(r,c);
}

inline ConstMatrixView MatrixView::operator()(const Range& r, const Range& c) const
{
  return ConstMatrixView::operator()(r,c);  
}

inline ConstVectorView MatrixView::operator()(const Range& r, Index c) const
{
  return ConstMatrixView::operator()(r,c);
}

inline ConstVectorView MatrixView::operator()(Index r, const Range& c) const
{
  return ConstMatrixView::operator()(r,c);
}

inline Numeric& MatrixView::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 );
}

inline MatrixView MatrixView::operator()(const Range& r, const Range& c)
{
  return MatrixView(mdata, mrr, mcr, r, c);
}

inline VectorView MatrixView::operator()(const Range& r, Index c)
{
  // Check that c is valid:
  assert( 0 <= c );
  assert( c < mcr.mextent );

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

inline VectorView MatrixView::operator()(Index r, const Range& c)
{
  // Check that r is valid:
  assert( 0 <= r );
  assert( r <  mrr.mextent );

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

inline ConstIterator2D MatrixView::begin() const
{
  return ConstMatrixView::begin();
}

inline ConstIterator2D MatrixView::end() const
{
  return ConstMatrixView::end();
}

inline Iterator2D MatrixView::begin()
{
  return Iterator2D(VectorView(mdata+mrr.mstart, mcr),
                    mrr.mstride);
}

inline Iterator2D MatrixView::end()
{
  return Iterator2D( VectorView(mdata + mrr.mstart + (mrr.mextent)*mrr.mstride,
                                mcr),
                     mrr.mstride );
}

inline MatrixView& MatrixView::operator=(const ConstMatrixView& m)
{
  // Check that sizes are compatible:
  assert(mrr.mextent==m.mrr.mextent);
  assert(mcr.mextent==m.mcr.mextent);

  copy( m.begin(), m.end(), begin() );
  return *this;
}

inline MatrixView& MatrixView::operator=(const MatrixView& m)
{
  // Check that sizes are compatible:
  assert(mrr.mextent==m.mrr.mextent);
  assert(mcr.mextent==m.mcr.mextent);

  copy( m.begin(), m.end(), begin() );
  return *this;
}

inline MatrixView& MatrixView::operator=(const Matrix& m)
{
  // Check that sizes are compatible:
  assert(mrr.mextent==m.mrr.mextent);
  assert(mcr.mextent==m.mcr.mextent);

  copy( m.begin(), m.end(), begin() );
  return *this;
}

inline MatrixView& MatrixView::operator=(const ConstVectorView& v)
{
  // Check that sizes are compatible:
  assert( mrr.mextent==v.nelem() );
  assert( mcr.mextent==1         );
  //  dummy = ConstMatrixView(v.mdata,v.mrange,Range(v.mrange.mstart,1));;
  ConstMatrixView dummy(v);
  copy( dummy.begin(), dummy.end(), begin() );
  return *this;
}

inline MatrixView& MatrixView::operator=(Numeric x)
{
  copy( x, begin(), end() );
  return *this;
}

inline MatrixView& MatrixView::operator*=(Numeric x)
{
  const Iterator2D er=end();
  for ( Iterator2D r=begin(); r!=er ; ++r )
    {
      const Iterator1D ec = r->end();
      for ( Iterator1D c = r->begin(); c!=ec ; ++c )
        *c *= x;
    }
  return *this;
}

inline MatrixView& MatrixView::operator/=(Numeric x)
{
  const Iterator2D er=end();
  for ( Iterator2D r=begin(); r!=er ; ++r )
    {
      const Iterator1D ec = r->end();
      for ( Iterator1D c = r->begin(); c!=ec ; ++c )
        *c /= x;
    }
  return *this;
}

inline MatrixView& MatrixView::operator+=(Numeric x)
{
  const Iterator2D er=end();
  for ( Iterator2D r=begin(); r!=er ; ++r )
    {
      const Iterator1D ec = r->end();
      for ( Iterator1D c = r->begin(); c!=ec ; ++c )
        *c += x;
    }
  return *this;
}

inline MatrixView& MatrixView::operator-=(Numeric x)
{
  const Iterator2D er=end();
  for ( Iterator2D r=begin(); r!=er ; ++r )
    {
      const Iterator1D ec = r->end();
      for ( Iterator1D c = r->begin(); c!=ec ; ++c )
        *c -= x;
    }
  return *this;
}

inline MatrixView& MatrixView::operator*=(const ConstMatrixView& x)
{
  assert(nrows()==x.nrows());
  assert(ncols()==x.ncols());
  ConstIterator2D  sr = x.begin();
  Iterator2D        r = begin();
  const Iterator2D er = end();
  for ( ; r!=er ; ++r,++sr )
    {
      ConstIterator1D  sc = sr->begin(); 
      Iterator1D        c = r->begin();
      const Iterator1D ec = r->end();
      for ( ; c!=ec ; ++c,++sc )
        *c *= *sc;
    }
  return *this;
}

inline MatrixView& MatrixView::operator/=(const ConstMatrixView& x)
{
  assert(nrows()==x.nrows());
  assert(ncols()==x.ncols());
  ConstIterator2D  sr = x.begin();
  Iterator2D        r = begin();
  const Iterator2D er = end();
  for ( ; r!=er ; ++r,++sr )
    {
      ConstIterator1D  sc = sr->begin(); 
      Iterator1D        c = r->begin();
      const Iterator1D ec = r->end();
      for ( ; c!=ec ; ++c,++sc )
        *c /= *sc;
    }
  return *this;
}

inline MatrixView& MatrixView::operator+=(const ConstMatrixView& x)
{
  assert(nrows()==x.nrows());
  assert(ncols()==x.ncols());
  ConstIterator2D  sr = x.begin();
  Iterator2D        r = begin();
  const Iterator2D er = end();
  for ( ; r!=er ; ++r,++sr )
    {
      ConstIterator1D  sc = sr->begin(); 
      Iterator1D        c = r->begin();
      const Iterator1D ec = r->end();
      for ( ; c!=ec ; ++c,++sc )
        *c += *sc;
    }
  return *this;
}

inline MatrixView& MatrixView::operator-=(const ConstMatrixView& x)
{
  assert(nrows()==x.nrows());
  assert(ncols()==x.ncols());
  ConstIterator2D  sr = x.begin();
  Iterator2D        r = begin();
  const Iterator2D er = end();
  for ( ; r!=er ; ++r,++sr )
    {
      ConstIterator1D  sc = sr->begin(); 
      Iterator1D        c = r->begin();
      const Iterator1D ec = r->end();
      for ( ; c!=ec ; ++c,++sc )
        *c -= *sc;
    }
  return *this;
}

inline MatrixView& MatrixView::operator*=(const ConstVectorView& x)
{
  assert(nrows()==x.nelem());
  assert(ncols()==1);
  ConstIterator1D  sc = x.begin(); 
  Iterator2D        r = begin();
  const Iterator2D er = end();
  for ( ; r!=er ; ++r,++sc )
    {
      Iterator1D        c = r->begin();
      *c *= *sc;
    }
  return *this;
}

inline MatrixView& MatrixView::operator/=(const ConstVectorView& x)
{
  assert(nrows()==x.nelem());
  assert(ncols()==1);
  ConstIterator1D  sc = x.begin(); 
  Iterator2D        r = begin();
  const Iterator2D er = end();
  for ( ; r!=er ; ++r,++sc )
    {
      Iterator1D        c = r->begin();
      *c /= *sc;
    }
  return *this;
}

inline MatrixView& MatrixView::operator+=(const ConstVectorView& x)
{
  assert(nrows()==x.nelem());
  assert(ncols()==1);
  ConstIterator1D  sc = x.begin(); 
  Iterator2D        r = begin();
  const Iterator2D er = end();
  for ( ; r!=er ; ++r,++sc )
    {
      Iterator1D        c = r->begin();
      *c += *sc;
    }
  return *this;
}

inline MatrixView& MatrixView::operator-=(const ConstVectorView& x)
{
  assert(nrows()==x.nelem());
  assert(ncols()==1);
  ConstIterator1D  sc = x.begin(); 
  Iterator2D        r = begin();
  const Iterator2D er = end();
  for ( ; r!=er ; ++r,++sc )
    {
      Iterator1D        c = r->begin();
      *c -= *sc;
    }
  return *this;
}

inline MatrixView::MatrixView() :
  ConstMatrixView()
{
  // Nothing to do here.
}

inline MatrixView::MatrixView(Numeric *data,
                            const Range& rr,
                            const Range& cr) :
  ConstMatrixView(data, rr, cr)
{
  // Nothing to do here.
}

inline MatrixView::MatrixView(Numeric *data,
                            const Range& pr, const Range& pc,
                            const Range& nr, const Range& nc) :
  ConstMatrixView(data,pr,pc,nr,nc)
{
  // Nothing to do here.
}

inline void copy(ConstIterator2D origin,
                 const ConstIterator2D& end,
                 Iterator2D target)
{
  for ( ; origin!=end ; ++origin,++target )
    {
      ConstIterator1D       o = origin->begin();
      const ConstIterator1D e = origin->end();
      Iterator1D            t = target->begin();
      for ( ; o!=e ; ++o,++t )
        *t = *o;
    }
}

inline void copy(Numeric x,
                 Iterator2D target,
                 const Iterator2D& end)
{
  for ( ; target!=end ; ++target )
    {
      Iterator1D       t = target->begin();
      const Iterator1D e = target->end();
      for ( ; t!=e ; ++t )
        *t = x;
    }
}


// Functions for Matrix:
// ---------------------

inline Matrix::Matrix() :
  MatrixView::MatrixView()
{
  // Nothing to do here. However, note that the default constructor
  // for MatrixView has been called in the initializer list. That is
  // crucial, otherwise internal range objects will not be properly
  // initialized. 
}

inline Matrix::Matrix(Index r, Index c) :
  MatrixView( new Numeric[r*c],
             Range(0,r,c),
             Range(0,c))
{
  // Nothing to do here.
}

inline Matrix::Matrix(Index r, Index c, Numeric fill) :
  MatrixView( new Numeric[r*c],
              Range(0,r,c),
              Range(0,c))
{
  // Here we can access the raw memory directly, for slightly
  // increased efficiency:
  for ( Numeric *x=mdata; x<mdata+r*c; ++x )
    *x = fill;
}

inline Matrix::Matrix(const ConstMatrixView& m) :
  MatrixView( new Numeric[m.nrows()*m.ncols()],
             Range( 0, m.nrows(), m.ncols() ),
             Range( 0, m.ncols() ) )
{
  copy(m.begin(),m.end(),begin());
}

inline Matrix::Matrix(const Matrix& m) :
  MatrixView( new Numeric[m.nrows()*m.ncols()],
             Range( 0, m.nrows(), m.ncols() ),
             Range( 0, m.ncols() ) )
{
  // There is a catch here: If m is an empty matrix, then it will have
  // 0 colunns. But this is used to initialize the stride of the row
  // Range! Thus, this method has to be consistent with the behaviour
  // of Range::Range. For now, Range::Range allows also stride 0.
  copy(m.begin(),m.end(),begin());
}

inline Matrix& Matrix::operator=(const Matrix& m)
{
  //  cout << "Matrix copy: m = " << m.nrows() << " " << m.ncols() << "\n";
  //  cout << "             n = " << nrows() << " " << ncols() << "\n";

  // None of the extents can be zero for a valid matrix, so we just
  // have to check one.
  if ( 0 == mrr.mextent )
    {
      // Adjust if previously empty.
      resize( m.mrr.mextent, m.mcr.mextent ); 
    }
  else
    {
      // Check that sizes are compatible:
      assert( mrr.mextent==m.mrr.mextent );
      assert( mcr.mextent==m.mcr.mextent );
    }

  copy( m.begin(), m.end(), begin() );
  return *this;
}

inline Matrix& Matrix::operator=(Numeric x)
{
  copy( x, begin(), end() );
  return *this;
}

inline Matrix& Matrix::operator=(const ConstVectorView& v)
{
  // Check that sizes are compatible:
  assert( mrr.mextent==v.nelem() );
  assert( mcr.mextent==1         );
  //  dummy = ConstMatrixView(v.mdata,v.mrange,Range(v.mrange.mstart,1));;
  ConstMatrixView dummy(v);
  copy( dummy.begin(), dummy.end(), begin() );
  return *this;
}

inline void Matrix::resize(Index r, Index c)
{
  assert( 0<=r );
  assert( 0<=c );

  if ( mrr.mextent!=r || mcr.mextent!=c )
    {
      delete mdata;
      mdata = new Numeric[r*c];

      mrr.mstart = 0;
      mrr.mextent = r;
      mrr.mstride = c;

      mcr.mstart = 0;
      mcr.mextent = c;
      mcr.mstride = 1;
    }
}

inline Matrix::~Matrix()
{
//   cout << "Destroying a Matrix:\n"
//        << *this << "\n........................................\n";
  delete [] mdata;
}


// Some general Matrix Vector functions:

inline Numeric operator*(const ConstVectorView& a, const ConstVectorView& b)
{
  // Check dimensions:
  assert( a.nelem() == b.nelem() );

  const ConstIterator1D ae = a.end();
  ConstIterator1D       ai = a.begin();
  ConstIterator1D       bi = b.begin();

  Numeric res = 0;
  for ( ; ai!=ae ; ++ai, ++bi )
    res += (*ai) * (*bi);

  return res;
}

inline void mult( VectorView y,
                  const ConstMatrixView& M,
                  const ConstVectorView& x )
{
  // Check dimensions:
  assert( y.nelem() == M.nrows() );
  assert( M.ncols() == x.nelem() );

  // Let's get a 1D iterator for y:
  const Iterator1D ye = y.end();
  Iterator1D       yi = y.begin();

  // ... a 2D iterator for M:
  ConstIterator2D  mi = M.begin();

  // ... and 1D iterators pointing to the start and end of x:
  const ConstIterator1D xs = x.begin();
  const ConstIterator1D xe = x.end();

  // This walks through the rows of y and M:
  for ( ; yi!=ye ; ++yi, ++mi )
    {
      // Compute the scalar product between this row of M and x: 
      Numeric dummy = 0;
      
      // Iterators for row of M:
      ConstIterator1D       ri = mi->begin();

      // ... and for x (always initialized to the start of x):
      ConstIterator1D       xi = xs;

      for ( ; xi!=xe ; ++xi, ++ri )
        {
          dummy += (*ri) * (*xi);
        }

      *yi = dummy;
    }
}

inline void mult( MatrixView A,
                  const ConstMatrixView& B,
                  const ConstMatrixView& C )
{
  // Check dimensions:
  assert( A.nrows() == B.nrows() );
  assert( A.ncols() == C.ncols() );
  assert( B.ncols() == C.nrows() );

  // Let's get the transpose of C, so that we can use 2D iterators to
  // access the columns (= rows of the transpose).
  ConstMatrixView CT = transpose(C);

  const Iterator2D ae = A.end();
  Iterator2D       ai = A.begin();
  ConstIterator2D  bi = B.begin();

  // This walks through the rows of A and B:
  for ( ; ai!=ae ; ++ai, ++bi )
    {
      const Iterator1D ace = ai->end();
      Iterator1D       aci = ai->begin();
      ConstIterator2D  cti = CT.begin();

      // This walks through the columns of A with a 1D iterator, and
      // at the same time through the rows of CT, which are the columns of
      // C, with a 2D iterator:
      for ( ; aci!=ace ; ++aci, ++cti )
        {
          // The operator * is used to compute the scalar product
          // between rows of B and rows of C.transpose().
          *aci = (*bi) * (*cti);
        }
    }
}

inline ConstMatrixView transpose(ConstMatrixView m)
{
  return ConstMatrixView(m.mdata, m.mcr, m.mrr);
}

inline MatrixView transpose(MatrixView m)
{
  return MatrixView(m.mdata, m.mcr, m.mrr);
}

inline void transform( VectorView y,
                       double (&my_func)(double),
                       ConstVectorView x )
{
  // Check dimensions:
  assert( y.nelem()==x.nelem() );

  const ConstIterator1D xe = x.end();
  ConstIterator1D       xi = x.begin();
  Iterator1D            yi = y.begin();
  for ( ; xi!=xe; ++xi, ++yi )
    *yi = my_func(*xi);
}

inline void transform( MatrixView y,
                       double (&my_func)(double),
                       ConstMatrixView x )
{
  // Check dimensions:
  assert( y.nrows()==x.nrows() );
  assert( y.ncols()==x.ncols() );

  const ConstIterator2D rxe = x.end();
  ConstIterator2D        rx = x.begin();
  Iterator2D             ry = y.begin();
  for ( ; rx!=rxe; ++rx, ++ry )
    {
      const ConstIterator1D cxe = rx->end();
      ConstIterator1D        cx = rx->begin();
      Iterator1D             cy = ry->begin();
      for ( ; cx!=cxe; ++cx, ++cy )
        *cy = my_func(*cx);
    }
}

inline Numeric max(const ConstVectorView& x)
{
  // Initial value for max:
  Numeric max = x[0];

  const ConstIterator1D xe = x.end();
  ConstIterator1D       xi = x.begin();

  for ( ; xi!=xe ; ++xi )
    {
      if ( *xi > max )
        max = *xi;
    }

  return max;
}

inline Numeric max(const ConstMatrixView& x)
{
  // Initial value for max:
  Numeric max = x(0,0);

  const ConstIterator2D rxe = x.end();
  ConstIterator2D        rx = x.begin();

  for ( ; rx!=rxe ; ++rx )
    {
      const ConstIterator1D cxe = rx->end();
      ConstIterator1D        cx = rx->begin();

      for ( ; cx!=cxe ; ++cx )
        if ( *cx > max )
          max = *cx;
    }
  
  return max;
}

inline Numeric min(const ConstVectorView& x)
{
  // Initial value for min:
  Numeric min = x[0];

  const ConstIterator1D xe = x.end();
  ConstIterator1D       xi = x.begin();

  for ( ; xi!=xe ; ++xi )
    {
      if ( *xi < min )
        min = *xi;
    }

  return min;
}

inline Numeric min(const ConstMatrixView& x)
{
  // Initial value for min:
  Numeric min = x(0,0);

  const ConstIterator2D rxe = x.end();
  ConstIterator2D        rx = x.begin();

  for ( ; rx!=rxe ; ++rx )
    {
      const ConstIterator1D cxe = rx->end();
      ConstIterator1D        cx = rx->begin();

      for ( ; cx!=cxe ; ++cx )
        if ( *cx < min )
          min = *cx;
    }
  
  return min;
}


#endif    // matpackI_h

---