b2matrix_packed.H

//------------------------------------------------------------------------
// b2matrix_packed.H --
//
//     A framework for generic linear algebraic (matrix) computations.
//
// written by Mathias Doreille
//            Thomas Blome <thomas.blome@dlr.de>
//
// Copyright (c) 2004-2012 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_PACKED_H__
#define __B2MATRIX_PACKED_H__
 
#include "b2linear_algebra_def.H"
 
namespace b2000::b2linalg {
 
struct Mlower_packed {
    using base = Mpacked_ref;
    using const_base = Mpacked_st_constref;
    using copy = Mlower_packed;
};
 
template <typename T>
class Matrix<T, Mlower_packed> {
public:
    Matrix() = default;
 
    Matrix(size_t s_) : s(s_), m(s_ * (s_ + 1) / 2) {}
 
    Matrix(const Matrix& m_) : s(m_.s), m(m_.m) {}
 
    template <typename T1>
    Matrix(const Matrix<T1, Mlower_packed>& m_)
        : s(m_.s.begin(), m_.s.end()), m(m_.m.begin(), m_.m.end()) {}
 
    template <typename T1>
    Matrix& operator=(const Matrix<T1, Mlower_packed>& m_) {
        s = m_.s;
        m.assign(m_.m.begin(), m_.m.end());
        return *this;
    }
 
    template <typename STORAGE>
    Matrix& operator=(const Matrix<T, STORAGE>& m_) {
        if (m_.size1() != m_.size2()) { Exception() << "The matrix is not symmetric." << THROW; }
        s = m_.size1();
        m.resize(s * (s + 1) / 2);
        Matrix<T, Mpacked_ref> ref(*this);
        ref = m_;
        return *this;
    }
 
    Matrix& operator=(const Matrix& m_) {
        s = m_.s;
        m = m_.m;
        return *this;
    }
 
    bool is_null() const { return this == &null; }
 
    void set_zero() { std::fill(m.begin(), m.end(), 0); }
 
    std::pair<size_t, size_t> size() const { return std::pair<size_t, size_t>(s, s); }
 
    size_t size1() const { return s; }
 
    size_t size2() const { return s; }
 
    size_t esize() const { return s * (s + 1) / 2; }
 
    void resize(size_t s_) {
        s = s_;
        m.resize(s_ * (s_ + 1) / 2);
    }
 
    void resize(size_t s1, size_t s2) {
        if (s1 != s2) { Exception() << THROW; }
        s = s1;
        m.resize(s1 * (s1 + 1) / 2);
    }
 
    const T& operator()(size_t i, size_t j) const {
        if (i < j) {
            return m[(2 * s - i - 1) * i / 2 + j];
        } else {
            return m[(2 * s - j - 1) * j / 2 + i];
        }
    }
 
    T& operator()(size_t i, size_t j) {
        if (i < j) {
            return m[(2 * s - i - 1) * i / 2 + j];
        } else {
            return m[(2 * s - j - 1) * j / 2 + i];
        }
    }
 
    static Matrix null;
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        l << "lower packed matrix ";
        l.write(m.s * (m.s + 1) / 2, &m.m[0], 1);
        return l;
    }
 
private:
    size_t s{};
    std::vector<T> m{};
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mlower_packed> Matrix<T, Mlower_packed>::null;
 
struct Mlower_packed_ref {
    using base = Mpacked_ref;
    using const_base = Mpacked_st_constref;
    using copy = Mlower_packed;
};
 
template <typename T>
class Matrix<T, Mlower_packed_ref> {
public:
    Matrix() = default;
 
    Matrix(size_t s_, T* m_) : s(s_), m(m_) {}
 
    Matrix(const Matrix& m_) : s(m_.s), m(m_.m) {}
 
    Matrix(Matrix<T, Mlower_packed>& m_) : s(m_.s), m(&m_.m[0]) {}
 
    bool is_null() const { return m == 0; }
 
    void set_zero() { std::fill(m, m + esize(), 0); }
 
    Matrix& operator=(const Matrix& m_) {
        if (s != m_.s) { Exception() << "The two matrix do not have the same size." << THROW; }
        std::copy(m_.m, m_.m + s, m);
        return *this;
    }
 
