b2matrix_rectangle.H

//------------------------------------------------------------------------
// b2matrix_rectangle.H --
//
//     A framework for generic linear algebraic (matrix) computations.
//
// written by Mathias Doreille
//            Thomas Blome <thomas.blome@dlr.de>
//
// Copyright (c) 2004-2012,2017 SMR Engineering & Development SA
//                              2502 Bienne, Switzerland
//
// Copyright (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
//                    Linder Höhe, 51147 Köln
//
// All Rights Reserved.  Proprietary source code.  The contents of
// this file may not be disclosed to third parties, copied or
// duplicated in any form, in whole or in part, without the prior
// written permission of SMR.
//------------------------------------------------------------------------
 
#ifndef __B2MATRIX_RECTANGLE_H__
#define __B2MATRIX_RECTANGLE_H__
 
#include <algorithm>
#include <vector>
 
#include "b2linear_algebra_def.H"
 
namespace b2000::b2linalg {
 
struct Mrectangle {
    using base = Mrectangle_increment_ref;
    using const_base = Mrectangle_increment_st_constref;
    using copy = Mrectangle;
};
 
template <typename T>
class Matrix<T, Mrectangle> {
public:
    Matrix() = default;
 
    explicit Matrix(size_t s1_, size_t s2_) : s1(s1_), s2(s2_), m(s1_ * s2_) {}
 
    Matrix(size_t s1_, size_t s2_, const T& t) : s1(s1_), s2(s2_), m(s1_ * s2_, t) {}
 
    Matrix(const Matrix& m_) : s1(m_.s1), s2(m_.s2), m(m_.m) {}
 
    template <typename T1>
    Matrix(const Matrix<T1, Mrectangle>& m_)
        : s1(m_.s1.begin(), m_.s1.end()),
          s2(m_.s2.begin(), m_.s2.end()),
          m(m_.m.begin(), m_.m.end()) {}
 
    template <typename T1>
    Matrix& operator=(const Matrix<T1, Mrectangle>& m_) {
        s1 = m_.s1;
        s2 = m_.s2;
        m.assign(m_.m.begin(), m_.m.end());
        return *this;
    }
 
    template <typename STORAGE>
    Matrix(const Matrix<T, STORAGE>& m_)
        : s1(m_.size1()), s2(m_.size2()), m(m_.size1() * m_.size2()) {
        Matrix<T, Mrectangle_increment_ref>(*this) = m_;
    }
 
    template <typename T1, typename STORAGE>
    Matrix& operator=(const Matrix<T1, STORAGE>& m_) {
        s1 = m_.size1();
        s2 = m_.size2();
        m.resize(m_.size1() * m_.size2());
        for (size_t j = 0; j != size2(); ++j) {
            for (size_t i = 0; i != size1(); ++i) { (*this)(i, j) = m_(i, j); }
        }
        return *this;
    }
 
    double norm_inf() const {
        double res = 0;
        for (typename std::vector<T>::const_iterator i = m.begin(); i != m.end(); ++i) {
            res = std::max(res, norm(*i));
        }
        return res;
    }
 
    void copy_real(const Matrix<b2000::csda<T>, Mrectangle_constref>& m_) {
        s1 = m_.size1();
        s2 = m_.size2();
        m.resize(m_.size1() * m_.size2());
        for (size_t i = 0; i != m.size(); ++i) { m[i] = std::real(m_.m[i]); }
    }
 
    void copy_real(const Matrix<std::complex<T>, Mrectangle_constref>& m_) {
        s1 = m_.size1();
        s2 = m_.size2();
        m.resize(m_.size1() * m_.size2());
        for (size_t i = 0; i != m.size(); ++i) { m[i] = std::real(m_.m[i]); }
    }
 
    void copy_imag(const Matrix<b2000::csda<T>, Mrectangle_constref>& m_) {
        s1 = m_.size1();
        s2 = m_.size2();
        m.resize(m_.size1() * m_.size2());
        for (size_t i = 0; i != m.size(); ++i) { m[i] = std::imag(m_.m[i]); }
    }
 
    void copy_imag(const Matrix<std::complex<T>, Mrectangle_constref>& m_) {
        s1 = m_.size1();
        s2 = m_.size2();
        m.resize(m_.size1() * m_.size2());
        for (size_t i = 0; i != m.size(); ++i) { m[i] = std::imag(m_.m[i]); }
    }
 
    bool is_null() const { return this == &null; }
 
    void set_zero() { std::fill(m.begin(), m.end(), 0); }
 
    void swap(Matrix& m_) {
        std::swap(s1, m_.s1);
        std::swap(s2, m_.s2);
        m.swap(m_.m);
    }
 
    Matrix& operator=(const Matrix& m_) {
        s1 = m_.s1;
        s2 = m_.s2;
        m = m_.m;
        return *this;
    }
 
    template <typename STORAGE>
    Matrix& operator=(const Matrix<T, STORAGE>& m_) {
        resize(m_.size());
        Matrix<T, Mrectangle_increment_ref>(*this) = m_;
        return *this;
    }
 
    template <typename STORAGE>
    Matrix& operator+=(const Matrix<T, STORAGE>& m_) {
        *this = *this + m_;
        return *this;
    }
 
    template <typename STORAGE>
    Matrix& operator-=(const Matrix<T, STORAGE>& m_) {
        *this = *this - m_;
        return *this;
    }
 
    Matrix& operator*=(const T& scalar) {
        blas::scal(m.size(), scalar, &m[0], 1);
        return *this;
    }
 
    std::pair<size_t, size_t> size() const { return std::pair<size_t, size_t>(s1, s2); }
 
    size_t size1() const { return s1; }
 
