matpackI.cc

Go to the documentation of this file.

/* Copyright (C) 2001-2012 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. */

#include <cstring>
#include <cmath>
#include "matpackI.h"
#include "exceptions.h"

using std::setw;
using std::runtime_error;


// Define the global joker object:
extern const Joker joker = Joker();

static const Numeric RAD2DEG = 57.295779513082323;

// 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. */
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);
}

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

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

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
    }
}

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 ConstVectorView:
// ------------------------------

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

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

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

  return s;
}

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

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

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

ConstVectorView::operator ConstMatrixView() const
{
    // The old version (before 2013-01-18) of this was:
    //    return ConstMatrixView(mdata,mrange,Range(mrange.mstart,1));
    // Bus this was a bug! The problem is that the matrix index operator adds
    // the mstart from both row and columm range object to mdata
    
    return ConstMatrixView(mdata,mrange,Range(0,1));
}

ConstVectorView::ConstVectorView(const Numeric& a) :
  mrange(0,1), mdata(&const_cast<Numeric&>(a))
{
  // Nothing to do here.
}

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

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

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

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 << " " << setw(3) << *i;
    }

  return os;
}


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

VectorView::VectorView (const Vector&)
{
  throw runtime_error("Creating a VectorView from a const Vector is not allowed.");
  // This is not really a runtime error, but I don't want to start
  // producing direct output from inside matpack. And just exiting is
  // not so nice. 
  // If you see this error, there is a bug in the code, not in the
  // ARTS input.
}

VectorView::VectorView (Vector& v)
{
  mdata = v.mdata;
  mrange = v.mrange;
}

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

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

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

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

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

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

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;
}

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;
}

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;
}

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

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

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

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

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

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;
}

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;
}

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;
}

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;
}

VectorView::operator MatrixView()
{
    // The old version (before 2013-01-18) of this was:
    //    return ConstMatrixView(mdata,mrange,Range(mrange.mstart,1));
    // Bus this was a bug! The problem is that the matrix index operator adds
    // the mstart from both row and columm range object to mdata
    
    return MatrixView(mdata,mrange,Range(0,1));
}

const Numeric * VectorView::get_c_array() const
{
    if (mrange.mstart != 0 || mrange.mstride != 1)
        throw std::runtime_error("A VectorView can only be converted to a plain C-array if it's pointing to a continuous block of data");

    return mdata;
}

Numeric * VectorView::get_c_array()
{
    if (mrange.mstart != 0 || mrange.mstride != 1)
        throw std::runtime_error("A VectorView can only be converted to a plain C-array if it's pointing to a continuous block of data");
    
  return mdata;
}

VectorView::VectorView(Numeric& a) :
  ConstVectorView(a)
{
  // Nothing to do here.
}

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

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

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

void copy(ConstIterator1D origin,
                 const ConstIterator1D& end,
                 Iterator1D target)
{
  if (origin.mstride == 1 && target.mstride == 1)
    memcpy((void *)target.mx,
           (void *)origin.mx,
           sizeof(Numeric) * (end.mx - origin.mx));
  else
    for ( ; origin!=end ; ++origin,++target )
      *target = *origin;
}

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


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

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

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

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:
  const Numeric *stop = mdata+n;
  for ( Numeric *x=mdata; x<stop; ++x )
    *x = fill;
}

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;
    }
}

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

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

Vector::Vector(const std::vector<Numeric>& v) :
  VectorView( new Numeric[v.size()],
              Range(0,v.size()))
{
    std::vector<Numeric>::const_iterator vec_it_end = v.end();
    Iterator1D this_it = this->begin();
    for (std::vector<Numeric>::const_iterator vec_it = v.begin();
         vec_it != vec_it_end;
         ++vec_it, ++this_it)
        *this_it = *vec_it;
}


Vector& Vector::operator=(Vector v)
{
  swap(*this, v);
  return *this;
}


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

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;
// }

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


void swap(Vector& v1, Vector& v2)
{
  std::swap(v1.mrange, v2.mrange);
  std::swap(v1.mdata, v2.mdata);
}


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


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

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

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

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

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);
}

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);
}

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

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

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

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

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.
}

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:
// -------------------------

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

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

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

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

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);
}

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);
}

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

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

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

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

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;
}

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;
}

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;
}

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;
}

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

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;
}

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;
}

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;
}

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;
}

const Numeric *MatrixView::get_c_array() const
{
    if (mrr.mstart != 0 || mrr.mstride != mcr.mextent
        || mcr.mstart != 0 || mcr.mstride != 1)
        throw std::runtime_error("A MatrixView can only be converted to a plain C-array if it's pointing to a continuous block of data");

    return mdata;
}