    template <typename STORAGE>
    Matrix& operator=(const Matrix<T, STORAGE>& m_) {
        Matrix<T, Mpacked_ref> ref(*this);
        ref = m_;
        return *this;
    }
 
    std::pair<size_t, size_t> size() const { return std::pair<size_t, size_t>(s, s); }
 
    size_t size1() const { return s; }
 
    size_t size2() const { return s; }
 
    size_t esize() const { return s * (s + 1) / 2; }
 
    const T& operator()(size_t i, size_t j) const {
        if (i < j) {
            return m[(2 * s - i - 1) * i / 2 + j];
        } else {
            return m[(2 * s - j - 1) * j / 2 + i];
        }
    }
 
    T& operator()(size_t i, size_t j) {
        if (i < j) {
            return m[(2 * s - i - 1) * i / 2 + j];
        } else {
            return m[(2 * s - j - 1) * j / 2 + i];
        }
    }
 
    static Matrix null;
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        l << "lower packed matrix ";
        l.write(m.s * (m.s + 1) / 2, m.m, 1);
        return l;
    }
 
private:
    size_t s{};
    T* m{};
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mlower_packed_ref> Matrix<T, Mlower_packed_ref>::null;
 
struct Mlower_packed_constref {
    using base = Mpacked_st_constref;
    using const_base = Mpacked_st_constref;
    using copy = Mlower_packed;
};
 
template <typename T>
class Matrix<T, Mlower_packed_constref> {
public:
    Matrix() = default;
 
    Matrix(size_t s_, const T* m_) : s(s_), m(m_) {}
 
    Matrix(const Matrix& m_) : s(m_.s), m(m_.m) {}
 
    Matrix(const Matrix<T, Mlower_packed>& m_) : s(m_.s), m(&m_.m[0]) {}
 
    Matrix(const Matrix<T, Mlower_packed_ref>& m_) : s(m_.s), m(m_.m) {}
 
    bool is_null() const { return m == 0; }
 
    std::pair<size_t, size_t> size() const { return std::pair<size_t, size_t>(s, s); }
 
    size_t size1() const { return s; }
 
    size_t size2() const { return s; }
 
    size_t esize() const { return s * (s + 1) / 2; }
 
    const T& operator()(size_t i, size_t j) const {
        if (i < j) {
            return m[(2 * s - i - 1) * i / 2 + j];
        } else {
            return m[(2 * s - j - 1) * j / 2 + i];
        }
    }
 
    static Matrix null;
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        l << "lower packed matrix ";
        l.write(m.s * (m.s + 1) / 2, m.m, 1);
        return l;
    }
 
private:
    size_t s{};
    const T* m{};
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mlower_packed_constref> Matrix<T, Mlower_packed_constref>::null;
 
struct Mupper_packed {
    using base = Mpacked_ref;
    using const_base = Mpacked_st_constref;
    using copy = Mupper_packed;
};
 
template <typename T>
class Matrix<T, Mupper_packed> {
public:
    Matrix() = default;
 
    Matrix(size_t s_) : s(s_), m(s_ * (s_ + 1) / 2) {}
 
    Matrix(const Matrix& m_) : s(m_.s), m(m_.m) {}
 
    template <typename T1>
    Matrix(const Matrix<T1, Mupper_packed>& m_)
        : s(m_.s.begin(), m_.s.end()), m(m_.m.begin(), m_.m.end()) {}
 
    template <typename T1>
    Matrix& operator=(const Matrix<T1, Mupper_packed>& m_) {
        s = m_.s;
        m.assign(m_.m.begin(), m_.m.end());
        return *this;
    }
 
    Matrix& operator=(const Matrix& m_) {
        s = m_.s;
        m = m_.m;
        return *this;
    }
 
    template <typename STORAGE>
    Matrix& operator=(const Matrix<T, STORAGE>& m_) {
        if (m_.size1() != m_.size2()) { Exception() << "The matrix is not symmetric." << THROW; }
        s = m_.size1();
        m.resize(s * (s + 1) / 2);
        Matrix<T, Mpacked_ref> ref(*this);
        ref = m_;
        return *this;
    }
 
