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
 
#include <algorithm>
#include <utility>
 
#include "array.h"
#include "matpack.h"
#include "matpack_concepts.h"
 
// Declare existance of some classes
class bofstream;
class Vector;
class Matrix;
struct Sparse;
class Tensor3;
class Tensor4;
class Tensor5;
class Tensor6;
class Tensor7;
 
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 const 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 struct 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:
 
  [[nodiscard]] constexpr Index get_start() const noexcept { return mstart; }
  [[nodiscard]] constexpr Index get_extent() const noexcept { return mextent; }
  [[nodiscard]] constexpr Index get_stride() const 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;
  };
 
  friend std::ostream& operator<<(std::ostream& os, const Range& r);
 
  constexpr void swap(Range& other) noexcept {
    using std::swap;
    swap(mstart, other.mstart);
    swap(mextent, other.mextent);
    swap(mstride, other.mstride);
  }
 
  friend constexpr void swap(Range& a, Range& b) noexcept { a.swap(b); }
 
 private:
  Index mstart;
  Index mextent;
  Index mstride;
};
 
template <size_t N>
struct Shape {
  std::array<Index, N> data;
  bool operator<=>(const Shape& other) const = default;
  friend std::ostream& operator<<(std::ostream& os, const Shape& shape) {
    os << shape.data[0];
    for (size_t i = 1; i < N; i++) os << 'x' << shape.data[i];
    return os;
  }
};
 
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;
    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;
    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:
  static constexpr bool matpack_type{true};
 
  constexpr ConstVectorView(const ConstVectorView&) = default;
  constexpr ConstVectorView(ConstVectorView&&) = default;
  ConstVectorView& operator=(const ConstVectorView&) = default;
  ConstVectorView& operator=(ConstVectorView&&) = default;
 
  // Typedef for compatibility with STL
  using const_iterator = ConstIterator1D;
 
  // Member functions:
 
  [[nodiscard]] bool empty() const noexcept { return mrange.mextent == 0; };
 
  [[nodiscard]] Index nelem() const noexcept { return mrange.mextent; }
 
  [[nodiscard]] Index size() const noexcept { return mrange.mextent; }
 
  [[nodiscard]] Shape<1> shape() const { return {nelem()}; }
 
  [[nodiscard]] 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);
  }
 
  [[nodiscard]] 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:
 
  [[nodiscard]] ConstIterator1D begin() const ARTS_NOEXCEPT;
 
  [[nodiscard]] 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 std::ostream& operator<<(std::ostream& os, const ConstVectorView& v);
 
  ConstVectorView(const Numeric& a) ARTS_NOEXCEPT;
 
  [[nodiscard]] Index selem() const noexcept {return mrange.mstart;}
 
  [[nodiscard]] Index delem() const noexcept {return mrange.mstride;}
 
  [[nodiscard]] Numeric* get_c_array() const noexcept {return mdata;}
 
 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
  using iterator = Iterator1D;
 
#define GETFUN(n) *(mdata + mrange.mstart + n * mrange.mstride)
  Numeric& operator[](Index n) ARTS_NOEXCEPT {  // Check if index is valid:
    ARTS_ASSERT(0 <= n);
    ARTS_ASSERT(n < mrange.mextent);
    return GETFUN(n);
  }
 
  Numeric& get(Index n) ARTS_NOEXCEPT { return GETFUN(n); }
#undef GETFUN
 
  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;
 
  ~VectorView() override = 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;
    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(ConstVectorView x, Index stride) ARTS_NOEXCEPT
      : msv(std::move(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;
    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);
 
  explicit Vector(const matpack::vector_like_not_vector auto &init)
      : Vector(matpack::column_size(init)) {
    *this = 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) ARTS_NOEXCEPT;
 
  Vector& operator=(std::initializer_list<Numeric> v);
 
  Vector &operator=(const matpack::vector_like_not_vector auto &init) {
    if (const auto s = matpack::shape<Index, 1>(init); shape().data not_eq s) resize(s[0]);
 
    for (Index i = 0; i < size(); i++)
      operator[](i) = init[i];
 
    return *this;
  }
 
  Vector& operator=(const Array<Numeric>& x);
 
  Vector& operator=(Numeric x);
 
  void resize(Index n);
 
  friend void swap(Vector& v1, Vector& v2) noexcept;
 
  ~Vector() noexcept override;
 
  template <class F>
  void transform_elementwise(F&& func) {
    std::transform(mdata, mdata + size(), mdata, func);
  }
};
 
