matpackI.h

Go to the documentation of this file.

/* Copyright (C) 2001-2019 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
 
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wshadow"
#pragma GCC diagnostic ignored "-Wconversion"
#pragma GCC diagnostic ignored "-Wfloat-conversion"
#include <Eigen/Dense>
#pragma GCC diagnostic pop
 
#include "array.h"
#include "matpack.h"
 
// Declare existance of some classes
class bofstream;
 
// Declaration of Eigen types
typedef Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic> StrideType;
typedef Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
    MatrixType;
typedef Eigen::Map<MatrixType, 0, StrideType> MatrixViewMap;
typedef Eigen::Map<const MatrixType, 0, StrideType> ConstMatrixViewMap;
typedef Eigen::Matrix<Numeric, 4, 4, Eigen::RowMajor> Matrix4x4Type;
typedef Eigen::Map<Matrix4x4Type, 0, StrideType> Matrix4x4ViewMap;
typedef Eigen::Map<const Matrix4x4Type, 0, StrideType> ConstMatrix4x4ViewMap;
 
class Joker {
  // Nothing here.
};
 
// Declare existence of the global joker object:
extern const Joker joker;
 
// Declare the existence of class ConstMatrixView:
class ConstIterator1D;
 
// 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:
 
  constexpr Range(Index start, Index extent, Index stride = 1) ARTS_NOEXCEPT
    : mstart(start), mextent(extent), mstride(stride) {
      // Start must be >= 0:
      ARTS_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.
      ARTS_ASSERT(0 <= mextent);
      // Stride can be anything except 0.
      // SAB 2001-09-21: Allow 0 stride.
      //  ARTS_ASSERT( 0!=mstride);
    }
 
  constexpr Range(Index start, Joker, Index stride = 1) ARTS_NOEXCEPT
  : mstart(start), mextent(-1), mstride(stride) {
    // Start must be >= 0:
    ARTS_ASSERT(0 <= mstart);
  }
 
  constexpr Range(Joker, Index stride = 1) noexcept : mstart(0), mextent(-1), mstride(stride) {
    // Nothing to do here.
  };
 
  constexpr Range(Index max_size, const Range& r) ARTS_NOEXCEPT
  : mstart(r.mstart), mextent(r.mextent), mstride(r.mstride) {
    // Start must be >= 0:
    ARTS_ASSERT(0 <= mstart);
    // ... and < max_size:
    ARTS_ASSERT(mstart < max_size);
    
    // Stride must be != 0:
    ARTS_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;
      ARTS_ASSERT(0 <= fin);
      ARTS_ASSERT(fin < max_size);
      #endif
    }
  }
 
  constexpr Range(const Range& p, const Range& n) ARTS_NOEXCEPT
      : 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:
    ARTS_ASSERT(p.mstart <= mstart);
    // and <= prev_fin, except for Joker:
    ARTS_ASSERT(mstart <= prev_fin || mextent == -1);
 
    // Resulting stride must be != 0:
    ARTS_ASSERT(0 != mstride, "Problem: mstride==", mstride, " for mstart==",mstart, " and mextent==", mextent);
 
    // 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;
      ARTS_ASSERT(p.mstart <= fin);
      ARTS_ASSERT(fin <= prev_fin);
  #endif
    }
  };
 
  // 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 class ConstComplexVectorView;
  friend class ComplexVectorView;
  friend class ComplexVector;
  friend class ConstComplexMatrixView;
  friend class ComplexMatrixView;
  friend class ComplexMatrix;
  friend class ComplexIterator2D;
  friend class ConstComplexIterator2D;
 
  friend void mult_general(VectorView,
                           const ConstMatrixView&,
                           const ConstVectorView&) ARTS_NOEXCEPT;
 
  // Member functions:
 
  constexpr Index get_start() const ARTS_NOEXCEPT { return mstart; }
  constexpr Index get_extent() const ARTS_NOEXCEPT { return mextent; }
  constexpr Index get_stride() const ARTS_NOEXCEPT { return mstride; }
 
  constexpr Range operator()(const Range r) const ARTS_NOEXCEPT {
    return (r.mextent < 0) ? (mextent < 0) ? Range(mstart + r.mstart * mstride,
                                                   joker,
                                                   r.mstride * mstride)
                                           : Range(mstart + r.mstart * mstride,
                                                   mextent,
                                                   r.mstride * mstride)
                           : Range(mstart + r.mstart * mstride,
                                   r.mextent,
                                   r.mstride * mstride);
  }
 
  constexpr Index operator()(const Index i) const ARTS_NOEXCEPT { return mstart + i * mstride; };
 
 private:
  Index mstart;
  Index mextent;
  Index mstride;
};
 
class Iterator1D {
 public:
  using difference_type = Index;
  using value_type = Numeric;
  using pointer = Numeric*;
  using reference = Numeric&;
  using iterator_category = std::random_access_iterator_tag;
 
  Iterator1D() = default;
 
  Iterator1D(Numeric* x, Index stride) ARTS_NOEXCEPT
      : mx(x), mstride(stride) { /* Nothing to do here. */
  }
 
  // Operators:
 
  Iterator1D& operator++() ARTS_NOEXCEPT {
    mx += mstride;
    return *this;
  }
 
  Numeric& operator*() const ARTS_NOEXCEPT { return *mx; }
 
  bool operator!=(const Iterator1D& other) const ARTS_NOEXCEPT {
    if (mx != other.mx)
      return true;
    else
      return false;
  }
 
  bool operator==(const Iterator1D& other) const ARTS_NOEXCEPT {
    return !operator!=(other);
  }
 
  Index operator-(const Iterator1D& other) const ARTS_NOEXCEPT {
    return (Index)(mx - other.mx) / mstride;
  }
 
  friend void copy(ConstIterator1D origin,
                   const ConstIterator1D& end,
                   Iterator1D target);
 
 private:
  Numeric* mx{nullptr};
  Index mstride{0};
};
 