    bool is_null() const { return this == &null; }
 
    void set_zero() { std::fill(m.begin(), m.end(), 0); }
 
    std::pair<size_t, size_t> size() const { return std::pair<size_t, size_t>(s, s); }
 
    size_t size1() const { return s; }
 
    size_t size2() const { return s; }
 
    size_t esize() const { return s * (s + 1) / 2; }
 
    void resize(size_t s1, size_t s2) {
        if (s1 != s2) { Exception() << THROW; }
        s = s1;
        m.resize(s1 * (s1 + 1) / 2);
    }
 
    const T& operator()(size_t i, size_t j) const {
        if (j < i) {
            return m[i * (i + 1) / 2 + j];
        } else {
            return m[j * (j + 1) / 2 + i];
        }
    }
 
    T& operator()(size_t i, size_t j) {
        if (j < i) {
            return m[i * (i + 1) / 2 + j];
        } else {
            return m[j * (j + 1) / 2 + i];
        }
    }
 
    static Matrix null;
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        l << "upper packed matrix ";
        l.write(m.s * (m.s + 1) / 2, &m.m[0], 1);
        return l;
    }
 
private:
    size_t s{};
    std::vector<T> m{};
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mupper_packed> Matrix<T, Mupper_packed>::null;
 
struct Mupper_packed_ref {
    using base = Mpacked_ref;
    using const_base = Mpacked_st_constref;
    using copy = Mupper_packed;
};
 
template <typename T>
class Matrix<T, Mupper_packed_ref> {
public:
    Matrix() = default;
 
    Matrix(size_t s_, T* m_) : s(s_), m(m_) {}
 
    Matrix(const Matrix& m_) : s(m_.s), m(m_.m) {}
 
    Matrix(Matrix<T, Mupper_packed>& m_) : s(m_.s), m(&m_.m[0]) {}
 
    bool is_null() const { return m == 0; }
 
    void set_zero() { std::fill(m, m + esize(), 0); }
 
    Matrix& operator=(const Matrix& m_) {
        if (s != m_.s) { Exception() << "The two matrix do not have the same size." << THROW; }
        std::copy(m_.m, m_.m + s, m);
        return *this;
    }
 
    template <typename STORAGE>
    Matrix& operator=(const Matrix<T, STORAGE>& m_) {
        Matrix<T, Mpacked_ref> ref(*this);
        ref = m_;
        return *this;
    }
 
    std::pair<size_t, size_t> size() const { return std::pair<size_t, size_t>(s, s); }
 
    size_t size1() const { return s; }
 
    size_t size2() const { return s; }
 
    size_t esize() const { return s * (s + 1) / 2; }
 
    const T& operator()(size_t i, size_t j) const {
        if (j < i) {
            return m[i * (i + 1) / 2 + j];
        } else {
            return m[j * (j + 1) / 2 + i];
        }
    }
 
    T& operator()(size_t i, size_t j) {
        if (j < i) {
            return m[i * (i + 1) / 2 + j];
        } else {
            return m[j * (j + 1) / 2 + i];
        }
    }
 
    static Matrix null;
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        l << "upper packed matrix ";
        l.write(m.s * (m.s + 1) / 2, m.m, 1);
        return l;
    }
 
private:
    size_t s{};
    T* m{};
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mupper_packed_ref> Matrix<T, Mupper_packed_ref>::null;
 
struct Mupper_packed_constref {
    using base = Mpacked_st_constref;
    using const_base = Mpacked_st_constref;
    using copy = Mupper_packed;
};
 
template <typename T>
class Matrix<T, Mupper_packed_constref> {
public:
    Matrix() = default;
 
    Matrix(size_t s_, const T* m_) : s(s_), m(m_) {}
 
    Matrix(const Matrix& m_) : s(m_.s), m(m_.m) {}
 
    Matrix(const Matrix<T, Mupper_packed>& m_) : s(m_.s), m(&m_.m[0]) {}
 
    Matrix(const Matrix<T, Mupper_packed_ref>& m_) : s(m_.s), m(m_.m) {}
 
    bool is_null() const { return m == 0; }
 
    std::pair<size_t, size_t> size() const { return std::pair<size_t, size_t>(s, s); }
 