Numeric *MatrixView::get_c_array()
{
    if (mrr.mstart != 0 || mrr.mstride != mcr.mextent
        || mcr.mstart != 0 || mcr.mstride != 1)
        throw std::runtime_error("A MatrixView can only be converted to a plain C-array if it's pointing to a continuous block of data");

    return mdata;
}

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;
}

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;
}

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;
}

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;
}

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;
}

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;
}

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;
}

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;
}

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

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

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.
}

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;
    }
}

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:
// ---------------------

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. 
}

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

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:
  const Numeric *stop = mdata+r*c;
  for ( Numeric *x=mdata; x<stop; ++x )
    *x = fill;
}

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());
}

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());
}


Matrix& Matrix::operator=(Matrix m)
{
  swap(*this, m);
  return *this;
}

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


Matrix& Matrix::operator=(const ConstVectorView& v)
{
  resize( v.nelem(), 1 ); 
  ConstMatrixView dummy(v);
  copy( dummy.begin(), dummy.end(), begin() );
  return *this;
}

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;
    }
}


void swap(Matrix& m1, Matrix& m2)
{
  std::swap(m1.mrr, m2.mrr);
  std::swap(m1.mcr, m2.mcr);
  std::swap(m1.mdata, m2.mdata);
}


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


// Some general Matrix Vector functions:

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;
}

void mult( VectorView y,
           const ConstMatrixView& M,
           const ConstVectorView& x )
{
  // Check dimensions:
  assert( y.mrange.mextent == M.mrr.mextent );
  assert( M.mcr.mextent == x.mrange.mextent );
  assert( M.mcr.mextent != 0 && M.mrr.mextent != 0);

  // Let's first find the pointers to the starting positions
  Numeric *mdata   = M.mdata + M.mcr.mstart + M.mrr.mstart;
  Numeric *xdata   = x.mdata + x.mrange.mstart;
  Numeric *yelem   = y.mdata + y.mrange.mstart;

  Index i = M.mrr.mextent;
  while (i--)
    {
      Numeric *melem = mdata;
      Numeric *xelem = xdata; // Reset xelem to first element of source vector

      // Multiply first element of matrix row with first element of source
      // vector. This is done outside the loop to avoid initialization of the
      // target vector's element with zero (saves one assignment)
      *yelem = *melem * *xelem;

      Index j = M.mcr.mextent;   // --j (instead of j-- like in the outer loop)
      while (--j)                // is important here because we only want
        {                        // mextent-1 cycles
          melem += M.mcr.mstride;
          xelem += x.mrange.mstride;
          *yelem += *melem * *xelem;
        }

      mdata += M.mrr.mstride;     // Jump to next matrix row
      yelem += y.mrange.mstride;  // Jump to next element in target vector
    }
}


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);
        }
    }
}


void cross3(VectorView c, const ConstVectorView& a, const ConstVectorView& b)
{
  assert( a.nelem() == 3 );
  assert( b.nelem() == 3 );
  assert( c.nelem() == 3 );

  c[0] = a[1]*b[2] - a[2]*b[1];
  c[1] = a[2]*b[0] - a[0]*b[2];
  c[2] = a[0]*b[1] - a[1]*b[0];
}

Numeric vector_angle(ConstVectorView a, ConstVectorView b)
{
    assert( a.nelem() == b.nelem() );
    Numeric arg = (a*b) / sqrt(a*a) / sqrt(b*b);

    // Due to numerical inaccuracy, arg might be slightly larger than 1.
    // We catch those cases to avoid spurious returns of NaNs
    return fabs(arg) > 1. ? 0. : acos(arg) * RAD2DEG;
}


void proj(Vector& c, ConstVectorView a, ConstVectorView b)
{
    assert( a.nelem() == b.nelem() );
    assert( a.nelem() == c.nelem() );

    const Numeric C = (a*b) / (a*a);
        c  = a;
        c *= C;
};


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

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

MatrixView transpose(Vector v)
{
  return transpose((MatrixView)v);
}

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);
}

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);
    }
}

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;
}

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;
}

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;
}

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;
}

Numeric mean(const ConstVectorView& x)
{
  // Initial value for mean:
  Numeric mean = 0;

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

  for ( ; xi!=xe ; ++xi ) mean += *xi;

  mean /= (Numeric)x.nelem();
  
  return mean;
}

Numeric mean(const ConstMatrixView& x)
{
  // Initial value for mean:
  Numeric mean = 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 ) mean += *cx;
    }
  
  mean /= (Numeric)(x.nrows()*x.ncols());

  return mean;
}

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

  // Check that sizes are compatible:
  assert(mrange.mextent==v.nelem());

  // Iterators for Array:
  Array<Numeric>::const_iterator i=v.begin();
  const Array<Numeric>::const_iterator e=v.end();
  // Iterator for Vector:
  Iterator1D target = begin();

  for ( ; i!=e ; ++i,++target )
    *target = *i;

  return *this;
}

// Helper function for debugging
#ifndef NDEBUG

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

#endif