class ConstIterator1D {
 public:
  using difference_type = Index;
  using value_type = const Numeric;
  using pointer = const Numeric*;
  using reference = const Numeric&;
  using iterator_category = std::random_access_iterator_tag;
 
  ConstIterator1D() = default;
 
  ConstIterator1D(const Numeric* x, Index stride) ARTS_NOEXCEPT
      : mx(x), mstride(stride) { /* Nothing to do here. */
  }
 
  // Operators:
  ConstIterator1D& operator++() ARTS_NOEXCEPT {
    mx += mstride;
    return *this;
  }
 
  const Numeric& operator*() const ARTS_NOEXCEPT { return *mx; }
 
  bool operator!=(const ConstIterator1D& other) const ARTS_NOEXCEPT {
    if (mx != other.mx)
      return true;
    else
      return false;
  }
 
  bool operator==(const ConstIterator1D& other) const ARTS_NOEXCEPT {
    return !operator!=(other);
  }
 
  Index operator-(const ConstIterator1D& other) const ARTS_NOEXCEPT {
    return (Index)(mx - other.mx) / mstride;
  }
 
  friend void copy(ConstIterator1D origin,
                   const ConstIterator1D& end,
                   Iterator1D target);
 
 private:
  const Numeric* mx{nullptr};
  Index mstride{0};
};
 
// Declare the vector class:
class Vector;
 
// Declare the MatrixView class
class MatrixView;
 
class ConstVectorView {
 public:
  constexpr ConstVectorView(const ConstVectorView&) = default;
  constexpr ConstVectorView(ConstVectorView&&) = default;
  ConstVectorView& operator=(const ConstVectorView&) = default;
  ConstVectorView& operator=(ConstVectorView&&) = default;
 
  // Typedef for compatibility with STL
  typedef ConstIterator1D const_iterator;
 
  // Member functions:
 
  bool empty() const ARTS_NOEXCEPT;
 
  Index nelem() const ARTS_NOEXCEPT;
  
  Index size() const ARTS_NOEXCEPT;
 
  Numeric sum() const ARTS_NOEXCEPT;
 
  // Const index operators:
  Numeric operator[](Index n) const ARTS_NOEXCEPT {  // Check if index is valid:
    ARTS_ASSERT(0 <= n);
    ARTS_ASSERT(n < mrange.mextent);
    return get(n);
  }
 
  Numeric get(Index n) const ARTS_NOEXCEPT {
    return *(mdata + mrange.mstart + n * mrange.mstride);
  }
 