    size_t size1() const { return s; }
 
    size_t size2() const { return s; }
 
    size_t esize() const { return s * (s + 1) / 2; }
 
    const T& operator()(size_t i, size_t j) const {
        if (j < i) {
            return m[i * (i + 1) / 2 + j];
        } else {
            return m[j * (j + 1) / 2 + i];
        }
    }
 
    static Matrix null;
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        l << "upper packed matrix ";
        l.write(m.s * (m.s + 1) / 2, m.m, 1);
        return l;
    }
 
private:
    size_t s{};
    const T* m{};
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mupper_packed_constref> Matrix<T, Mupper_packed_constref>::null;
 
struct Mpacked {
    using base = Mpacked_ref;
    using const_base = Mpacked_st_constref;
    using copy = Mupper_packed;
};
 
template <typename T>
class Matrix<T, Mpacked> {
public:
    Matrix() = default;
 
    explicit Matrix(size_t s_, bool upper_ = false) : s(s_), m(s_ * (s_ + 1) / 2), upper(upper_) {}
 
    Matrix(const Matrix& m_) : s(m_.s), m(m_.m), upper(m_.upper) {}
 
    template <typename T1>
    Matrix(const Matrix<T1, Mpacked>& m_)
        : s(m_.s.begin(), m_.s.end()), m(m_.m.begin(), m_.m.end()), upper(m_.upper) {}
 
    template <typename T1>
    Matrix& operator=(const Matrix<T1, Mpacked>& m_) {
        s = m_.s;
        m.assign(m_.m.begin(), m_.m.end());
        upper = m_.upper;
        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 set_zero() { std::fill(m.begin(), m.end(), 0); }
 
    bool is_null() const { return this == &null; }
 
    Matrix& operator=(const Matrix& m_) {
        s = m_.s;
        m = m_.m;
        upper = m_.upper;
        return *this;
    }
 
    template <typename STORAGE>
    Matrix& operator=(const Matrix<T, STORAGE>& m_) {
        s = m_.size1();
        m.resize(s * (s + 1) / 2);
        Matrix<T, Mpacked_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>(s, s); }
 
    size_t size1() const { return s; }
 
    size_t size2() const { return s; }
 
    size_t esize() const { return s * (s + 1) / 2; }
 
    bool is_upper() const { return upper; }
 
    void resize(size_t s_, bool upper_ = false) {
        s = s_;
        m.resize(s_ * (s_ + 1) / 2);
        upper = upper_;
    }
 
    void resize(size_t s1, size_t s2, bool upper_ = false) {
        if (s1 != s2) { Exception() << THROW; }
        s = s1;
        m.resize(s1 * (s1 + 1) / 2);
        upper = upper_;
    }
 
    const T& operator()(size_t i, size_t j) const {
        if (j < i) {
            if (upper) {
                return m[i * (i + 1) / 2 + j];
            } else {
                return m[(2 * s - j - 1) * j / 2 + i];
            }
        } else if (upper) {
            return m[j * (j + 1) / 2 + i];
        } else {
            return m[(2 * s - i - 1) * i / 2 + j];
        }
    }
 
    T& operator()(size_t i, size_t j) {
        if (j < i) {
            if (upper) {
                return m[i * (i + 1) / 2 + j];
            } else {
                return m[(2 * s - j - 1) * j / 2 + i];
            }
        } else if (upper) {
            return m[j * (j + 1) / 2 + i];
        } else {
            return m[(2 * s - i - 1) * i / 2 + j];
        }
    }
 
    Matrix operator=(const Matrix<T, Sum<Mpacked_st_constref, Mpacked_st_constref> >& m_) {
        if (m_.m1.upper != m_.m2.upper) { UnimplementedError() << THROW; }
        resize(m_.size1(), m_.size2(), m_.m1.upper);
        for (size_t i = 0; i != m.size(); ++i) {
            m[i] = m_.scale * (m_.m1.scale * m_.m1.m[i] + m_.m2.scale * m_.m2.m[i]);
        }
        return *this;
    }
 
    void swap(Matrix& m_) {
        std::swap(s, m_.s);
        m.swap(m_.m);
        std::swap(upper, m_.upper);
    }
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        if (m.upper) {
            l << "upper packed matrix ";
        } else {
            l << "lower packed matrix ";
        }
        l.write(m.s * (m.s + 1) / 2, &m.m[0], 1);
        return l;
    }
 