// Declare class Matrix:
class Matrix;
 
class ConstMatrixView {
 public:
  static constexpr bool matpack_type{true};
 
  constexpr ConstMatrixView(const ConstMatrixView&) = default;
  constexpr ConstMatrixView(ConstMatrixView&&) = default;
  ConstMatrixView& operator=(const ConstMatrixView&) = default;
  ConstMatrixView& operator=(ConstMatrixView&&) = default;
 
  // Typedef for compatibility with STL
  using const_iterator = ConstIterator2D;
 
  // Member functions:
  [[nodiscard]] Index selem() const noexcept { return mrr.mstart + mcr.mstart; }
  [[nodiscard]] Index nrows() const noexcept { return mrr.mextent; }
  [[nodiscard]] Index ncols() const noexcept { return mcr.mextent; }
  [[nodiscard]] Index drows() const noexcept { return mrr.mstride; }
  [[nodiscard]] Index dcols() const noexcept { return mcr.mstride; }
 
  // Total size
  [[nodiscard]] Index size() const noexcept { return nrows() * ncols(); }
  [[nodiscard]] bool empty() const noexcept { return size() == 0; }
 
  [[nodiscard]] Shape<2> shape() const { return {nrows(), ncols()}; }
 
  // 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);
  }
 
  [[nodiscard]] 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:
  [[nodiscard]] ConstIterator2D begin() const ARTS_NOEXCEPT;
  [[nodiscard]] ConstIterator2D end() const ARTS_NOEXCEPT;
 
  // View on diagonal vector
  [[nodiscard]] 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 std::ostream& operator<<(std::ostream& os, const ConstMatrixView& v);
 
  // Conversion to a plain C-array
  [[nodiscard]] Numeric* get_c_array() const noexcept {return mdata;}
 
 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
  using iterator = Iterator2D;
 
#define GETFUN(r, c) \
  *(mdata + mrr.mstart + r * mrr.mstride + mcr.mstart + c * mcr.mstride)
  // 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 GETFUN(r, c);
  }
 
  Numeric& get(Index r, Index c) ARTS_NOEXCEPT { return GETFUN(r, c); }
#undef GETFUN
 
  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;
 
  ~MatrixView() override = 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;
  }
 
  explicit Matrix(const matpack::matrix_like_not_matrix auto &init)
      : Matrix(matpack::row_size(init), matpack::column_size(init)) {
    *this = init;
  }
 
  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) ARTS_NOEXCEPT;
  Matrix& operator=(Numeric x);
  Matrix& operator=(const ConstVectorView& v);
 
  Matrix &operator=(const matpack::matrix_like_not_matrix auto &init) {
    if (const auto s = matpack::shape<Index, 2>(init); 
    shape().data not_eq s)
      resize(s[0], s[1]);
 
    auto [I, J] = shape().data;
    for (Index i = 0; i < I; i++)
      for (Index j = 0; j < J; j++)
        operator()(i, j) = init(i, j);
 
    return *this;
  }
 
  // Resize function:
  void resize(Index r, Index c);
 
  // Swap function:
  friend void swap(Matrix& m1, Matrix& m2) noexcept;
 
  // Destructor:
  ~Matrix() noexcept override;
 
  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; }
 
  template <class F>
  void transform_elementwise(F&& func) {
    std::transform(mdata, mdata + size(), mdata, func);
  }
};
 
// 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;
 
// Helper function for debugging
#ifndef NDEBUG
 
Numeric debug_matrixview_get_elem(MatrixView& mv, Index r, Index c);
 
#endif
 
using ArrayOfVector = Array<Vector>;
 
using ArrayOfArrayOfVector = Array<ArrayOfVector>;
 
using ArrayOfMatrix = Array<Matrix>;
 
using ArrayOfArrayOfMatrix = Array<ArrayOfMatrix>;
 
#endif  // matpackI_h