  ConstVectorView operator[](const Range& r) const ARTS_NOEXCEPT;
 
  friend Numeric operator*(const ConstVectorView& a, const ConstVectorView& b) ARTS_NOEXCEPT;
 
  // Functions returning iterators:
 
  ConstIterator1D begin() const ARTS_NOEXCEPT;
 
  ConstIterator1D end() const ARTS_NOEXCEPT;
 
  operator ConstMatrixView() const;
 
  virtual ~ConstVectorView() = default;
 
  // 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 class ConstComplexVectorView;
  friend int poly_root_solve(Matrix& roots, Vector& coeffs);
  friend void mult(VectorView, const ConstMatrixView&, const ConstVectorView&);
  friend void mult(VectorView, const Sparse&, ConstVectorView);
  friend void transpose_mult(VectorView, const Sparse&, ConstVectorView);
  friend void mult_general(VectorView,
                           const ConstMatrixView&,
                           const ConstVectorView&) ARTS_NOEXCEPT;
  friend void lubacksub(VectorView,
                        ConstMatrixView,
                        ConstVectorView,
                        const ArrayOfIndex&);
  friend void diagonalize(MatrixView, VectorView, VectorView, ConstMatrixView);
 
  friend ConstMatrixViewMap MapToEigen(const ConstVectorView&);
  friend ConstMatrixViewMap MapToEigenCol(const ConstVectorView&);
  friend MatrixViewMap MapToEigen(VectorView&);
  friend MatrixViewMap MapToEigenCol(VectorView&);
 
  ConstVectorView(const Numeric& a) ARTS_NOEXCEPT;
 
 protected:
  // Constructors:
  ConstVectorView() = default;
 
  ConstVectorView(Numeric* data, const Range& range) ARTS_NOEXCEPT;
 
  ConstVectorView(Numeric* data, const Range& p, const Range& n) ARTS_NOEXCEPT;
 
  // Data members:
  // -------------
  Range mrange{0, 0};
  Numeric* mdata{nullptr};
};
 
class VectorView : public ConstVectorView {
 public:
  // Make const methods visible from base class
  using ConstVectorView::begin;
  using ConstVectorView::end;
  using ConstVectorView::operator[];
  using ConstVectorView::get;
 
  constexpr VectorView(const VectorView&) = default;
 
  VectorView(const Vector&);
 
  VectorView(Vector& v) ARTS_NOEXCEPT;
 
  // Typedef for compatibility with STL
  typedef Iterator1D iterator;
 
  Numeric& operator[](Index n) ARTS_NOEXCEPT {  // Check if index is valid:
    ARTS_ASSERT(0 <= n);
    ARTS_ASSERT(n < mrange.mextent);
    return get(n);
  }
 
  Numeric& get(Index n) ARTS_NOEXCEPT {
    return *(mdata + mrange.mstart + n * mrange.mstride);
  }
 
  VectorView operator[](const Range& r) ARTS_NOEXCEPT;
 
  // Iterators:
 
  Iterator1D begin() ARTS_NOEXCEPT;
 
  Iterator1D end() ARTS_NOEXCEPT;
 
  // 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) ARTS_NOEXCEPT;
 
  VectorView operator/=(Numeric x) ARTS_NOEXCEPT;
 
  VectorView operator+=(Numeric x) ARTS_NOEXCEPT;
 
  VectorView operator-=(Numeric x) ARTS_NOEXCEPT;
 
  VectorView operator*=(const ConstVectorView& x) ARTS_NOEXCEPT;
 
  VectorView operator/=(const ConstVectorView& x) ARTS_NOEXCEPT;
 
  VectorView operator+=(const ConstVectorView& x) ARTS_NOEXCEPT;
 
  VectorView operator-=(const ConstVectorView& x) ARTS_NOEXCEPT;
 
  operator MatrixView() ARTS_NOEXCEPT;
 
  const Numeric* get_c_array() const ARTS_NOEXCEPT;
 
  Numeric* get_c_array() ARTS_NOEXCEPT;
 
  virtual ~VectorView() = default;
 
  // 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;
  friend class ComplexVectorView;
 
  VectorView(Numeric& a) ARTS_NOEXCEPT;
 
 protected:
  // Constructors:
  VectorView() = default;
 