    static Matrix null;
 
private:
    size_t s{};
    std::vector<T> m{};
    bool upper{};
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mpacked> Matrix<T, Mpacked>::null;
 
struct Mpacked_ref {
    using base = Mpacked_ref;
    using const_base = Mpacked_st_constref;
    using copy = Mlower_packed;
};
 
template <typename T>
class Matrix<T, Mpacked_ref> {
public:
    Matrix() = default;
 
    Matrix(size_t s_, T* m_, bool upper_) : s(s_), m(m_), upper(upper_) {}
 
    Matrix(const Matrix& m_) : s(m_.s), m(m_.m), upper(m_.upper) {}
 
    Matrix(Matrix<T, Mupper_packed>& m_) : s(m_.s), m(&m_.m[0]), upper(true) {}
 
    Matrix(Matrix<T, Mupper_packed_ref>& m_) : s(m_.s), m(m_.m), upper(true) {}
 
    Matrix(Matrix<T, Mlower_packed>& m_) : s(m_.s), m(&m_.m[0]), upper(false) {}
 
    Matrix(Matrix<T, Mlower_packed_ref>& m_) : s(m_.s), m(m_.m), upper(false) {}
 
    Matrix(Matrix<T, Mpacked>& m_) : s(m_.s), m(&m_.m[0]), upper(m_.upper) {}
 
    bool is_null() const { return m == 0; }
 
    void set_zero() { std::fill(m, m + esize(), 0); }
 
    Matrix& operator=(const Matrix& m_) {
        if (s != m_.s) { Exception() << "The two matrix do not have the same size." << THROW; }
        if (upper == m_.upper) {
            std::copy(m_.m, m_.m + s, m);
        } else {
            UnimplementedError() << THROW;
        }
        return *this;
    }
 
    Matrix& operator=(const Matrix<T, Mrectangle_increment_st_constref>& m_) {
        if (s != m_.s1 || s != m_.s2) {
            Exception() << "The two matrix do not have the same size." << THROW;
        }
        const T* ms1 = m_.m;
        const T* ms2 = ms1;
        T* md = m;
        const double scale = 0.5 * m_.scale;
        if (upper) {
            for (size_t j = 0; j != s; ++j, ++ms2) {
                const T* ms2_tmp = ms2;
                size_t i = 0;
                for (; i <= j; ++i, ++md, ms2_tmp += m_.ldm) { *md = scale * (ms1[i] + *ms2_tmp); }
                ms1 += m_.ldm;
            }
        } else {
            for (size_t j = 0; j != s; ++j, ms2 += m_.ldm + 1) {
                const T* ms2_tmp = ms2;
                for (size_t i = j; i != s; ++i, ++md, ms2_tmp += m_.ldm) {
                    *md = scale * (ms1[i] + *ms2_tmp);
                }
                ms1 += m_.ldm;
            }
        }
        return *this;
    }
 
    Matrix& operator=(const Matrix<T, Mpacked_st_constref>& m_) {
        if (s != m_.s) { Exception() << "The two matrix do not have the same size." << THROW; }
        if (upper == m_.upper) {
            size_t ss = s * (s + 1) / 2;
            std::copy(m_.m, m_.m + ss, m);
            if (m_.scale != T(1)) { blas::scal(ss, T(m_.scale), m, 1); }
        } else {
            UnimplementedError() << THROW;
        }
        return *this;
    }
 
    Matrix& operator=(const Matrix<
                      T, MMProd<
                               MMProd<Mrectangle_increment_st_constref, Mpacked_st_constref>,
                               Mrectangle_increment_st_constref> >& m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrix do not have the same size." << THROW;
        }
        if (m_.m1.m1.m != m_.m2.m || m_.m1.m1.s1 != m_.m2.s1 || m_.m1.m1.s2 != m_.m2.s2
            || m_.m1.m1.ldm != m_.m2.ldm || m_.m1.m1.trans == m_.m2.trans) {
            UnimplementedError() << "Not implemented (the product is not a symmetric matrix)."
                                 << THROW;
        }
 