    size_t size2() const { return s2; }
 
    size_t esize() const { return s1 * s2; }
 
    void resize(std::pair<size_t, size_t> s) {
        s1 = s.first;
        s2 = s.second;
        m.resize(s1 * s2);
    }
 
    void resize(size_t s1_, size_t s2_) {
        s1 = s1_;
        s2 = s2_;
        m.resize(s1 * s2);
    }
 
    const T& operator()(size_t i, size_t j) const { return m[i + j * s1]; }
 
    T& operator()(size_t i, size_t j) { return m[i + j * s1]; }
 
    Vector<T, Vdense_ref> operator[](size_t i) { return Vector<T, Vdense_ref>(s1, &m[i * s1]); }
 
    Vector<T, Vdense_constref> operator[](size_t i) const {
        return Vector<T, Vdense_constref>(s1, &m[i * s1]);
    }
 
    Matrix<T, Mrectangle_ref> operator()(const Interval& i) {
        return Matrix<T, Mrectangle_ref>(s1, i.size(), &m[i[0] * s1]);
    }
 
    Matrix<T, Mrectangle_constref> operator()(const Interval& i) const {
        return Matrix<T, Mrectangle_constref>(s1, i.size(), &m[i[0] * s1]);
    }
 
    Matrix<T, Mrectangle_increment_ref> operator()(const Interval& i1, const Interval& i2) {
        return Matrix<T, Mrectangle_increment_ref>(
              i1.size(), i2.size(), &m[i1[0] + i2[0] * s1], s1);
    }
 
    Matrix<T, Mrectangle_increment_constref> operator()(
          const Interval& i1, const Interval& i2) const {
        return Matrix<T, Mrectangle_increment_constref>(
              i1.size(), i2.size(), &m[i1[0] + i2[0] * s1], s1);
    }
 
    Vector<T, Vdense_ref> operator()(const Interval& i1, size_t i2) {
        return Vector<T, Vdense_ref>(i1.size(), &m[i1[0] + i2 * s1]);
    }
 
    Vector<T, Vdense_constref> operator()(const Interval& i1, size_t i2) const {
        return Vector<T, Vdense_constref>(i1.size(), &m[i1[0] + i2 * s1]);
    }
 
    Matrix& operator=(const Matrix<T, Mpacked>& m_) {
        s1 = m_.s;
        s2 = s1;
        m.reserve(s1 * s2);
 
        std::vector<T> cache(s2);
 
        if (m_.upper) {
            auto it1e{std::end(m)};
            auto it2b{std::next(std::cend(m_.m), -s1)};
            auto it2e{std::cend(m_.m)};
 
            for (size_t j = 0; j < s2; ++j) {
                it1e = m.insert(it1e, std::begin(cache), std::next(std::begin(cache), j));
                it1e = m.insert(it1e, it2b, it2e);
 
                it2e = it2b;
                it2b = it2b - (s1 - j - 1);
 
                if (j != s2 - 1) {
                    for (size_t i = 0; i <= j; ++i) { cache[j - i] = m_(s1 - j - 2, s2 - i - 1); }
                }
            }
        } else {
            auto it1e{std::end(m)};
            auto it2b{std::cbegin(m_.m)};
            auto it2e{std::next(std::cbegin(m_.m), s1)};
 
            for (size_t j = 0; j < s2; ++j) {
                m.insert(it1e, std::begin(cache), std::next(std::begin(cache), j));
                it1e = std::end(m);
                m.insert(it1e, it2b, it2e);
 
                it1e = std::end(m);
                it2b = it2e;
                it2e = it2e + (s1 - j - 1);
 
                if (j != s2 - 1) {
                    for (size_t i = 0; i <= j; ++i) { cache[i] = m_(j + 1, i); }
                }
            }
        }
 
        return *this;
    }
 
    Matrix& operator=(const Matrix<
                      T, MMProd<
                               Minverse_constref<Msym_compressed_col_update_inv>,
                               Mrectangle_increment_st_constref> >& m_) {
        if (size() != m_.size()) { resize(m_.size()); }
        if (m_.scale != T(1) || m_.trans) { UnimplementedError() << THROW; }
        m_.m1.m.resolve(
              s1, s2, m_.m2.m, m_.m2.ldm, &m[0], s1, ' ', m_.m1.null_space,
              m_.m1.max_null_space_vector);
        return *this;
    }
 
    Matrix& operator=(const Matrix<
                      T, MMProd<
                               Minverse_constref<Mcompressed_col_update_inv>,
                               Mrectangle_increment_st_constref> >& m_) {
        if (size() != m_.size()) { resize(m_.size()); }
        if (m_.scale != T(1) || m_.trans) { UnimplementedError() << THROW; }
        m_.m1.m.resolve(
              s1, s2, m_.m2.m, m_.m2.ldm, &m[0], s1, ' ', m_.m1.null_space,
              m_.m1.max_null_space_vector);
        return *this;
    }
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        l << "rectangular matrix ";
        l.write(m.s1, m.s2, &m.m[0], m.s1);
        return l;
    }
 
    // TODO: Not optimised.
    // https://en.wikipedia.org/wiki/In-place_matrix_transposition
    void transposed() {
        if (s1 == s2) {
            for (size_t j = 0; j != s2; ++j) {
                for (size_t i = 0; i != j; ++i) { std::swap((*this)(i, j), (*this)(j, i)); }
            }
        } else if (s1 * s2 > 0) {
            std::vector<T> mt;
            mt.reserve(m.size());
            const size_t nm1 = s1 * s2 - 1;
            for (size_t a = 0; a != s1 * s2; ++a) {
                const size_t pa = (a == nm1 ? nm1 : s1 * a % nm1);
                mt.push_back(m[pa]);
            }
            std::swap(m, mt);
            std::swap(s1, s2);
        }
    }
 