  VectorView(Numeric* data, const Range& range) ARTS_NOEXCEPT;
 
  VectorView(Numeric* data, const Range& p, const Range& n) ARTS_NOEXCEPT;
};
 
class Iterator2D {
 public:
  // Constructors:
  Iterator2D() = default;
 
  Iterator2D(const VectorView& x, Index stride) ARTS_NOEXCEPT
      : msv(x), mstride(stride) { /* Nothing to do here. */
  }
 
  // Operators:
  Iterator2D& operator++() ARTS_NOEXCEPT {
    msv.mdata += mstride;
    return *this;
  }
 
  bool operator!=(const Iterator2D& other) const ARTS_NOEXCEPT {
    if (msv.mdata + msv.mrange.mstart !=
        other.msv.mdata + other.msv.mrange.mstart)
      return true;
    else
      return false;
  }
 
  VectorView* operator->() ARTS_NOEXCEPT { return &msv; }
 
  VectorView& operator*() ARTS_NOEXCEPT { return msv; }
 
 private:
  VectorView msv;
  Index mstride{0};
};
 
class ConstIterator2D {
 public:
  // Constructors:
  ConstIterator2D() = default;
 
  ConstIterator2D(const ConstVectorView& x, Index stride) ARTS_NOEXCEPT
      : msv(x), mstride(stride) { /* Nothing to do here. */
  }
 
  // Operators:
  ConstIterator2D& operator++() ARTS_NOEXCEPT {
    msv.mdata += mstride;
    return *this;
  }
 
  bool operator!=(const ConstIterator2D& other) const ARTS_NOEXCEPT {
    if (msv.mdata + msv.mrange.mstart !=
        other.msv.mdata + other.msv.mrange.mstart)
      return true;
    else
      return false;
  }
 
  const ConstVectorView* operator->() const ARTS_NOEXCEPT { return &msv; }
 
  const ConstVectorView& operator*() const ARTS_NOEXCEPT { return msv; }
 
 private:
  ConstVectorView msv;
  Index mstride{0};
};
 
class Vector : public VectorView {
 public:
  // Constructors:
  Vector() = default;
 
  Vector(std::initializer_list<Numeric> init);
  
  Vector(const Eigen::VectorXd& init);
 
  explicit Vector(Index n);
 
  Vector(Index n, Numeric fill);
 
  Vector(Numeric start, Index extent, Numeric stride);
 
  Vector(const ConstVectorView& v);
 
  Vector(const Vector& v);
  Vector(Vector&& v) noexcept : VectorView(std::forward<VectorView>(v)) {
    v.mdata = nullptr;
  }
 
  Vector(const std::vector<Numeric>&);
  
  Vector(Numeric* d, const Range& r0) ARTS_NOEXCEPT
  : VectorView(d, r0) {
    ARTS_ASSERT(r0.get_extent() >= 0, "Must have size. Has: ", r0.get_extent());
  }
 
  // Assignment operators:
 
  Vector& operator=(const Vector& v);
 
  Vector& operator=(Vector&& v) noexcept;
 
  Vector& operator=(std::initializer_list<Numeric> v);
 
  Vector& operator=(const Array<Numeric>& x);
 
  Vector& operator=(Numeric x);
 
  void resize(Index n);
 
  friend void swap(Vector& v1, Vector& v2);
 
  virtual ~Vector();
};
 
// Declare class Matrix:
class Matrix;
 
class ConstMatrixView {
 public:
  constexpr ConstMatrixView(const ConstMatrixView&) = default;
  constexpr ConstMatrixView(ConstMatrixView&&) = default;
  ConstMatrixView& operator=(const ConstMatrixView&) = default;
  ConstMatrixView& operator=(ConstMatrixView&&) = default;
 
  // Typedef for compatibility with STL
  typedef ConstIterator2D const_iterator;
 
  // Member functions:
  bool empty() const ARTS_NOEXCEPT;
  Index nrows() const ARTS_NOEXCEPT;
  Index ncols() const ARTS_NOEXCEPT;
 
  // Const index operators:
  Numeric operator()(Index r, Index c) const ARTS_NOEXCEPT {  // Check if indices are valid:
    ARTS_ASSERT(0 <= r);
    ARTS_ASSERT(0 <= c);
    ARTS_ASSERT(r < mrr.mextent);
    ARTS_ASSERT(c < mcr.mextent);
 
    return get(r, c);
  }
 