        std::vector<T> tmp(m_.m1.m2.s);
        T* ptr = &m[0];
        const char m_m1_m2_uplo = m_.m1.m2.upper ? 'U' : 'L';
        const T alpha = m_.m1.m1.scale * m_.m1.m2.scale * m_.m2.scale * m_.scale;
        if (upper) {
            if (m_.m2.trans) {
                for (size_t j = 0; j != s; ++j) {
                    blas::spmv(
                          m_m1_m2_uplo, tmp.size(), 1, m_.m1.m2.m, &m_.m2.m[j], m_.m2.ldm, 0,
                          &tmp[0], 1);
                    blas::gemv(
                          'N', j + 1, tmp.size(), alpha, &m_.m2.m[0], m_.m2.ldm, &tmp[0], 1, T(0),
                          ptr, 1);
                    ptr += j + 1;
                }
            } else {
                for (size_t j = 0; j != s; ++j) {
                    blas::spmv(
                          m_m1_m2_uplo, tmp.size(), 1, m_.m1.m2.m, &m_.m2.m[j * m_.m2.ldm], 1, 0,
                          &tmp[0], 1);
                    blas::gemv(
                          'T', tmp.size(), j + 1, alpha, &m_.m2.m[0], m_.m2.ldm, &tmp[0], 1, T(0),
                          ptr, 1);
                    ptr += j + 1;
                }
            }
        } else if (m_.m2.trans) {
            for (size_t j = 0; j != s; ++j) {
                blas::spmv(
                      m_m1_m2_uplo, tmp.size(), 1, m_.m1.m2.m, &m_.m2.m[j], m_.m2.ldm, 0, &tmp[0],
                      1);
                blas::gemv(
                      'N', s - j, tmp.size(), alpha, &m_.m2.m[j], m_.m2.ldm, &tmp[0], 1, T(0), ptr,
                      1);
                ptr += s - j;
            }
        } else {
            for (size_t j = 0; j != s; ++j) {
                blas::spmv(
                      m_m1_m2_uplo, tmp.size(), 1, m_.m1.m2.m, &m_.m2.m[j * m_.m2.ldm], 1, 0,
                      &tmp[0], 1);
                blas::gemv(
                      'T', tmp.size(), s - j, alpha, &m_.m2.m[j * m_.m2.ldm], m_.m2.ldm, &tmp[0], 1,
                      T(0), ptr, 1);
                ptr += s - j;
            }
        }
        return *this;
    }
 
    Matrix& operator=(
          const Matrix<
                T, MMProd<
                         MMProd<Mrectangle_increment_st_constref, Msym_compressed_col_st_constref>,
                         Mrectangle_increment_st_constref> >& m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrix do not have the same size." << THROW;
        }
        if (m_.m1.m1.m != m_.m2.m || m_.m1.m1.s1 != m_.m2.s1 || m_.m1.m1.s2 != m_.m2.s2
            || m_.m1.m1.ldm != m_.m2.ldm || m_.m1.m1.trans == m_.m2.trans) {
            UnimplementedError() << "Not implemented (the product is not a symmetric matrix)."
                                 << THROW;
        }
 
        Vector<T> v1;
        T* m_ptr = m;
        const T scale = m_.scale * m_.m1.scale * m_.m1.m1.scale;
        for (size_t j = 0; j != m_.m2.s2; ++j) {
            v1 = m_.m1.m2 * m_.m2[j];
            if (upper) {
                Vector<T, Vdense_ref> v2(j + 1, m_ptr);
                v2 = scale * m_.m1.m1(Interval(0, j + 1), Interval(0, v1.size())) * v1;
                m_ptr += j + 1;
            } else {
                Vector<T, Vdense_ref> v2(s - j, m_ptr);
                v2 = scale * m_.m1.m1(Interval(j, s), Interval(0, v1.size())) * v1;
                m_ptr += s - j;
            }
        }
        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;
    }
 
    std::pair<size_t, size_t> size() const { return std::pair<size_t, size_t>(s, s); }
 