    static Matrix null;
 
private:
    size_t s1{};
    size_t s2{};
    std::vector<T> m;
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mrectangle> Matrix<T, Mrectangle>::null;
 
struct Mrectangle_ref {
    using base = Mrectangle_ref;
    using const_base = Mrectangle_increment_st_constref;
    using copy = Mrectangle;
};
 
template <typename T>
class Matrix<T, Mrectangle_ref> {
public:
    Matrix() = default;
 
    Matrix(size_t s1_, size_t s2_, T* m_) : s1(s1_), s2(s2_), m(m_) {}
 
    Matrix(const Matrix& m_) : s1(m_.s1), s2(m_.s2), m(m_.m) {}
 
    Matrix(Matrix<T, Mrectangle>& m_) : s1(m_.s1), s2(m_.s2), m(&m_.m[0]) {}
 
    explicit Matrix(Vector<T, Vdense_ref> v) : s1(v.s), s2(1), m(v.v) {}
 
    Matrix(size_t s1_, size_t s2_, Vector<T, Vdense_ref> v) : s1(s1_), s2(s2_), m(v.v) {
        if (v.s != s1_ * s2_) { Exception() << THROW; }
    }
 
    bool is_null() const { return m == 0; }
 
    void set_zero() { std::fill_n(m, s1 * s2, 0); }
 
    template <typename STORAGE>
    Matrix& operator+=(const Matrix<T, STORAGE>& m_) {
        *this = *this + m_;
        return *this;
    }
 
    template <typename STORAGE>
    Matrix& operator-=(const Matrix<T, STORAGE>& m_) {
        *this = *this - m_;
        return *this;
    }
 
    Matrix& operator=(const Matrix& m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
        std::copy(m_.m, m_.m + s1 * s2, m);
        return *this;
    }
 
    template <typename STORAGE>
    Matrix& operator=(const Matrix<T, STORAGE>& m_) {
        Matrix<T, Mrectangle_increment_ref>(*this) = m_;
        return *this;
    }
 
    std::pair<size_t, size_t> size() const { return std::pair<size_t, size_t>(s1, s2); }
 
    size_t size1() const { return s1; }
 
    size_t size2() const { return s2; }
 
    size_t esize() const { return s1 * s2; }
 
    const T& operator()(size_t i, size_t j) const { return m[i + j * s1]; }
 
    T& operator()(size_t i, size_t j) { return m[i + j * s1]; }
 
    Vector<T, Vdense_ref> operator[](size_t i) { return Vector<T, Vdense_ref>(s1, m + i * s1); }
 
    Vector<T, Vdense_constref> operator[](size_t i) const {
        return Vector<T, Vdense_constref>(s1, m + i * s1);
    }
 
    Matrix<T, Mrectangle_ref> operator()(const Interval& i) {
        return Matrix<T, Mrectangle_ref>(s1, i.size(), m + (i[0] * s1));
    }
 
    Matrix<T, Mrectangle_constref> operator()(const Interval& i) const {
        return Matrix<T, Mrectangle_constref>(s1, i.size(), m + (i[0] * s1));
    }
 
    Matrix<T, Mrectangle_increment_ref> operator()(const Interval& i1, const Interval& i2) {
        return Matrix<T, Mrectangle_increment_ref>(
              i1.size(), i2.size(), m + i1[0] + i2[0] * s1, s1);
    }
 
    Matrix<T, Mrectangle_increment_constref> operator()(
          const Interval& i1, const Interval& i2) const {
        return Matrix<T, Mrectangle_increment_constref>(
              i1.size(), i2.size(), m + i1[0] + i2[0] * s1, s1);
    }
 
    static Matrix null;
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        l << "rectangular matrix ";
        l.write(m.s1, m.s2, m.m, m.s1);
        return l;
    }
 
private:
    size_t s1{};
    size_t s2{};
    T* m{};
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mrectangle_ref> Matrix<T, Mrectangle_ref>::null;
 
struct Mrectangle_constref {
    using base = Mrectangle_constref;
    using const_base = Mrectangle_increment_st_constref;
    using copy = Mrectangle;
};
 
template <typename T>
class Matrix<T, Mrectangle_constref> {
public:
    Matrix() = default;
 
    Matrix(size_t s1_, size_t s2_, const T* m_) : s1(s1_), s2(s2_), m(m_) {}
 
    Matrix(const Matrix& m_) : s1(m_.s1), s2(m_.s2), m(m_.m) {}
 
    Matrix(const Matrix<T, Mrectangle>& m_) : s1(m_.s1), s2(m_.s2), m(&m_.m[0]) {}
 
    Matrix(const Matrix<T, Mrectangle_ref>& m_) : s1(m_.s1), s2(m_.s2), m(m_.m) {}
 
    explicit Matrix(const Vector<T, Vdense_constref>& v) : s1(v.s), s2(1), m(v.v) {}
 
    Matrix(size_t s1_, size_t s2_, const Vector<T, Vdense_constref>& v) : s1(s1_), s2(s2_), m(v.v) {
        if (v.s != s1_ * s2_) { Exception() << THROW; }
    }
 
    bool is_null() const { return m == 0; }
 
    std::pair<size_t, size_t> size() const { return std::pair<size_t, size_t>(s1, s2); }
 
    size_t size1() const { return s1; }
 
    size_t size2() const { return s2; }
 
    size_t esize() const { return s1 * s2; }
 
    const T& operator()(size_t i, size_t j) const { return m[i + j * s1]; }
 