  Numeric get(Index r, Index c) const ARTS_NOEXCEPT {
    return *(mdata + mrr.mstart + r * mrr.mstride + mcr.mstart +
             c * mcr.mstride);
  }
 
  ConstMatrixView operator()(const Range& r, const Range& c) const ARTS_NOEXCEPT;
  ConstVectorView operator()(const Range& r, Index c) const ARTS_NOEXCEPT;
  ConstVectorView operator()(Index r, const Range& c) const ARTS_NOEXCEPT;
 
  // Functions returning iterators:
  ConstIterator2D begin() const ARTS_NOEXCEPT;
  ConstIterator2D end() const ARTS_NOEXCEPT;
 
  // View on diagonal vector
  ConstVectorView diagonal() const ARTS_NOEXCEPT;
 
  virtual ~ConstMatrixView() = default;
 
  // 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 class ConstComplexMatrixView;
  friend ConstMatrixView transpose(ConstMatrixView m) ARTS_NOEXCEPT;
  friend int poly_root_solve(Matrix& roots, Vector& coeffs);
  friend void mult(VectorView, const ConstMatrixView&, const ConstVectorView&);
  friend void mult(MatrixView, const ConstMatrixView&, const ConstMatrixView&);
  friend void mult(MatrixView, const Sparse&, const ConstMatrixView&);
  friend void mult(MatrixView, const ConstMatrixView&, const Sparse&);
  friend void mult_general(VectorView,
                           const ConstMatrixView&,
                           const ConstVectorView&) ARTS_NOEXCEPT;
  friend void mult_general(MatrixView,
                           const ConstMatrixView&,
                           const ConstMatrixView&) ARTS_NOEXCEPT;
  friend void ludcmp(Matrix&, ArrayOfIndex&, ConstMatrixView);
  friend void lubacksub(VectorView,
                        ConstMatrixView,
                        ConstVectorView,
                        const ArrayOfIndex&);
  friend void inv(MatrixView, ConstMatrixView);
  friend void diagonalize(MatrixView, VectorView, VectorView, ConstMatrixView);
 
  friend ConstMatrixViewMap MapToEigen(const ConstMatrixView&);
  friend MatrixViewMap MapToEigen(MatrixView&);
 
  friend ConstMatrix4x4ViewMap MapToEigen4x4(const ConstMatrixView&);
  friend Matrix4x4ViewMap MapToEigen4x4(MatrixView&);
 
 protected:
  // Constructors:
  ConstMatrixView() = default;
  ConstMatrixView(Numeric* data, const Range& r, const Range& c) ARTS_NOEXCEPT;
  ConstMatrixView(Numeric* data,
                  const Range& pr,
                  const Range& pc,
                  const Range& nr,
                  const Range& nc) ARTS_NOEXCEPT;
 
  // Data members:
  // -------------
  Range mrr{0, 0, 1};
  Range mcr{0, 0, 1};
  Numeric* mdata{nullptr};
};
 
class MatrixView : public ConstMatrixView {
 public:
  // Make const methods visible from base class
  using ConstMatrixView::begin;
  using ConstMatrixView::end;
  using ConstMatrixView::operator();
  using ConstMatrixView::get;
 
  constexpr MatrixView(const MatrixView&) = default;
 
  // Typedef for compatibility with STL
  typedef Iterator2D iterator;
 
  // Index Operators:
  Numeric& operator()(Index r, Index c) ARTS_NOEXCEPT {  // Check if indices are valid:
    ARTS_ASSERT(0 <= r);
    ARTS_ASSERT(0 <= c);
    ARTS_ASSERT(r < mrr.mextent);
    ARTS_ASSERT(c < mcr.mextent);
 
    return get(r, c);
  }
 
  Numeric& get(Index r, Index c) ARTS_NOEXCEPT {
    return *(mdata + mrr.mstart + r * mrr.mstride + mcr.mstart +
             c * mcr.mstride);
  }
 
  MatrixView operator()(const Range& r, const Range& c) ARTS_NOEXCEPT;
  VectorView operator()(const Range& r, Index c) ARTS_NOEXCEPT;
  VectorView operator()(Index r, const Range& c) ARTS_NOEXCEPT;
 