    size_t size1() const { return s; }
 
    size_t size2() const { return s; }
 
    size_t esize() const { return s * (s + 1) / 2; }
 
    const T& operator()(size_t i, size_t j) const {
        if (j < i) {
            if (upper) {
                return m[i * (i + 1) / 2 + j];
            } else {
                return m[(2 * s - j - 1) * j / 2 + i];
            }
        } else if (upper) {
            return m[j * (j + 1) / 2 + i];
        } else {
            return m[(2 * s - i - 1) * i / 2 + j];
        }
    }
 
    T& operator()(size_t i, size_t j) {
        if (j < i) {
            if (upper) {
                return m[i * (i + 1) / 2 + j];
            } else {
                return m[(2 * s - j - 1) * j / 2 + i];
            }
        } else if (upper) {
            return m[j * (j + 1) / 2 + i];
        } else {
            return m[(2 * s - i - 1) * i / 2 + j];
        }
    }
 
    static Matrix null;
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        if (m.upper) {
            l << "upper packed matrix ";
        } else {
            l << "lower packed matrix ";
        }
        l.write(m.s * (m.s + 1) / 2, m.m, 1);
        return l;
    }
 
private:
    size_t s{};
    T* m{};
    bool upper{};
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mpacked_ref> Matrix<T, Mpacked_ref>::null;
 
struct Mpacked_st_constref {
    using base = Mpacked_st_constref;
    using const_base = Mpacked_st_constref;
    using copy = Mlower_packed;
};
 
template <typename T>
class Matrix<T, Mpacked_st_constref> {
public:
    Matrix() = default;
 
    Matrix(size_t s_, const T* m_, bool scale_, bool upper_)
        : s(s_), m(m_), scale(scale_), upper(upper_) {}
 
    Matrix(const Matrix& m_) : s(m_.s), m(m_.m), scale(m_.scale), upper(m_.upper) {}
 
    Matrix(const Matrix<T, Mupper_packed>& m_) : s(m_.s), m(&m_.m[0]), scale(1), upper(true) {}
 
    Matrix(const Matrix<T, Mupper_packed_ref>& m_) : s(m_.s), m(m_.m), scale(1), upper(true) {}
 
    Matrix(const Matrix<T, Mupper_packed_constref>& m_) : s(m_.s), m(m_.m), scale(1), upper(true) {}
 
    Matrix(const Matrix<T, Mlower_packed>& m_) : s(m_.s), m(&m_.m[0]), scale(1), upper(false) {}
 
    Matrix(const Matrix<T, Mlower_packed_ref>& m_) : s(m_.s), m(m_.m), scale(1), upper(false) {}
 
    Matrix(const Matrix<T, Mlower_packed_constref>& m_)
        : s(m_.s), m(m_.m), scale(1), upper(false) {}
 
    Matrix(const Matrix<T, Mpacked>& m_) : s(m_.s), m(&m_.m[0]), scale(1), upper(m_.upper) {}
 
    Matrix(const Matrix<T, Mpacked_ref>& m_) : s(m_.s), m(m_.m), scale(1), upper(m_.upper) {}
 
    bool is_null() const { return m == 0; }
 
    std::pair<size_t, size_t> size() const { return std::pair<size_t, size_t>(s, s); }
 
    size_t size1() const { return s; }
 
    size_t size2() const { return s; }
 
    size_t esize() const { return s * (s + 1) / 2; }
 
    const T& operator()(size_t i, size_t j) const {
        if (j < i) {
            if (upper) {
                return m[i * (i + 1) / 2 + j];
            } else {
                return m[(2 * s - j - 1) * j / 2 + i];
            }
        } else if (upper) {
            return m[j * (j + 1) / 2 + i];
        } else {
            return m[(2 * s - i - 1) * i / 2 + j];
        }
    }
 
    Matrix& operator*(T scale_) {
        scale *= scale_;
        return *this;
    }
 
    static Matrix null;
 
private:
    size_t s{};
    const T* m{};
    T scale{};
    bool upper{};
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mpacked_st_constref> Matrix<T, Mpacked_st_constref>::null;
 
}  // namespace b2000::b2linalg
 
#endif