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 <Eigen/Dense>
#include <cassert>
#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)
    : 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);
    }

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

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

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

  constexpr 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, except for Joker:
    assert(mstart <= prev_fin || mextent == -1);

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

  // 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&);

  // Member functions:

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

  constexpr Range operator()(const Range r) const {
    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 { 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)
      : mx(x), mstride(stride) { /* Nothing to do here. */
  }

  // Operators:

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

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

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

#ifdef __GLIBCXX__
  bool operator==(const Iterator1D& other) const { return !operator!=(other); }

  Index operator-(const Iterator1D& other) const {
    return (Index)(mx - other.mx) / mstride;
  }
#endif

  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)
      : mx(x), mstride(stride) { /* Nothing to do here. */
  }

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

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

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

  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;

  Index nelem() const;

  Numeric sum() const;

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

  Numeric get(Index n) const {
    return *(mdata + mrange.mstart + n * mrange.mstride);
  }

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

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

  // Functions returning iterators:

  ConstIterator1D begin() const;

  ConstIterator1D end() const;

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

 protected:
  // Constructors:
  ConstVectorView() = default;

  ConstVectorView(Numeric* data, const Range& range);

  ConstVectorView(Numeric* data, const Range& p, const Range& n);

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

  // Typedef for compatibility with STL
  typedef Iterator1D iterator;

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

  Numeric& get(Index n) {
    return *(mdata + mrange.mstart + n * mrange.mstride);
  }

  VectorView operator[](const Range& r);

  // Iterators:

  Iterator1D begin();

  Iterator1D end();

  // Assignment operators:

  VectorView& operator=(const ConstVectorView& v);

  VectorView& operator=(const VectorView& v);

  VectorView& operator=(const Vector& v);

  VectorView& operator=(const Array<Numeric>& v);

  VectorView& operator=(Numeric x);

  // Other operators:

  VectorView operator*=(Numeric x);

  VectorView operator/=(Numeric x);

  VectorView operator+=(Numeric x);

  VectorView operator-=(Numeric x);

  VectorView operator*=(const ConstVectorView& x);

  VectorView operator/=(const ConstVectorView& x);

  VectorView operator+=(const ConstVectorView& x);

  VectorView operator-=(const ConstVectorView& x);

  operator MatrixView();

  const Numeric* get_c_array() const;

  Numeric* get_c_array();

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

 protected:
  // Constructors:
  VectorView() = default;

  VectorView(Numeric* data, const Range& range);

  VectorView(Numeric* data, const Range& p, const Range& n);
};

class Iterator2D {
 public:
  // Constructors:
  Iterator2D() = default;

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

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

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

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

  VectorView& operator*() { return msv; }

 private:
  VectorView msv;
  Index mstride{0};
};

class ConstIterator2D {
 public:
  // Constructors:
  ConstIterator2D() = default;

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

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

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

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

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

 private:
  ConstVectorView msv;
  Index mstride{0};
};

class Vector : public VectorView {
 public:
  // Constructors:
  Vector() = default;

  Vector(std::initializer_list<Numeric> 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>&);

  // Assignment operators:

  Vector& operator=(const Vector& v);

  Vector& operator=(Vector&& v) noexcept;

  Vector& operator=(std::initializer_list<Numeric> v);

  Vector& operator=(const Array<Numeric>& v);

  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;
  Index nrows() const;
  Index ncols() const;

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

    return get(r, c);
  }

  Numeric get(Index r, Index c) const {
    return *(mdata + mrr.mstart + r * mrr.mstride + mcr.mstart +
             c * mcr.mstride);
  }

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

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

  // View on diagonal vector
  ConstVectorView diagonal() const;

  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);
  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&);
  friend void mult_general(MatrixView,
                           const ConstMatrixView&,
                           const ConstMatrixView&);
  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);
  ConstMatrixView(Numeric* data,
                  const Range& pr,
                  const Range& pc,
                  const Range& nr,
                  const Range& nc);

  // 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) {  // Check if indices are valid:
    assert(0 <= r);
    assert(0 <= c);
    assert(r < mrr.mextent);
    assert(c < mcr.mextent);

    return get(r, c);
  }

  Numeric& get(Index r, Index c) {
    return *(mdata + mrr.mstart + r * mrr.mstride + mcr.mstart +
             c * mcr.mstride);
  }

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

  // Functions returning iterators:
  Iterator2D begin();
  Iterator2D end();

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

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

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

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

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

  virtual ~MatrixView() = 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);
  friend MatrixView transpose(MatrixView m);
  friend void mult(MatrixView, const ConstMatrixView&, const ConstMatrixView&);

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

class Matrix : public MatrixView {
 public:
  // Constructors:
  Matrix() = default;
  Matrix(Index r, Index c);
  Matrix(Index r, Index c, Numeric fill);
  Matrix(const ConstMatrixView& v);
  Matrix(const Matrix& v);
  Matrix(Matrix&& v) noexcept : MatrixView(std::forward<MatrixView>(v)) {
    v.mdata = nullptr;
  }

  // Assignment operators:
  Matrix& operator=(const Matrix& x);
  Matrix& operator=(Matrix&& x) 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();

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

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

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

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

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

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

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

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

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

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

void cross3(VectorView c, const ConstVectorView& a, const ConstVectorView& b);

Numeric vector_angle(ConstVectorView a, ConstVectorView b);

void proj(Vector& c, ConstVectorView a, ConstVectorView b);

ConstMatrixView transpose(ConstMatrixView m);

MatrixView transpose(MatrixView m);

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

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

Numeric max(const ConstVectorView& x);

Numeric max(const ConstMatrixView& x);

Numeric min(const ConstVectorView& x);

Numeric min(const ConstMatrixView& x);

Numeric mean(const ConstVectorView& x);

Numeric mean(const ConstMatrixView& x);

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

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