  // Functions returning iterators:
  Iterator2D begin() ARTS_NOEXCEPT;
  Iterator2D end() ARTS_NOEXCEPT;
 
  // Assignment operators:
  MatrixView& operator=(const ConstMatrixView& m);
  MatrixView& operator=(const MatrixView& m);
  MatrixView& operator=(const Matrix& m);
  MatrixView& operator=(const ConstVectorView& m);
  MatrixView& operator=(Numeric x);
 
  // Other operators:
  MatrixView& operator*=(Numeric x) ARTS_NOEXCEPT;
  MatrixView& operator/=(Numeric x) ARTS_NOEXCEPT;
  MatrixView& operator+=(Numeric x) ARTS_NOEXCEPT;
  MatrixView& operator-=(Numeric x) ARTS_NOEXCEPT;
 
  MatrixView& operator*=(const ConstMatrixView& x) ARTS_NOEXCEPT;
  MatrixView& operator/=(const ConstMatrixView& x) ARTS_NOEXCEPT;
  MatrixView& operator+=(const ConstMatrixView& x) ARTS_NOEXCEPT;
  MatrixView& operator-=(const ConstMatrixView& x) ARTS_NOEXCEPT;
 
  MatrixView& operator*=(const ConstVectorView& x) ARTS_NOEXCEPT;
  MatrixView& operator/=(const ConstVectorView& x) ARTS_NOEXCEPT;
  MatrixView& operator+=(const ConstVectorView& x) ARTS_NOEXCEPT;
  MatrixView& operator-=(const ConstVectorView& x) ARTS_NOEXCEPT;
 
  // Conversion to a plain C-array
  const Numeric* get_c_array() const ARTS_NOEXCEPT;
  Numeric* get_c_array() ARTS_NOEXCEPT;
 
  virtual ~MatrixView() = default;
 
  // Friends:
  friend class VectorView;
  friend class Iterator3D;
  friend class Tensor3View;
  friend class Tensor4View;
  friend class Tensor5View;
  friend class Tensor6View;
  friend class Tensor7View;
  friend class ComplexMatrixView;
  friend ConstMatrixView transpose(ConstMatrixView m) ARTS_NOEXCEPT;
  friend MatrixView transpose(MatrixView m) ARTS_NOEXCEPT;
  friend void mult(MatrixView, const ConstMatrixView&, const ConstMatrixView&);
 
 protected:
  // Constructors:
  MatrixView() = default;
  MatrixView(Numeric* data, const Range& rr, const Range& cr) ARTS_NOEXCEPT;
  MatrixView(Numeric* data,
             const Range& pr,
             const Range& pc,
             const Range& nr,
             const Range& nc) ARTS_NOEXCEPT;
};
 
class Matrix : public MatrixView {
 public:
  // Constructors:
  Matrix() = default;
  Matrix(Index r, Index c);
  Matrix(Index r, Index c, Numeric fill);
  Matrix(const ConstMatrixView& m);
  Matrix(const Matrix& m);
  Matrix(Matrix&& v) noexcept : MatrixView(std::forward<MatrixView>(v)) {
    v.mdata = nullptr;
  }
  
  Matrix(Numeric* d, const Range& r0, const Range& r1) ARTS_NOEXCEPT
  : MatrixView(d, r0, r1) {
    ARTS_ASSERT(r0.get_extent() >= 0, "Must have size. Has: ", r0.get_extent());
    ARTS_ASSERT(r1.get_extent() >= 0, "Must have size. Has: ", r1.get_extent());
  }
 
  // Assignment operators:
  Matrix& operator=(const Matrix& m);
  Matrix& operator=(Matrix&& m) noexcept;
  Matrix& operator=(Numeric x);
  Matrix& operator=(const ConstVectorView& v);
 
  // Resize function:
  void resize(Index r, Index c);
 
  // Swap function:
  friend void swap(Matrix& m1, Matrix& m2);
 
  // Destructor:
  virtual ~Matrix();
  
  // Total size
  Index size() const noexcept {return nrows() * ncols();}
  