    Vector<T, Vdense_constref> operator[](size_t i) const {
        return Vector<T, Vdense_constref>(s1, m + i * s1);
    }
 
    Matrix<T, Mrectangle_constref> operator()(const Interval& i) const {
        return Matrix<T, Mrectangle_constref>(s1, i.size(), m + (i[0] * s1));
    }
 
    Matrix<T, Mrectangle_increment_constref> operator()(
          const Interval& i1, const Interval& i2) const {
        return Matrix<T, Mrectangle_increment_constref>(
              i1.size(), i2.size(), m + i1[0] + i2[0] * s1, s1);
    }
 
    static Matrix null;
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        l << "rectangular matrix ";
        l.write(m.s1, m.s2, m.m, m.s1);
        return l;
    }
 
private:
    size_t s1{};
    size_t s2{};
    const T* m{};
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mrectangle_constref> Matrix<T, Mrectangle_constref>::null;
 
struct Mrectangle_increment_ref {
    using base = Mrectangle_increment_ref;
    using const_base = Mrectangle_increment_st_constref;
    using copy = Mrectangle;
};
 
template <typename T>
class Matrix<T, Mrectangle_increment_ref> {
public:
    Matrix() = default;
 
    Matrix(size_t s1_, size_t s2_, T* m_, size_t ldm_) : s1(s1_), s2(s2_), m(m_), ldm(ldm_) {}
 
    Matrix(const Matrix& m_) : s1(m_.s1), s2(m_.s2), m(m_.m), ldm(m_.ldm) {}
 
    Matrix(Matrix<T, Mrectangle>& m_) : s1(m_.s1), s2(m_.s2), m(&m_.m[0]), ldm(m_.s1) {}
 
    Matrix(Matrix<T, Mrectangle_ref> m_) : s1(m_.s1), s2(m_.s2), m(m_.m), ldm(m_.s1) {}
 
    Matrix(size_t s1_, size_t s2_, Vector<T, Vincrement_ref> v)
        : s1(s1_), s2(s2_), m(v.v), ldm(s1_) {
        if (v.s != s1_ * s2_) { Exception() << THROW; }
    }
 
    bool is_null() const { return m == 0; }
 
    void set_zero() {
        T* ptr = m;
        for (size_t j = 0; j != s2; ++j, ptr += ldm) { std::fill_n(ptr, s1, 0.0); }
    }
 
    std::pair<size_t, size_t> size() const { return std::pair<size_t, size_t>(s1, s2); }
 
    size_t size1() const { return s1; }
 
    size_t size2() const { return s2; }
 
    size_t esize() const { return s1 * s2; }
 
    const T& operator()(size_t i, size_t j) const { return m[i + j * ldm]; }
 
    T& operator()(size_t i, size_t j) { return m[i + j * ldm]; }
 
    Vector<T, Vdense_ref> operator[](size_t i) { return Vector<T, Vdense_ref>(s1, m + i * ldm); }
 
    Vector<T, Vdense_constref> operator[](size_t i) const {
        return Vector<T, Vdense_constref>(s1, m + i * ldm);
    }
 
    Matrix operator()(const Interval& i1, const Interval& i2) {
        return Matrix(i1.size(), i2.size(), m + i1[0] + i2[0] * ldm, s1);
    }
 
    Matrix<T, Mrectangle_increment_constref> operator()(
          const Interval& i1, const Interval& i2) const {
        return Matrix<T, Mrectangle_increment_constref>(
              i1.size(), i2.size(), m + i1[0] + i2[0] * ldm, s1);
    }
 
    Matrix& operator*=(T s) {
        T* ptr = m;
        for (size_t j = 0; j != s2; ++j) {
            T* ptr_end = ptr + s1;
            for (; ptr != ptr_end; ++ptr) { *ptr *= s; }
            ptr += ldm - s1;
        }
        return *this;
    }
 
    template <typename STORAGE>
    Matrix& operator+=(const Matrix<T, STORAGE>& m_) {
        *this = *this + m_;
        return *this;
    }
 
    template <typename STORAGE>
    Matrix& operator-=(const Matrix<T, STORAGE>& m_) {
        *this = *this - m_;
        return *this;
    }
 
    Matrix& operator=(const Matrix& m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
        if (m == m_.m && ldm == m_.ldm) { return *this; }
        if (ldm == s1 && m_.ldm == m_.s1) {
            T* m1 = m;
            const T* m2 = m_.m;
            const T* m2_end = m2 + s1 * s2;
            while (m2 != m2_end) { *(m1++) = *(m2++); }
        } else {
            for (size_t i = 0; i != s2; ++i) { (*this)[i] = m_[i]; }
        }
        return *this;
    }
 