  template <std::size_t dim0>
  Vector reduce_rank() && ARTS_NOEXCEPT {
    static_assert(dim0 < 2, "Bad Dimension, Out-of-Bounds");
    
    Range r0(0, dim0 == 0 ? nrows() : ncols());
    
    Vector out(mdata, r0);
    ARTS_ASSERT(size() == out.size(), "Can only reduce size on same size input. "
      "The sizes are (input v. output): ", size(), " v. ", out.size());
    mdata = nullptr;
    return out;
  }
 
  Numeric* get_raw_data() ARTS_NOEXCEPT { return mdata; }
};
 
// Function declarations:
// ----------------------
 
void copy(ConstIterator1D origin,
          const ConstIterator1D& end,
          Iterator1D target);
 
void copy(Numeric x, Iterator1D target, const Iterator1D& end) ARTS_NOEXCEPT;
 
void copy(ConstIterator2D origin,
          const ConstIterator2D& end,
          Iterator2D target);
 
void copy(Numeric x, Iterator2D target, const Iterator2D& end) ARTS_NOEXCEPT;
 
void mult(VectorView y, const ConstMatrixView& M, const ConstVectorView& x);
 
void mult_general(MatrixView A,
                  const ConstMatrixView& B,
                  const ConstMatrixView& C) ARTS_NOEXCEPT;
 
void mult(MatrixView A, const ConstMatrixView& B, const ConstMatrixView& C);
 
void mult_general(MatrixView A,
                  const ConstMatrixView& B,
                  const ConstMatrixView& C) ARTS_NOEXCEPT;
 
void cross3(VectorView c, const ConstVectorView& a, const ConstVectorView& b) ARTS_NOEXCEPT;
 
Numeric vector_angle(ConstVectorView a, ConstVectorView b);
 
void proj(Vector& c, ConstVectorView a, ConstVectorView b) ARTS_NOEXCEPT;
 
ConstMatrixView transpose(ConstMatrixView m) ARTS_NOEXCEPT;
 
MatrixView transpose(MatrixView m) ARTS_NOEXCEPT;
 
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) ARTS_NOEXCEPT;
 
Numeric max(const ConstMatrixView& x) ARTS_NOEXCEPT;
 
Numeric min(const ConstVectorView& x) ARTS_NOEXCEPT;
 
Numeric min(const ConstMatrixView& x) ARTS_NOEXCEPT;
 
Numeric mean(const ConstVectorView& x) ARTS_NOEXCEPT;
 
Numeric nanmean(const ConstVectorView& x) ARTS_NOEXCEPT;
 
Numeric mean(const ConstMatrixView& x) ARTS_NOEXCEPT;
 
Numeric operator*(const ConstVectorView& a, const ConstVectorView& b) ARTS_NOEXCEPT;
 
std::ostream& operator<<(std::ostream& os, const ConstVectorView& v);
 
std::ostream& operator<<(std::ostream& os, const ConstMatrixView& v);
 
std::ostream& operator<<(std::ostream& os, const Range& r);
 
// Converts constant matrix to constant eigen map
ConstMatrixViewMap MapToEigen(const ConstMatrixView& A);
// Converts constant vector to constant eigen row-view
ConstMatrix4x4ViewMap MapToEigen4x4(const ConstMatrixView& A);
// Converts constant vector to constant eigen row-view
ConstMatrixViewMap MapToEigen(const ConstVectorView& A);
// Converts constant vector to constant eigen row-view
ConstMatrixViewMap MapToEigenRow(const ConstVectorView& A);
// Converts constant vector to constant eigen column-view
ConstMatrixViewMap MapToEigenCol(const ConstVectorView& A);
// Converts matrix to eigen map
MatrixViewMap MapToEigen(MatrixView& A);
// Converts vector to eigen map row-view
Matrix4x4ViewMap MapToEigen4x4(MatrixView& A);
// Converts vector to eigen map row-view
MatrixViewMap MapToEigen(VectorView& A);
// Converts vector to eigen map row-view
MatrixViewMap MapToEigenRow(VectorView& A);
// Converts vector to eigen map column-view
MatrixViewMap MapToEigenCol(VectorView& A);
 
// Helper function for debugging
#ifndef NDEBUG
 
Numeric debug_matrixview_get_elem(MatrixView& mv, Index r, Index c);
 
#endif
 
#endif  // matpackI_h