    Matrix& operator=(const Matrix<T, Mrectangle_constref>& m_) {
        if (size1() != m_.size1() || size2() != m_.size2()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
        const T* ptrin = m_.m;
        T* ptrout = m;
        for (size_t j = 0; j != s2; ++j, ptrout += ldm - s1) {
            for (size_t i = 0; i != s1; ++i, ++ptrout, ++ptrin) { *ptrout = *ptrin; }
        }
        return *this;
    }
 
    Matrix& operator=(const Matrix<T, Mpacked_st_constref>& m_) {
        if (size1() != m_.size1() || size2() != m_.size2()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
        for (size_t j = 0; j != size2(); ++j) {
            for (size_t i = 0; i != size1(); ++i) { (*this)(i, j) = m_(i, j); }
        }
        return *this;
    }
 
    Matrix& operator=(const Matrix<T, Mrectangle_increment_st_constref>& m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
        if (m == m_.m && ldm == m_.ldm && m_.trans == false && m_.scale == T(1)) { return *this; }
        if (m_.trans == false && ldm == s1 && m_.ldm == m_.s1) {
            T* m1 = m;
            const T* m2 = m_.m;
            const T* m2_end = m2 + s1 * s2;
            while (m2 != m2_end) { *(m1++) = m_.scale * *(m2++); }
        } else {
            for (size_t i = 0; i != s2; ++i) { (*this)[i] = m_[i]; }
        }
        return *this;
    }
 
    Matrix& operator=(
          const Matrix<T, Sum<Mrectangle_increment_st_constref, Mrectangle_increment_st_constref> >&
                m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
        if (!m_.m1.trans && !m_.m2.trans && ldm == s1 && m_.m1.ldm == s1 && m_.m2.ldm == s1
            && m_.m1.scale == T(1) && m_.m2.scale == T(1)) {
            const size_t i_end = s1 * s2;
            for (size_t i = 0; i != i_end; ++i) { m[i] = m_.m1.m[i] + m_.m2.m[i]; }
        } else {
            for (size_t i = 0; i != s2; ++i) { (*this)[i] = m_.scale * (m_.m1[i] + m_.m2[i]); }
        }
        return *this;
    }
 
    Matrix& operator=(
          const Matrix<
                T, MMProd<Mrectangle_increment_st_constref, Mrectangle_increment_st_constref> >&
                m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
        if (m_.trans) {
            blas::gemm(
                  m_.m2.trans ? 'n' : 't', m_.m1.trans ? 'n' : 't', s1, s2, m_.m2.size1(),
                  m_.scale * m_.m1.scale * m_.m2.scale, m_.m2.m, m_.m2.ldm, m_.m1.m, m_.m1.ldm, 0,
                  m, ldm);
        } else {
            blas::gemm(
                  m_.m1.trans ? 't' : 'n', m_.m2.trans ? 't' : 'n', s1, s2, m_.m1.size2(),
                  m_.scale * m_.m1.scale * m_.m2.scale, m_.m1.m, m_.m1.ldm, m_.m2.m, m_.m2.ldm, 0,
                  m, ldm);
        }
        return *this;
    }
 
    Matrix& operator=(
          const Matrix<
                T,
                Sum<Mrectangle_increment_st_constref,
                    MMProd<Mrectangle_increment_st_constref, Mrectangle_increment_st_constref> > >&
                m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
        if (m_.trans) {
            (*this) = transposed(m_.scale * m_.m1);
        } else {
            (*this) = m_.scale * m_.m1;
        }
        if (m_.m2.trans ^ m_.trans) {
            blas::gemm(
                  m_.m2.m2.trans ? 'n' : 't', m_.m2.m1.trans ? 'n' : 't', s1, s2, m_.m2.m2.size1(),
                  m_.scale * m_.m2.scale * m_.m2.m1.scale * m_.m2.m2.scale, m_.m2.m2.m,
                  m_.m2.m2.ldm, m_.m2.m1.m, m_.m2.m1.ldm, m_.scale * m_.m1.scale, m, ldm);
        } else {
            blas::gemm(
                  m_.m2.m1.trans ? 't' : 'n', m_.m2.m2.trans ? 't' : 'n', s1, s2, m_.m2.m1.size2(),
                  m_.scale * m_.m2.scale * m_.m2.m1.scale * m_.m2.m2.scale, m_.m2.m1.m,
                  m_.m2.m1.ldm, m_.m2.m2.m, m_.m2.m2.ldm, m_.scale * m_.m1.scale, m, ldm);
        }
        return *this;
    }
 
    Matrix& operator=(
          const Matrix<T, MMProd<Mcompressed_col_st_constref, Mrectangle_increment_st_constref> >&
                m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
        for (size_t j = 0; j != s2; ++j) { (*this)[j] = m_.m1 * m_.m2[j]; }
        return *this;
    }
 
    Matrix& operator=(
          const Matrix<
                T, Sum<Mrectangle_increment_st_constref,
                       MMProd<Mcompressed_col_st_constref, Mrectangle_increment_st_constref> > >&
                m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
 
        for (size_t i = 0; i != s2; ++i) {
            (*this)[i] = m_.scale * m_.m1[i] + (m_.scale * m_.m2.scale) * (m_.m2.m1 * m_.m2.m2[i]);
        }
        return *this;
    }
 
    template <typename ATYPE, typename BTYPE>
    Matrix& operator=(const Matrix<T, Eigenvector<ATYPE, BTYPE> >& eigenvector) {
        if (size() != eigenvector.size()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
        const_cast<Matrix<T, Eigenvector<ATYPE, BTYPE> >&>(eigenvector).resolve(*this);
        return *this;
    }
 
    template <typename ATYPE, typename BTYPE>
    Matrix& operator=(const Matrix<T, PolyEigenvector<ATYPE, BTYPE> >& eigenvector) {
        if (size() != eigenvector.size()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
        const_cast<Matrix<T, PolyEigenvector<ATYPE, BTYPE> >&>(eigenvector).resolve(*this);
        return *this;
    }
 
    Matrix& operator=(
          const Matrix<
                T, MMProd<Msym_compressed_col_st_constref, Mrectangle_increment_st_constref> >&
                m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
 
        for (size_t i = 0; i != s2; ++i) { (*this)[i] = m_.scale * (m_.m1 * m_.m2[i]); }
        return *this;
    }
 
    Matrix& operator=(
          const Matrix<
                T,
                Sum<Mrectangle_increment_st_constref,
                    MMProd<Msym_compressed_col_st_constref, Mrectangle_increment_st_constref> > >&
                m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
 
        for (size_t i = 0; i != s2; ++i) {
            (*this)[i] = m_.scale * m_.m1[i] + (m_.scale * m_.m2.scale) * (m_.m2.m1 * m_.m2.m2[i]);
        }
        return *this;
    }
 
    template <typename M_FORMAT>
    Matrix& operator+=(
          const Matrix<
                T, MMProd<Mstructured_constref<M_FORMAT>, Mrectangle_increment_st_constref> >& m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
 
        if (m_.scale == T(1) && m_.trans == false && m_.m2.trans == false) {
            for (size_t k = 0; k != s2; ++k) {
                Vector<T, Vindex_scale_constref> tmp1(
                      m_.m1.m.size1(), m_.m2.m + k * m_.m2.ldm, m_.m1.index, m_.m2.scale);
                Vector<T, Vindex_ref> tmp2(m_.m1.m.size1(), &(*this)(0, k), m_.m1.index);
                tmp2 += m_.m1.m * tmp1;
            }
        } else {
            UnimplementedError() << THROW;
        }
        return *this;
    }
 
    template <typename M_FORMAT>
    Matrix& operator=(
          const Matrix<
                T, Sum<Mrectangle_increment_st_constref,
                       MMProd<Mstructured_constref<M_FORMAT>, Mrectangle_increment_st_constref> > >&
                m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
 
        if (m_.m1.m == m && m_.m1.ldm == ldm && m_.m1.scale == T(1) && m_.m1.trans == false
            && m_.scale == T(1) && m_.trans == false && m_.m2.m2.trans == false) {
            for (size_t k = 0; k != s2; ++k) {
                Vector<T, Vindex_scale_constref> tmp1(
                      m_.m2.m1.m.size1(), m_.m2.m2.m + k * m_.m2.m2.ldm, m_.m2.m1.index,
                      m_.m2.m2.scale);
                Vector<T, Vindex_ref> tmp2(m_.m2.m1.m.size1(), &(*this)(0, k), m_.m2.m1.index);
                tmp2 = Vector<T, Vindex_scale_constref>(tmp2) + m_.m2.m1.m * tmp1;
            }
        } else {
            UnimplementedError() << THROW;
        }
        return *this;
    }
 
    template <typename M_FORMAT>
    Matrix& operator=(
          const Matrix<T, Sum<Mrectangle_increment_st_constref, Mstructured_constref<M_FORMAT> > >&
                m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
 
        if (m_.m1.m != m) {
            (*this) = m_.scale * m_.m1;
        } else if (m_.m1.ldm != ldm) {
            UnimplementedError() << THROW;
        } else if (m_.scale != T(1)) {
            (*this) *= m_.scale;
        }
        for (size_t j = 0; j != m_.m2.m.size2(); ++j) {
            const size_t jj = (m_.m2.index2 ? m_.m2.index2 : m_.m2.index)[j];
            for (size_t i = 0; i != m_.m2.m.size1(); ++i) {
                (*this)(m_.m2.index[i], jj) += m_.m2.m(i, j) * m_.scale;
            }
        }
        return *this;
    }
 
    Matrix& operator=(
          const Matrix<T, MMProd<Mrectangle_inv, Mrectangle_increment_st_constref> >& m_) {
        if (m_.scale != T(1)) { UnimplementedError() << THROW; }
        *this = m_.m2;
        m_.m1.resolve(*this);
        return *this;
    }
 
    Matrix& operator=(
          const Matrix<T, MMProd<Mpacked_st_constref, Mrectangle_increment_st_constref> >& m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
        const T alpha = m_.m1.scale * m_.m2.scale * m_.scale;
        const char m1_up = m_.m1.upper ? 'U' : 'L';
        if (m_.m2.trans) {
            for (size_t i = 0; i != s1; ++i) {
                blas::spmv(
                      m1_up, m_.m1.s, alpha, m_.m1.m, &m_.m2.m[i], m_.m2.ldm, 0, &m[i * ldm], 1);
            }
        } else {
            for (size_t j = 0; j != s2; ++j) {
                blas::spmv(
                      m1_up, m_.m1.s, alpha, m_.m1.m, &m_.m2.m[j * m_.m2.ldm], 1, 0, &m[j * ldm],
                      1);
            }
        }
        return *this;
    }
 
    Matrix& operator=(
          const Matrix<T, MMProd<Mcompressed_col_st_constref, Mcompressed_col_st_constref> >& m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
        if (m_.m1.trans || m_.m2.trans) { UnimplementedError() << THROW; }
        const T scale = m_.scale * m_.m1.scale * m_.m2.scale;
        T* ptr = m;
        for (size_t j = 0; j != s2; ++j, ptr += ldm) {
            std::fill_n(ptr, s1, T(0));
            size_t k = m_.m2.si[j];
            const size_t k_end = m_.m2.si[j + 1];
            for (; k != k_end; ++k) {
                const size_t kk = m_.m2.index[k];
                const T kv = scale * m_.m2.m[k];
                size_t i = m_.m1.si[kk];
                const size_t i_end = m_.m1.si[kk + 1];
                for (; i != i_end; ++i) { ptr[m_.m1.index[i]] += kv * m_.m1.m[i]; }
            }
        }
        return *this;
    }
 
    void in_place_invert() {
        if (s1 != s2) { Exception() << THROW; }
        int info;
        int* ipiv = TemporaryBuffer<int>::get(s1);
        lapack::getrf(s1, s1, m, ldm, ipiv, info);
        if (info != 0) { Exception() << THROW; }
        lapack::getri(s1, m, ldm, ipiv, TemporaryBuffer<double>::get(s1 * s1), s1 * s1, info);
        if (info != 0) { Exception() << THROW; }
    }
 
    size_t diagonalise_base_vector_on_observation_matrix(
          const Index& observation_matrix, const double tol = 1e-13) {
        if (observation_matrix.size() > s2) { Exception() << THROW; }
        Matrix<T, Mrectangle> u(observation_matrix.size(), s2);
        {
            T* orig = m;
            T* dest = &u(0, 0);
            for (size_t j = 0; j != s2; ++j, orig += ldm) {
                for (size_t i = 0; i != observation_matrix.size(); ++i, ++dest) {
                    *dest = orig[observation_matrix[i]];
                }
            }
        }
        Vector<T, Vdense> sv(u.size1());
        Matrix<T, Mrectangle> vt(s2, s2);
        int info;
        int lwork;
        {
            T lwork_tmp;
            lapack::gesvd(
                  'O', 'A', u.size1(), s2, &u(0, 0), u.size1(), &sv[0], NULL, 1, &vt(0, 0), s2,
                  &lwork_tmp, -1, info);
            lwork = int(lwork_tmp);
        }
        if (info != 0) { Exception() << THROW; }
        {
            std::vector<T> work(lwork);
            lapack::gesvd(
                  'O', 'A', u.size1(), s2, &u(0, 0), u.size1(), &sv[0], NULL, 1, &vt(0, 0), s2,
                  &work[0], lwork, info);
            if (info < 0) { Exception() << THROW; }
            if (info > 0) {
                Exception() << "could not compute the svd decomposition of the matrix." << THROW;
            }
            u.resize(u.size1(), u.size1());
        }
 
        Matrix<T, Mrectangle> tmp;
        tmp = (*this) * transposed(vt);
        size_t rank;
        for (rank = 0; rank != sv.size(); ++rank) {
            if (std::abs(sv[rank]) < tol) { break; }
            tmp[rank] *= 1.0 / sv[rank];
        }
        Interval inter(0, u.size1());
        (*this)(Interval(0, s1), inter) = tmp(Interval(0, s1), inter) * transposed(u(inter, inter));
        (*this)(Interval(0, s1), Interval(u.size1(), s2)) =
              tmp(Interval(0, s1), Interval(u.size1(), s2));
 
        return s2 + rank - observation_matrix.size();
    }
 
    static Matrix null;
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        l << "rectangular matrix ";
        l.write(m.s1, m.s2, m.m, m.ldm);
        return l;
    }
 
private:
    size_t s1{};
    size_t s2{};
    T* m{};
    size_t ldm{};
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mrectangle_increment_ref> Matrix<T, Mrectangle_increment_ref>::null;
 
struct Mrectangle_increment_constref {
    using base = Mrectangle_increment_constref;
    using const_base = Mrectangle_increment_st_constref;
    using copy = Mrectangle;
};
 
template <typename T>
class Matrix<T, Mrectangle_increment_constref> {
public:
    Matrix() = default;
 
    Matrix(size_t s1_, size_t s2_, const T* m_, size_t ldm_)
        : s1(s1_), s2(s2_), m(m_), ldm(std::max(size_t(1), ldm_)) {}
 
    Matrix(const Matrix& m_) : s1(m_.s1), s2(m_.s2), m(m_.m), ldm(m_.ldm) {}
 
    Matrix(const Matrix<T, Mrectangle>& m_)
        : s1(m_.s1), s2(m_.s2), m(&m_.m[0]), ldm(std::max(size_t(1), m_.s1)) {}
 
    Matrix(const Matrix<T, Mrectangle_ref>& m_)
        : s1(m_.s1), s2(m_.s2), m(m_.m), ldm(std::max(size_t(1), m_.s1)) {}
 
    Matrix(const Matrix<T, Mrectangle_constref>& m_)
        : s1(m_.s1), s2(m_.s2), m(m_.m), ldm(std::max(size_t(1), m_.s1)) {}
 
    Matrix(const Matrix<T, Mrectangle_increment_ref>& m_)
        : s1(m_.s1), s2(m_.s2), m(m_.m), ldm(m_.ldm) {}
 
    Matrix(size_t s1_, size_t s2_, const Vector<T, Vincrement_constref>& v)
        : s1(s1_), s2(s2_), m(v.v), ldm(std::max(size_t(1), s1_)) {
        if (v.s != s1 * s2) { Exception() << THROW; }
    }
 
    bool is_null() const { return m == 0; }
 
    std::pair<size_t, size_t> size() const { return std::pair<size_t, size_t>(s1, s2); }
 
    size_t size1() const { return s1; }
 
    size_t size2() const { return s2; }
 
    size_t esize() const { return s1 * s2; }
 
    const T& operator()(size_t i, size_t j) const { return m[i + j * ldm]; }
 
    Vector<T, Vdense_constref> operator[](size_t i) const {
        return Vector<T, Vdense_constref>(s1, m + i * ldm);
    }
 
    Matrix<T, Mrectangle_increment_constref> operator()(
          const Interval& i1, const Interval& i2) const {
        return Matrix<T, Mrectangle_increment_constref>(
              i1.size(), i2.size(), m + i1[0] + i2[0] * ldm, ldm);
    }
 
    static Matrix null;
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        l << "rectangular matrix ";
        l.write(m.s1, m.s2, m.m, m.ldm);
        return l;
    }
 
private:
    size_t s1{};
    size_t s2{};
    const T* m{};
    size_t ldm{};
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mrectangle_increment_constref> Matrix<T, Mrectangle_increment_constref>::null;
 
struct Mrectangle_increment_st_constref {
    using base = Mrectangle_increment_st_constref;
    using const_base = Mrectangle_increment_st_constref;
    using copy = Mrectangle;
};
 
template <typename T>
class Matrix<T, Mrectangle_increment_st_constref> {
public:
    Matrix() = default;
 
    Matrix(size_t s1_, size_t s2_, const T* m_, size_t ldm_, T scale_, bool transpose_ = false)
        : s1(s1_),
          s2(s2_),
          m(m_),
          ldm(std::max(size_t(1), ldm_)),
          scale(scale_),
          trans(transpose_) {}
 
    Matrix(const Matrix& m_)
        : s1(m_.s1), s2(m_.s2), m(m_.m), ldm(m_.ldm), scale(m_.scale), trans(m_.trans) {}
 
    Matrix(const Matrix<T, Mrectangle>& m_)
        : s1(m_.s1),
          s2(m_.s2),
          m(&m_.m[0]),
          ldm(std::max(size_t(1), m_.s1)),
          scale(1),
          trans(false) {}
 
    Matrix(const Matrix<T, Mrectangle_ref>& m_)
        : s1(m_.s1), s2(m_.s2), m(m_.m), ldm(std::max(size_t(1), m_.s1)), scale(1), trans(false) {}
 
    Matrix(const Matrix<T, Mrectangle_constref>& m_)
        : s1(m_.s1), s2(m_.s2), m(m_.m), ldm(std::max(size_t(1), m_.s1)), scale(1), trans(false) {}
 
    Matrix(const Matrix<T, Mrectangle_increment_ref>& m_)
        : s1(m_.s1), s2(m_.s2), m(m_.m), ldm(m_.ldm), scale(1), trans(false) {}
 
    Matrix(const Matrix<T, Mrectangle_increment_constref>& m_)
        : s1(m_.s1), s2(m_.s2), m(m_.m), ldm(m_.ldm), scale(1), trans(false) {}
 
    Matrix(size_t s1_, size_t s2_, const Vector<T, Vincrement_scale_constref>& v)
        : s1(s1_), s2(s2_), m(v.v), ldm(std::max(size_t(1), s1_)), scale(v.scale), trans(false) {
        if (v.s != s1 * s2) { Exception() << THROW; }
    }
 
    bool is_null() const { return m == 0; }
 
    std::pair<size_t, size_t> size() const {
        if (trans) {
            return std::pair<size_t, size_t>(s2, s1);
        } else {
            return std::pair<size_t, size_t>(s1, s2);
        }
    }
 
    size_t size1() const {
        if (trans) {
            return s2;
        } else {
            return s1;
        }
    }
 
    size_t size2() const {
        if (trans) {
            return s1;
        } else {
            return s2;
        }
    }
 
    size_t esize() const { return s1 * s2; }
 
    T operator()(size_t i, size_t j) const {
        if (trans) {
            return scale * m[j + i * ldm];
        } else {
            return scale * m[i + j * ldm];
        }
    }
 
    Vector<T, Vincrement_scale_constref> operator[](size_t i) const {
        if (trans) {
            return Vector<T, Vincrement_scale_constref>(s2, m + i, ldm, scale);
        } else {
            return Vector<T, Vincrement_scale_constref>(s1, m + i * ldm, 1, scale);
        }
    }
 
    Matrix<T, Mrectangle_increment_st_constref> operator()(
          const Interval& i1, const Interval& i2) const {
        if (trans) {
            return Matrix<T, Mrectangle_increment_st_constref>(
                  i2.size(), i1.size(), m + i2[0] + i1[0] * ldm, ldm, scale, trans);
        } else {
            return Matrix<T, Mrectangle_increment_st_constref>(
                  i1.size(), i2.size(), m + i1[0] + i2[0] * ldm, ldm, scale, trans);
        }
    }
 
    Matrix& operator*(T s) {
        scale *= s;
        return *this;
    }
 
    Matrix& transpose() {
        trans = trans ? false : true;
        return *this;
    }
 
    static Matrix null;
 
private:
    size_t s1{};
    size_t s2{};
    const T* m{};
    size_t ldm{};
    T scale{};
    bool trans{};
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mrectangle_increment_st_constref> Matrix<T, Mrectangle_increment_st_constref>::null;
 
struct Mrectangle_inv {
    using base = Mrectangle_inv;
    using const_base = Mrectangle_inv;
    using copy = Mrectangle_inv;
};
 
template <typename T>
class Matrix<T, Mrectangle_inv> {
public:
    Matrix(const Matrix<T, Mrectangle>& m_) : s(m_.s1), m(m_.m), ipiv(m_.s1) {
        if (m_.s1 != m_.s2) { UnimplementedError() << THROW; }
        int info;
        lapack::getrf(s, s, &m[0], s, &ipiv[0], info);
        if (info < 0) { Exception() << THROW; }
        if (info > 0) { Exception() << "The matrix is singular." << THROW; }
    }
 
    size_t size1() const { return s; }
 
    size_t size2() const { return s; }
 
    std::pair<size_t, size_t> size() const { return std::pair<size_t, size_t>(s, s); }
 
    void resolve(Vector<T, Vincrement_ref> v_) const {
        if (v_.incr != 1) { UnimplementedError() << THROW; }
        int info;
        lapack::getrs('N', s, 1, &m[0], s, &ipiv[0], v_.v, s, info);
        if (info < 0) { Exception() << THROW; }
    }
 
    void resolve(Matrix<T, Mrectangle_increment_ref> m_) const {
        int info;
        lapack::getrs('N', s, m_.size2(), &m[0], s, &ipiv[0], m_.m, m_.ldm, info);
        if (info < 0) { Exception() << THROW; }
    }
 
private:
    size_t s;
    std::vector<T> m;
    std::vector<int> ipiv;
    MVFRIEND;
};
 
}  // namespace b2000::b2linalg
 
#endif