b2matrix_sym_compressed.H

//------------------------------------------------------------------------
// b2matrix_sym_compressed.H --
//
//     A framework for generic linear algebraic (matrix) computations.
//
// written by Mathias Doreille
//
// Copyright (c) 2004-2012,2016
//     SMR Engineering & Development SA
//     2502 Bienne, Switzerland
//
// 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_SYM_COMPRESSED_H_
#define B2MATRIX_SYM_COMPRESSED_H_
 
#include <ostream>
#include <unordered_map>
 
#include "b2linear_algebra_def.H"
#include "b2sparse_solver.H"
#include "b2vector_index.H"
 
namespace b2000::b2linalg {
 
struct Msym_compressed_col {
    typedef Msym_compressed_col_ref base;
    typedef Msym_compressed_col_st_constref const_base;
    typedef Msym_compressed_col copy;
    typedef Msym_compressed_col_update_inv inverse;
};
 
template <typename T>
class Matrix<T, Msym_compressed_col> {
public:
    Matrix(size_t s1_ = 0, size_t s2_ = 0, size_t snn_ = 0)
        : s1(s1_), si(s2_ + 1), m(snn_), index(snn_) {}
 
    Matrix(const Matrix& m_) : s1(m_.s1), si(m_.si), m(m_.m), index(m_.index) {}
 
    template <typename T1>
    Matrix& operator=(const Matrix<T1, Msym_compressed_col>& m_) {
        s1 = m_.s1;
        si = m_.si;
        m.assign(m_.m.begin(), m_.m.end());
        index = m_.index;
        return *this;
    }
 
    template <typename T1>
    Matrix(const Matrix<T1, Msym_compressed_col>& m_)
        : s1(m_.s1),
          si(m_.si.begin(), m_.si.end()),
          m(m_.m.begin(), m_.m.end()),
          index(m_.index.begin(), m_.index.end()) {}
 
    bool is_null() const { return this == &null; }
 
    void set_zero() {
        std::fill(si.begin(), si.end(), 0);
        index.clear();
        m.clear();
    }
 
    void resize(size_t s1_, size_t s2_) {
        s1 = s1_;
        if (size2() < s2_) {
            si.insert(si.end(), s2_ - size2(), si.back());
        } else {
            si.resize(s2_ + 1);
        }
    }
 
    void resize(size_t s1_, size_t s2_, size_t snn = 0) {
        s1 = s1_;
        si.resize(s2_ + 1);
        index.resize(snn);
        m.resize(snn);
    }
 
    void swap(std::vector<size_t>& colind, std::vector<size_t>& rowind, std::vector<T>& value) {
        s1 = colind.size() - 1;
        si.swap(colind);
        index.swap(rowind);
        m.swap(value);
    }
 
    void get_values(size_t*& colind, size_t*& rowind, T*& value) {
        colind = &si[0];
        rowind = &index[0];
        value = &m[0];
    }
 
    std::pair<size_t, size_t> size() const {
        return std::pair<size_t, size_t>(s1, si.size() == 0 ? 0 : si.size() - 1);
    }
 
    size_t size1() const { return s1; }
 
    size_t size2() const { return si.size() == 0 ? 0 : si.size() - 1; }
 
    size_t nb_nonzero() const { return si.back(); }
 
    T operator()(size_t i, size_t j) const {
        if (i < j) { std::swap(i, j); }
        std::vector<size_t>::const_iterator s = index.begin() + si[j];
        std::vector<size_t>::const_iterator e = index.begin() + si[j + 1];
        std::vector<size_t>::const_iterator ii = std::lower_bound(s, e, i);
        if (ii < e && *ii == i) { return m[ii - index.begin()]; }
        return 0;
    }
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        l << "symmetric column compressed matrix of size (" << m.s1 << ", " << m.s1 << ") ";
        l.write(m.si.size(), &m.si[0], 1, "colind");
        l.write(m.index.size(), &m.index[0], 1, "rowind");
        l.write(m.m.size(), &m.m[0], 1, "value");
        return l;
    }
 
    static Matrix null;
 
private:
    size_t s1;
    std::vector<size_t> si;
    std::vector<T> m;
    std::vector<size_t> index;
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Msym_compressed_col> Matrix<T, Msym_compressed_col>::null;
 
struct Msym_compressed_col_ref {
    typedef Msym_compressed_col_ref base;
    typedef Msym_compressed_col_st_constref const_base;
    typedef Msym_compressed_col copy;
    typedef Msym_compressed_col_update_inv inverse;
};
 
template <typename T>
class Matrix<T, Msym_compressed_col_ref> {
public:
    Matrix() : s1(0), s2(0), si(0), index(0), m(0) {}
 
    Matrix(size_t s1_, size_t s2_, size_t snn_, size_t* si_, size_t* index_, T* m_)
        : s1(s1_), s2(s2_), si(si_), index(index_), m(m_) {}
 
    Matrix(const Matrix& m_) : s1(m_.s1), s2(m_.s2), si(m_.si), index(m_.index), m(m_.m) {}
 
    Matrix(Matrix<T, Msym_compressed_col>& m_)
        : s1(m_.s1), s2(m_.si.size() - 1), si(&m_.si[0]), index(&m_.index[0]), m(&m_.m[0]) {}
 
    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; }
 
    T operator()(size_t i, size_t j) const {
        if (i < j) { std::swap(i, j); }
        const size_t* s = index + si[j];
        const size_t* e = index + si[j + 1];
        const size_t* ii = std::lower_bound(s, e, i);
        if (ii < e && *ii == i) { return m[ii - index]; }
        return 0;
    }
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        l << "symmetric column compressed matrix of size (" << m.s1 << ", " << m.s1 << ") ";
        l.write(m.s2 + 1, m.si, 1, "colind");
        l.write(m.si[m.s2], m.index, 1, "rowind");
        l.write(m.si[m.s2], m.m, 1, "value");
        return l;
    }
 
    static Matrix null;
 
private:
    size_t s1;
    size_t s2;
    size_t* si;
    size_t* index;
    T* m;
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Msym_compressed_col_ref> Matrix<T, Msym_compressed_col_ref>::null;
 
struct Msym_compressed_col_constref {
    typedef Msym_compressed_col_st_constref base;
    typedef Msym_compressed_col_st_constref const_base;
    typedef Msym_compressed_col copy;
    typedef Msym_compressed_col_update_inv inverse;
};
 
template <typename T>
class Matrix<T, Msym_compressed_col_constref> {
public:
    Matrix() : s1(0), s2(0), si(0), index(0), m(0) {}
 
    Matrix(
          size_t s1_, size_t s2_, size_t snn_, const size_t* si_, const size_t* index_, const T* m_)
        : s1(s1_), s2(s2_), si(si_), index(index_), m(m_) {}
 
    Matrix(const Matrix& m_) : s1(m_.s1), s2(m_.s2), si(m_.si), index(m_.index), m(m_.m) {}
 
    Matrix(const Matrix<T, Msym_compressed_col_ref>& m_)
        : s1(m_.s1), s2(m_.s2), si(m_.si), index(m_.index), m(m_.m) {}
 
    Matrix(const Matrix<T, Msym_compressed_col>& m_)
        : s1(m_.s1), s2(m_.si.size() - 1), si(&m_.si[0]), index(&m_.index[0]), m(&m_.m[0]) {}
 
    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; }
 
    T operator()(size_t i, size_t j) const {
        if (i < j) { std::swap(i, j); }
        const size_t* s = index + si[j];
        const size_t* e = index + si[j + 1];
        const size_t* ii = std::lower_bound(s, e, i);
        if (ii < e && *ii == i) { return m[ii - index]; }
        return 0;
    }
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        l << "symmetric column compressed matrix of size (" << m.s1 << ", " << m.s1 << ") ";
        l.write(m.s2 + 1, m.si, 1, "colind");
        l.write(m.si[m.s2], m.index, 1, "rowind");
        l.write(m.si[m.s2], m.m, 1, "value");
        return l;
    }
 
    static Matrix null;
 
private:
    size_t s1;
    size_t s2;
    const size_t* si;
    const size_t* index;
    const T* m;
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Msym_compressed_col_constref> Matrix<T, Msym_compressed_col_constref>::null;
 
struct Msym_compressed_col_st_constref {
    typedef Msym_compressed_col_st_constref base;
    typedef Msym_compressed_col_st_constref const_base;
    typedef Msym_compressed_col copy;
    typedef Msym_compressed_col_update_inv inverse;
};
 
template <typename T>
class Matrix<T, Msym_compressed_col_st_constref> {
public:
    Matrix() : s1(0), s2(0), si(0), index(0), m(0), scale(0) {}
 
    Matrix(
          size_t s1_, size_t s2_, size_t snn_, const size_t* si_, const size_t* index_, const T* m_,
          T scale_)
        : s1(s1_), s2(s2_), si(si_), index(index_), m(m_), scale(scale_) {}
 
    Matrix(const Matrix& m_)
        : s1(m_.s1), s2(m_.s2), si(m_.si), index(m_.index), m(m_.m), scale(m_.scale) {}
 
    Matrix(const Matrix<T, Msym_compressed_col_ref>& m_)
        : s1(m_.s1), s2(m_.s2), si(m_.si), index(m_.index), m(m_.m), scale(1) {}
 
    Matrix(const Matrix<T, Msym_compressed_col_constref>& m_)
        : s1(m_.s1), s2(m_.s2), si(m_.si), index(m_.index), m(m_.m), scale(1) {}
 
    Matrix(const Matrix<T, Msym_compressed_col>& m_)
        : s1(m_.s1),
          s2(m_.si.size() - 1),
          si(&m_.si[0]),
          index(&m_.index[0]),
          m(&m_.m[0]),
          scale(1) {}
 
    Matrix(Matrix<T, Msym_compressed_col_update_inv>& m_) : scale(1) {
        m_.value_to_ccarray();
        s1 = m_.si.size() - 1;
        s2 = s1;
        si = &m_.si[0];
        index = &m_.index[0];
        m = &m_.m[0];
    }
 
    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; }
 
    T operator()(size_t i, size_t j) const {
        if (i < j) { std::swap(i, j); }
        const size_t* s = index + si[j];
        const size_t* e = index + si[j + 1];
        const size_t* ii = std::lower_bound(s, e, i);
        if (ii < e && *ii == i) { return scale * m[ii - index]; }
        return 0;
    }
 
    Matrix& operator*(T scale_) {
        scale *= scale_;
        return *this;
    }
 
    static Matrix null;
 
private:
    size_t s1;
    size_t s2;
    const size_t* si;
    const size_t* index;
    const T* m;
    T scale;
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Msym_compressed_col_st_constref> Matrix<T, Msym_compressed_col_st_constref>::null;
 
struct Msym_compressed_col_update_inv {
    typedef Msym_compressed_col_ref base;
    typedef Msym_compressed_col_st_constref const_base;
    typedef Msym_compressed_col_update_inv copy;
    typedef Msym_compressed_col_update_inv inverse;
};
 
template <typename T>
class Matrix<T, Msym_compressed_col_update_inv> {
public:
    Matrix(
          size_t s1_ = 0,
          SparseMatrixConnectivityType connectivity_ = sparse_matrix_connectivity_unknown,
          const Dictionary& dictionary_ = Dictionary::get_empty())
        : value(s1_), solver(0), connectivity(connectivity_), dictionary(&dictionary_) {}
 
    Matrix(const Matrix& m_)
        : si(m_.si),
          m(m_.m),
          index(m_.index),
          value(m_.value),
          solver(0),
          connectivity(m_.connectivity),
          dictionary(m_.dictionary) {}
 
    virtual ~Matrix() { delete solver; }
 
    void resize(
          size_t s, SparseMatrixConnectivityType connectivity_ = sparse_matrix_connectivity_unknown,
          const Dictionary& dictionary_ = Dictionary::get_empty()) {
        if (si.empty()) {
            value.resize(s);
        } else {
            UnimplementedError() << THROW;
        }
        if (connectivity_ != sparse_matrix_connectivity_unknown) { connectivity = connectivity_; }
        if (&dictionary_ != &Dictionary::get_empty()) { dictionary = &dictionary_; }
    }
 
    void resize(
          size_t s1, size_t s2,
          SparseMatrixConnectivityType connectivity_ = sparse_matrix_connectivity_unknown,
          const Dictionary& dictionary_ = Dictionary::get_empty()) {
        if (s1 != s2) { Exception() << THROW; }
        resize(s1);
        if (connectivity_ != sparse_matrix_connectivity_unknown) { connectivity = connectivity_; }
        if (&dictionary_ != &Dictionary::get_empty()) { dictionary = &dictionary_; }
    }
 
    void set_same_structure(const Matrix& m_) {
        const_cast<Matrix&>(m_).value_to_ccarray();
        si = m_.si;
        index = m_.index;
        m.resize(m_.m.size());
        std::fill(m.begin(), m.end(), 0);
        connectivity = m_.connectivity;
        dictionary = m_.dictionary;
    }
 
    virtual void set_zero() {
        value.clear();
        si.clear();
        index.clear();
        m.clear();
        delete solver;
        solver = 0;
    }
 
    virtual bool is_null() const { return this == &null; }
 
    virtual void set_zero_same_struct() {
        for (typename std::vector<std::vector<std::pair<size_t, T>>>::iterator i = value.begin();
             i != value.end(); ++i) {
            i->clear();
        }
        std::fill(m.begin(), m.end(), 0);
        if (solver) { solver->update_value(); }
    }
 
    std::pair<size_t, size_t> size() const {
        size_t s;
        if (si.empty()) {
            s = value.size();
        } else {
            s = si.size() - 1;
        }
        return std::pair<size_t, size_t>(s, s);
    }
 
    size_t get_nb_nonzero() const {
        if (si.empty()) {
            size_t r = 0;
            for (size_t i = 0; i != value.size(); ++i) { r += value[i].size(); }
        }
        return index.size();
    }
 
    size_t size1() const {
        if (si.empty()) { return value.size(); }
        return si.size() - 1;
    }
 
    size_t size2() const {
        if (si.empty()) { return value.size(); }
        return si.size() - 1;
    }
 
    void InitializeRow(size_t row, const std::map<size_t, T>& row_contributions) {
        value[row].reserve(row_contributions.size());
        auto pos{begin(row_contributions)};
 
        for (; pos != end(row_contributions); pos++) {
            value[row].push_back(std::make_pair(pos->first, pos->second));
        }
    }
 
    T operator()(size_t i, size_t j) const {
        if (i < j) { std::swap(i, j); }
        if (si.empty()) {
            typename std::vector<std::pair<size_t, T>>::const_iterator ii = std::lower_bound(
                  value[j].begin(), value[j].end(), std::pair<size_t, T>(i, 0),
                  CompareFirstOfPair());
 
            if (ii != value[j].end() && ii->first == i) { return ii->second; }
 
        } else {
            std::vector<size_t>::const_iterator s = index.begin() + si[j];
            std::vector<size_t>::const_iterator e = index.begin() + si[j + 1];
            std::vector<size_t>::const_iterator ii = std::lower_bound(s, e, i);
            if (ii < e && *ii == i) { return m[ii - index.begin()]; }
        }
        return 0;
    }
 
    T& operator()(size_t i, size_t j) {
        if (i < j) { std::swap(i, j); }
 
        if (si.empty()) {  
            typename std::vector<std::pair<size_t, T>>::iterator ii = std::lower_bound(
                  value[j].begin(), value[j].end(), std::pair<size_t, T>(i, T{}),
                  CompareFirstOfPair());
            if (ii != value[j].end() && ii->first == i) { return ii->second; }
        } else {  
            std::vector<size_t>::const_iterator s = index.begin() + si[j];
            std::vector<size_t>::const_iterator e = index.begin() + si[j + 1];
            std::vector<size_t>::const_iterator ii = std::lower_bound(s, e, i);
            if (ii < e && *ii == i) { return m[ii - index.begin()]; }
        }
 
        static T zero = 0;
        return zero;
    }
 
    Matrix& operator+=(const Matrix<T, Mstructured_constref<Mpacked_st_constref>>& m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrix do not have the same size, " << size() << " and "
                        << m_.size() << THROW;
        }
        std::vector<std::pair<size_t, T>> tmp;
        tmp.reserve(m_.m.s);
        for (size_t j = 0; j != m_.m.s; ++j) {
            tmp.clear();
            size_t index_j = m_.index[j];
            for (size_t i = 0; i != m_.m.s; ++i) {
                size_t index_i = m_.index[i];
                if (index_i >= index_j) {
                    tmp.push_back(std::pair<size_t, T>(index_i, m_.m(i, j)));
                }
            }
            std::sort(tmp.begin(), tmp.end(), CompareFirstOfPair());
            add_colomn(index_j, tmp.begin(), tmp.end());
        }
        return *this;
    }
 
    Matrix& operator+=(
          const Matrix<
                T,
                MMProd<
                      MMProd<
                            Mcompressed_col_st_constref, Mstructured_constref<Mpacked_st_constref>>,
                      Mcompressed_col_st_constref>>& m_) {
        if (size1() < m_.size1()) {
            Exception() << "The two matrix do not have the same size, " << size() << " and "
                        << m_.size() << THROW;
        }
        if (m_.m1.m1.trans == true || m_.m2.trans == false || m_.m1.m1.si != m_.m2.si
            || m_.m1.m1.index != m_.m2.index || m_.m1.m1.m != m_.m2.m) {
            UnimplementedError() << THROW;
        }
 
        const size_t* colind = m_.m2.si;
        const size_t* rowind = m_.m2.index;
        const T* value_ = m_.m2.m;
        const size_t input_size = m_.m1.m2.m.s;
        const size_t* input_dof_numbering = m_.m1.m2.index;
 
        bool new_output_matrix = false;
        std::map<size_t, size_t> tmp_rowind;
        const size_t* input_dof_numbering_begin = input_dof_numbering;
        const size_t* const input_dof_numbering_end = input_dof_numbering_begin + input_size;
        while (input_dof_numbering_begin != input_dof_numbering_end) {
            const size_t* rowind_begin = rowind + colind[*input_dof_numbering_begin];
            const size_t* const rowind_end = rowind + colind[*input_dof_numbering_begin + 1];
            const T* value_begin = value_ + colind[*input_dof_numbering_begin];
            std::pair<size_t, size_t> tmp_rowind_insert(
                  0, input_dof_numbering_begin - input_dof_numbering);
            while (rowind_begin != rowind_end) {
                tmp_rowind_insert.first = *rowind_begin;
                if (!tmp_rowind.insert(tmp_rowind_insert).second || *value_begin != T(1)) {
                    new_output_matrix = true;
                }
                ++rowind_begin;
                ++value_begin;
            }
            ++input_dof_numbering_begin;
        }
        const size_t output_size = tmp_rowind.size();
        std::vector<std::pair<size_t, size_t>> output_dof_numbering(
              tmp_rowind.begin(), tmp_rowind.end());
        std::vector<std::pair<size_t, T>> output_col(output_size);
        for (size_t i = 0; i != output_size; ++i) {
            output_col[i].first = output_dof_numbering[i].first;
        }
        const Matrix<T, Mpacked_st_constref>& input_matrix = m_.m1.m2.m;
        const T scale_ = m_.scale * m_.m1.scale * m_.m1.m1.scale * m_.m2.scale;
        if (!new_output_matrix) {
            for (size_t j = 0; j != output_size; ++j) {
                size_t jj = output_dof_numbering[j].second;
                for (size_t i = j; i != output_size; ++i) {
                    output_col[i].second =
                          scale_ * input_matrix(output_dof_numbering[i].second, jj);
                }
                add_colomn(output_dof_numbering[j].first, output_col.begin() + j, output_col.end());
            }
        } else {
            {
                std::map<size_t, size_t>::iterator ii = tmp_rowind.begin();
                for (size_t i = 0; i != output_size; ++i, ++ii) { ii->second = i; }
            }
            T* output_value = TemporaryBuffer<T>::get(output_size * input_size);
            T* output_value_l = output_value;
            for (size_t j = 0; j != input_size; ++j, output_value_l += output_size) {
                for (size_t i = 0; i != input_size; ++i) {
                    const size_t* rowind_begin = rowind + colind[input_dof_numbering[i]];
                    const size_t* const rowind_end = rowind + colind[input_dof_numbering[i] + 1];
                    const T* value_begin = value_ + colind[input_dof_numbering[i]];
                    const T v = scale_ * input_matrix(j, i);
                    while (rowind_begin != rowind_end) {
                        output_value_l[tmp_rowind[*rowind_begin++]] += *value_begin++ * v;
                    }
                }
            }
            for (size_t i = 0; i != output_size; ++i) {
                for (size_t j = 0; j != output_size; ++j) { output_col[j].second = 0; }
                for (size_t j = 0; j != input_size; ++j) {
                    const size_t* rowind_begin = rowind + colind[input_dof_numbering[j]];
                    const size_t* const rowind_end = rowind + colind[input_dof_numbering[j] + 1];
                    const T* value_begin = value_ + colind[input_dof_numbering[j]];
                    const T v = output_value[i + j * output_size];
                    while (rowind_begin != rowind_end) {
                        output_col[tmp_rowind[*rowind_begin++]].second += *value_begin++ * v;
                    }
                }
                add_colomn(output_dof_numbering[i].first, output_col.begin() + i, output_col.end());
            }
            std::fill_n(output_value, output_size * input_size, 0);
        }
        return *this;
    }
 
    Matrix& operator+=(
          const Matrix<
                T, MMProd<
                         MMProd<Mcompressed_col_st_constref, Msym_compressed_col_st_constref>,
                         Mcompressed_col_st_constref>>& m_) {
        if (size() < m_.size()) {
            Exception() << "The two matrix do not have the same size, " << size() << " and "
                        << m_.size() << THROW;
        }
        if (m_.m1.m1.trans == true || m_.m2.trans == false || m_.m1.m1.si != m_.m2.si
            || m_.m1.m1.index != m_.m2.index || m_.m1.m1.m != m_.m2.m) {
            UnimplementedError() << THROW;
        }
 
        Matrix<T, Mcompressed_col> M2;
        M2 = m_.m1.m2;
        Matrix<T, Mcompressed_col> M12;
        M12 = m_.m1.scale * (m_.m1.m1 * M2);
        M2 = m_.m2;
        Matrix<T, Mcompressed_col> M123;
        M123 = m_.scale * (M12 * M2);
 
        std::vector<std::pair<size_t, T>> res_tmp;
        size_t mi = M123.si[0];
        for (size_t j = 0; j != M123.s1; ++j) {
            const size_t mi_end = M123.si[j + 1];
            for (; mi != mi_end && M123.index[mi] < j; ++mi) { ; }
            if (mi != mi_end) {
                res_tmp.resize(mi_end - mi);
                for (size_t i = 0; mi != mi_end; ++i, ++mi) {
                    res_tmp[i].first = M123.index[mi];
                    res_tmp[i].second = M123.m[mi];
                }
                add_colomn(j, res_tmp.begin(), res_tmp.end());
            }
        }
        return *this;
    }
 
    Matrix& operator+=(
          const Matrix<T, MMProd<Mcompressed_col_st_constref, Mcompressed_col_st_constref>>& m_) {
        if (size1() < m_.size1() || size2() < m_.size2()) {
            Exception() << "The two matrix do not have the same size, " << size() << " and "
                        << m_.size() << THROW;
        }
 
        if (m_.trans || m_.m1.trans || !m_.m2.trans) { UnimplementedError() << THROW; }
        if (m_.m1.m != m_.m2.m) { Exception() << THROW; }
 
        Matrix<T, Mcompressed_col> m_m2;
        m_m2 = m_.m2;
 
        size_t* tmp_index = TemporaryBuffer<size_t>::get(m_.m1.size1());
        const size_t end_list_flag = m_.m1.size1();
        T* tmp_value = TemporaryBuffer<T>::get(m_.m1.size1());
 
        T scale = m_.scale * m_.m1.scale;
        const size_t* b_colind_begin = &m_m2.si[0];
        const size_t* const b_colind_end = b_colind_begin + m_m2.size2();
        size_t end_list = end_list_flag;
        size_t col = 0;
        while (b_colind_begin != b_colind_end) {
            const size_t* b_rowind_begin = &m_m2.index[*b_colind_begin];
            const T* b_value_begin = &m_m2.m[*b_colind_begin];
            const size_t* const b_rowind_end = &m_m2.index[*++b_colind_begin];
            while (b_rowind_begin != b_rowind_end) {
                const size_t* a_rowind_begin = m_.m1.index + m_.m1.si[*b_rowind_begin];
                const size_t* const a_rowind_end = m_.m1.index + m_.m1.si[*b_rowind_begin++ + 1];
                while (a_rowind_begin != a_rowind_end && *a_rowind_begin < col) {
                    ++a_rowind_begin;
                }
                const T* a_value_begin = m_.m1.m + (a_rowind_begin - m_.m1.index);
                const T b_value_begin_v = *b_value_begin++;
                while (a_rowind_begin != a_rowind_end) {
                    tmp_value[*a_rowind_begin] += *a_value_begin++ * b_value_begin_v;
                    if (!(std::numeric_limits<size_t>::max() - tmp_index[*a_rowind_begin])) {
                        tmp_index[*a_rowind_begin] = end_list;
                        end_list = *a_rowind_begin;
                    }
                    ++a_rowind_begin;
                }
            }
            std::vector<std::pair<size_t, T>> res_tmp;
            while (end_list != end_list_flag) {
                res_tmp.push_back(std::pair<size_t, T>(end_list, scale * tmp_value[end_list]));
                const size_t end_list_next = tmp_index[end_list];
                tmp_index[end_list] = std::numeric_limits<size_t>::max();
                tmp_value[end_list] = T(0);
                end_list = end_list_next;
            }
 
            std::sort(res_tmp.begin(), res_tmp.end(), CompareFirstOfPair());
            add_colomn(col++, res_tmp.begin(), res_tmp.end());
        }
        return *this;
    }
 
    Matrix& operator+=(const Matrix& m_) {
        value_to_ccarray();
        const_cast<Matrix&>(m_).value_to_ccarray();
        if (size1() != m_.size1()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
        if (!std::equal(si.begin(), si.end(), m_.si.begin())
            || !std::equal(
                  index.begin() + si.front(), index.begin() + si.back(),
                  m_.index.begin() + m_.si.front())) {
            UnimplementedError() << THROW;
        }
        blas::axpy(m.size(), 1, &m_.m[0], 1, &m[0], 1);
        return *this;
    }
 
    Matrix& operator-=(const Matrix& m_) {
        if (size1() != m_.size1()) {
            Exception() << "The two matrices do not have the same size." << THROW;
        }
        value_to_ccarray();
        const_cast<Matrix&>(m_).value_to_ccarray();
        if (!std::equal(si.begin(), si.end(), m_.si.begin())
            || !std::equal(
                  index.begin() + si.front(), index.begin() + si.back(),
                  m_.index.begin() + m_.si.front())) {
            if (si.back() == 0) {
                si = m_.si;
                index = m_.index;
                m = m_.m;
                blas::scal(m.size(), -1, &m[0], 1);
            } else {
                UnimplementedError() << THROW;
            }
        } else {
            blas::axpy(m.size(), -1, &m_.m[0], 1, &m[0], 1);
        }
        return *this;
    }
 
    Matrix& operator+=(const Matrix<T, Msym_compressed_col_st_constref>& m_) {
        if (m_.size1() != m_.size2()) { UnimplementedError() << THROW; }
        value_to_ccarray();
        if (!std::equal(si.begin(), si.end(), m_.si)
            || !std::equal(
                  index.begin() + si[0], index.begin() + si[size2()], m_.index + m_.si[0])) {
            if (si.back() == 0) {
                si.assign(m_.si, m_.si + m_.s2 + 1);
                index.assign(m_.index, m_.index + si.back());
                m.assign(m_.m, m_.m + si.back());
                blas::scal(m.size(), m_.scale, &m[0], 1);
            } else {
                size_t s_i = m_.si[0];
                size_t d_i = si[0];
                for (size_t j = 1; j != m_.s2 + 1; ++j) {
                    const size_t s_i_end = m_.si[j];
                    const size_t d_i_end = si[j];
                    for (; d_i != d_i_end && s_i != s_i_end; ++d_i) {
                        if (m_.index[s_i] == index[d_i]) { m[d_i] += m_.scale * m_.m[s_i++]; }
                    }
                    d_i = d_i_end;
                    if (s_i != s_i_end) {
                        if (s_i_end - s_i == 1 && m_.m[s_i] == T(0)) {
                            ++s_i;
                        } else {
                            Exception() << THROW;
                        }
                    }
                }
            }
        } else {
            blas::axpy(m.size(), m_.scale, &m_.m[0], 1, &m[0], 1);
        }
        return *this;
    }
 
    Matrix& operator=(
          const Matrix<T, Sum<Msym_compressed_col_st_constref, Msym_compressed_col_st_constref>>&
                m_) {
        if (m_.size1() != m_.size2()) { UnimplementedError() << THROW; }
        if (!std::equal(m_.m1.si, m_.m1.si + m_.m1.s2 + 1, m_.m2.si)
            || !std::equal(
                  m_.m1.index + m_.m1.si[0], m_.m1.index + m_.m1.si[m_.m1.s2],
                  m_.m2.index + m_.m2.si[0])) {
            UnimplementedError() << THROW;
        }
        value.clear();
        si.clear();
        si.insert(si.begin(), m_.m1.si, m_.m1.si + m_.m1.s2 + 1);
        index.clear();
        index.insert(index.begin(), m_.m1.index + m_.m1.si[0], m_.m1.index + m_.m1.si[m_.m1.s2]);
        m.resize(index.size());
        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;
    }
 
    Matrix& operator*=(T s) {
        value_to_ccarray();
        blas::scal(m.size(), s, &m[0], 1);
        return *this;
    }
 
    template <typename STORAGE>
    void scale(const Vector<T, STORAGE>& v_) {
        value_to_ccarray();
        size_t i = si[0];
        for (size_t j = 1; j != si.size(); ++j) {
            const size_t i_end = si[j];
            for (; i != i_end; ++i) { m[i] *= v_[j] * v_[index[i]]; }
        }
    }
 
    Vector<T, Vindex1_constref> get_diagonal() {
        value_to_ccarray();
        return Vector<T, Vindex1_constref>(si.size() - 1, &m[0], si.back(), &si[0]);
    }
 
    void get_diagonal(Vector<T>& diag) {
        value_to_ccarray();
        diag.resize(si.size() - 1);
        for (size_t j = 0; j != diag.size(); ++j) { diag[j] = m[si[j]]; }
    }
 
    Matrix<T, Msym_compressed_col_update_sub_ref> operator()(const Interval& i, const Interval& j) {
        return Matrix<T, Msym_compressed_col_update_sub_ref>(*this, i, j);
    }
 
    operator const Matrix<T, Msym_compressed_col_st_constref>() const {
        const_cast<Matrix*>(this)->value_to_ccarray();
        return Matrix<T, Msym_compressed_col_st_constref>(
              si.size() - 1, si.size() - 1, m.size(), &si[0], &index[0], &m[0], 1);
    }
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        const_cast<Matrix&>(m).value_to_ccarray();
        l << "symmetric column compressed matrix of size (" << m.size1() << ", " << m.size2()
          << ") ";
        l.write(m.si.size(), &m.si[0], 1, "colind");
        l << ", ";
        l.write(m.index.size(), &m.index[0], 1, "rowind");
        l << ", ";
        l.write(m.m.size(), &m.m[0], 1, "value");
        return l;
    }
 
    friend std::ostream& operator<<(std::ostream& out, const Matrix& m) {
        const_cast<Matrix&>(m).value_to_ccarray();
        out << "{" << std::setprecision(15);
        size_t i = m.si[0];
        for (size_t j = 0; j != m.si.size() - 1; ++j) {
            for (size_t i_end = m.si[j + 1]; i != i_end; ++i) {
                if (m.m[i] != 0) { out << "(" << m.index[i] << ", " << j << "):" << m.m[i] << ","; }
            }
        }
        out << "}";
 
        /*
        std::vector<size_t> sii(m.si);
        const size_t s2 = sii.size() - 1;
        out << "(";
        for (size_t i = 0;;) {
            out << "(";
            for (size_t j = 0;;) {
                if (sii[j] == m.si[j + 1] || m.index[sii[j]] != i)
                    out << "0.0";
                else
                    out << m.m[sii[j]++];
                if (++j != s2)
                    out << ", ";
                else
                    break;
            }
            out << ")";
            if (++i != m.s1)
                out << "," << std::endl;
            else
                break;
        }
        out << ")";
        */
        return out;
    }
 
    size_t get_null_space_size() {
        if (si.size() <= 1) { return 0; }
        if (solver == 0) { return 0; }
        Matrix* noconst_this = const_cast<Matrix<T, Msym_compressed_col_update_inv>*>(this);
        return noconst_this->solver->get_null_space_size();
    }
 
    void get_null_space(Matrix<T, Mrectangle_ref> nks) {
        if (si.size() <= 1) { return; }
        if (solver == 0) { return; }
        Matrix* noconst_this = const_cast<Matrix<T, Msym_compressed_col_update_inv>*>(this);
        noconst_this->solver->get_null_space(nks.size1(), nks.size2(), nks.m, nks.size1());
    }
 
    void remove_zero(const double tol = 0) {
        value_to_ccarray();
        size_t i = si[0];
        size_t i_out = i;
        for (size_t j = 1; j != si.size(); ++j) {
            for (; i != si[j]; ++i) {
                if (b2000::norm(m[i]) > tol) {
                    m[i_out] = m[i];
                    index[i_out] = index[i];
                    ++i_out;
                }
            }
            si[j] = i_out;
        }
        index.resize(i_out);
        m.resize(i_out);
    }
 
    static Matrix null;
 
private:
    bool value_to_ccarray() {
        if (value.empty()) {
            if (si.empty()) { si.push_back(0); }
            return false;
        }
        if (si.empty()) {
            si.reserve(value.size() + 1);
            si.push_back(0);
            typename std::vector<std::vector<std::pair<size_t, T>>>::const_iterator begin =
                  value.begin();
            typename std::vector<std::vector<std::pair<size_t, T>>>::const_iterator end =
                  value.end();
            size_t nnz = 0;
            for (; begin != end; ++begin) {
                nnz += begin->size();
                if (begin->empty()) { ++nnz; }
            }
            index.reserve(nnz);
            m.reserve(nnz);
            for (begin = value.begin(); begin != end; ++begin) {
                if (begin->empty()) {
                    index.push_back(begin - value.begin());
                    m.push_back(0);
                } else {
                    typename std::vector<std::pair<size_t, T>>::const_iterator i = begin->begin();
                    typename std::vector<std::pair<size_t, T>>::const_iterator i_end = begin->end();
                    for (; i != i_end; ++i) {
                        index.push_back(i->first);
                        m.push_back(i->second);
                    }
                }
                si.push_back(m.size());
            }
        } else {
            size_t nnz = 0;
            for (size_t j = 0; j != value.size(); ++j) {
                nnz += si[j + 1] - si[j] + value[j].size();
            }
            std::vector<size_t> si_tmp;
            si_tmp.reserve(si.size());
            si_tmp.push_back(0);
            std::vector<size_t> index_tmp;
            index_tmp.reserve(nnz);
            std::vector<T> m_tmp;
            m_tmp.reserve(nnz);
            for (size_t j = 0; j != value.size(); ++j) {
                size_t i = si[j];
                const size_t i_end = si[j + 1];
                if (!value[j].empty()) {
                    typename std::vector<std::pair<size_t, T>>::const_iterator iv =
                          value[j].begin();
                    typename std::vector<std::pair<size_t, T>>::const_iterator iv_end =
                          value[j].end();
                    for (;;) {
                        if (i == i_end || iv->first < index[i]) {
                            index_tmp.push_back(iv->first);
                            m_tmp.push_back(iv->second);
                            if (++iv == iv_end) { break; }
                        } else {
                            index_tmp.push_back(index[i]);
                            m_tmp.push_back(m[i]);
                            ++i;
                        }
                    }
                }
                index_tmp.insert(index_tmp.end(), index.begin() + i, index.begin() + i_end);
                m_tmp.insert(m_tmp.end(), m.begin() + i, m.begin() + i_end);
                si_tmp.push_back(m_tmp.size());
            }
            si.swap(si_tmp);
            index.swap(index_tmp);
            m.swap(m_tmp);
            delete solver;
            solver = 0;
            value.clear();
            return true;
        }
        value.clear();
        return false;
    }
 
    void resolve(
          size_t s, size_t nrhs, const T* b, size_t ldb, T* x, size_t ldx, char left_or_right,
          Matrix<T, Mrectangle>& null_space, ssize_t max_null_space_vector) const {
        if (s == 0) { return; }
        Matrix* noconst_this = const_cast<Matrix<T, Msym_compressed_col_update_inv>*>(this);
        noconst_this->value_to_ccarray();
        try {
            if (solver == 0) {
                noconst_this->solver =
                      LDLt_sparse_solver<T>::new_default(connectivity, *dictionary);
                noconst_this->solver->init(
                      si.size() - 1, index.size(), &si[0], &index[0], &m[0], connectivity,
                      *dictionary);
            }
            noconst_this->solver->resolve(s, nrhs, b, ldb, x, ldx, left_or_right);
        } catch (SingularMatrixError& e) {
            if (max_null_space_vector != 0) {
                const size_t nbnv =
                      max_null_space_vector > 0
                                  && ssize_t(e.null_space_size) > max_null_space_vector
                            ? max_null_space_vector
                            : e.null_space_size;
                null_space.resize(s, nbnv);
                noconst_this->solver->get_null_space(s, nbnv, &null_space(0, 0), s);
            }
            throw e;
        }
        const size_t nss = noconst_this->solver->get_null_space_size();
        if (max_null_space_vector != 0 && nss > 0) {
            const size_t nbnv = max_null_space_vector > 0 && ssize_t(nss) > max_null_space_vector
                                      ? max_null_space_vector
                                      : nss;
            null_space.resize(s, nbnv);
            noconst_this->solver->get_null_space(s, nbnv, &null_space(0, 0), s);
        }
    }
 
    void add_colomn(
          size_t col, typename std::vector<std::pair<size_t, T>>::const_iterator begin,
          typename std::vector<std::pair<size_t, T>>::const_iterator end) {
        if (begin == end) { return; }
 
        if (si.empty()) {
            std::vector<std::pair<size_t, T>>& vcol = value[col];
            if (vcol.empty() || begin->first > vcol.back().first) {
                vcol.insert(vcol.end(), begin, end);
            } else {
                std::vector<std::pair<size_t, T>> tmp;
                tmp.reserve(vcol.size() + (end - begin));
                typename std::vector<std::pair<size_t, T>>::const_iterator vbegin = vcol.begin();
                typename std::vector<std::pair<size_t, T>>::const_iterator vend = vcol.end();
                while (vbegin != vend && begin != end) {
                    if (vbegin->first < begin->first) {
                        tmp.push_back(*vbegin++);
                    } else if (vbegin->first > begin->first) {
                        tmp.push_back(*begin++);
                    } else {
                        tmp.push_back(
                              std::pair<size_t, T>(vbegin->first, vbegin->second + begin->second));
                        ++begin;
                        ++vbegin;
                    }
                }
                tmp.insert(tmp.end(), vbegin, vend);
                tmp.insert(tmp.end(), begin, end);
                vcol.swap(tmp);
            }
        } else {
            size_t colind = si[col];
            T* beginv = &m[colind];
            size_t* begini = &index[colind];
            size_t* begini_end = &index[si[col + 1]];
            while (begini != begini_end) {
                if (*begini == begin->first) {
                    *beginv += begin->second;
                    if (++begin == end) { break; }
                }
                ++begini;
                ++beginv;
            }
            if (begin == end) { return; }
            {
                if (value.empty()) { value.resize(si.size() - 1); }
                std::vector<std::pair<size_t, T>>& vcol = value[col];
                beginv = &m[colind];
                begini = &index[colind];
                if (vcol.empty() || begin->first > vcol.back().first) {
                    while (begin != end) {
                        if (begini == begini_end || *begini > begin->first) {
                            vcol.push_back(*begin++);
                        } else {
                            if (*begini == begin->first) { *beginv += (begin++)->second; }
                            ++begini;
                            ++beginv;
                        }
                    }
                } else {
                    std::vector<std::pair<size_t, T>> tmp;
                    tmp.reserve(vcol.size() + (end - begin));
                    typename std::vector<std::pair<size_t, T>>::const_iterator vbegin =
                          vcol.begin();
                    typename std::vector<std::pair<size_t, T>>::const_iterator vend = vcol.end();
                    while (begin != end) {
                        if (begini == begini_end || *begini > begin->first) {
                            while (vbegin != vend && vbegin->first < begin->first) {
                                tmp.push_back(*vbegin++);
                            }
                            if (vbegin != vend) {
                                if (vbegin->first > begin->first) {
                                    tmp.push_back(*begin);
                                } else {
                                    tmp.push_back(std::pair<size_t, T>(
                                          vbegin->first, vbegin->second + begin->second));
                                    ++vbegin;
                                }
                            } else {
                                tmp.push_back(*begin);
                            }
                            ++begin;
                        } else {
                            if (*begini == begin->first) { *beginv += (begin++)->second; }
                            ++begini;
                            ++beginv;
                        }
                    }
                    tmp.insert(tmp.end(), vbegin, vend);
                    vcol.swap(tmp);
                }
            }
        }
    }
 
    std::vector<size_t> si;     
    std::vector<T> m;           
    std::vector<size_t> index;  
    std::vector<std::vector<std::pair<size_t, T>>>
          value;  
    LDLt_sparse_solver<T>* solver;
    SparseMatrixConnectivityType connectivity;
    const Dictionary* dictionary;
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Msym_compressed_col_update_inv> Matrix<T, Msym_compressed_col_update_inv>::null;
 
struct Msym_compressed_col_update_inv_ext {
    typedef Msym_compressed_col_update_inv_ext const_base;
    typedef Msym_compressed_col_update_inv_ext copy;
    typedef Msym_compressed_col_update_inv_ext inverse;
};
 
template <typename T>
class Matrix<T, Msym_compressed_col_update_inv_ext>
    : public Matrix<T, Msym_compressed_col_update_inv> {
public:
    Matrix(
          size_t s1_ = 0, size_t size_ext_ = 0,
          SparseMatrixConnectivityType connectivity_ = sparse_matrix_connectivity_unknown,
          const Dictionary& dictionary_ = Dictionary::get_empty())
        : Matrix<T, Msym_compressed_col_update_inv>(s1_, connectivity_, dictionary_),
          size_ext(size_ext_),
          m_ab(2 * s1_ * size_ext_) {
        if (size_ext_ > 1) { UnimplementedError() << THROW; }
    }
 
    Matrix(const Matrix& m_) : size_ext(m_.size_ext) {}
 
    std::pair<size_t, size_t> size() const {
        return std::pair<size_t, size_t>(
              Matrix<T, Msym_compressed_col_update_inv>::size1() + size_ext,
              Matrix<T, Msym_compressed_col_update_inv>::size2() + size_ext);
    }
 
    size_t size1() const { return Matrix<T, Msym_compressed_col_update_inv>::size1() + size_ext; }
 
    size_t size2() const { return Matrix<T, Msym_compressed_col_update_inv>::size2() + size_ext; }
 
    void resize(
          size_t s, size_t s_ext,
          SparseMatrixConnectivityType connectivity_ = sparse_matrix_connectivity_unknown,
          const Dictionary& dictionary_ = Dictionary::get_empty()) {
        Matrix<T, Msym_compressed_col_update_inv>::resize(s, connectivity_, dictionary_);
        size_ext = s_ext;
        m_ab.resize(s * 2 * size_ext);
        if (size_ext != 1) { UnimplementedError() << THROW; }
    }
 
    void set_zero() { Matrix<T, Msym_compressed_col_update_inv>::set_zero(); }
 
    void set_zero_same_struct() {
        Matrix<T, Msym_compressed_col_update_inv>::set_zero_same_struct();
    }
 
    bool is_null() const { return this == &null; }
 
    T operator()(size_t i, size_t j) const {
        if (i < Matrix<T, Msym_compressed_col_update_inv>::size1()
            && j < Matrix<T, Msym_compressed_col_update_inv>::size2()) {
            return Matrix<T, Msym_compressed_col_update_inv>::operator()(i, j);
        }
        return 0;
    }
 
    Matrix<T, Msym_compressed_col_update_sub_ref> operator()(const Interval& i, const Interval& j) {
        return Matrix<T, Msym_compressed_col_update_sub_ref>(*this, i, j);
    }
 
    void resolve(
          size_t s, size_t nrhs, const T* b, size_t ldb, T* x, size_t ldx, const T* ma, const T* mb,
          const T* mc, char left_or_right, Matrix<T, Mrectangle>& null_space,
          ssize_t max_null_space_vector) const {
        Matrix* noconst_this = const_cast<Matrix<T, Msym_compressed_col_update_inv_ext>*>(this);
        if (noconst_this->value_to_ccarray()) {
            delete noconst_this->solver;
            noconst_this->solver = 0;
        }
        try {
            if (this->solver == 0) {
                noconst_this->solver = LDLt_sparse_solver<T>::new_default(
                      Matrix<T, Msym_compressed_col_update_inv>::connectivity,
                      *Matrix<T, Msym_compressed_col_update_inv>::dictionary);
                noconst_this->solver->init(
                      Matrix<T, Msym_compressed_col_update_inv>::si.size() - 1,
                      Matrix<T, Msym_compressed_col_update_inv>::index.size(),
                      &Matrix<T, Msym_compressed_col_update_inv>::si[0],
                      &Matrix<T, Msym_compressed_col_update_inv>::index[0],
                      &Matrix<T, Msym_compressed_col_update_inv>::m[0],
                      Matrix<T, Msym_compressed_col_update_inv>::connectivity,
                      *Matrix<T, Msym_compressed_col_update_inv>::dictionary);
            }
 
            if (mc != 0) {
                std::copy(ma, ma + s - 1, const_cast<T*>(&m_ab[0]));
                std::copy(mb, mb + s - 1, const_cast<T*>(&m_ab[s - 1]));
                noconst_this->solver->resolve(
                      s - 1, 2, &m_ab[0], s - 1, const_cast<T*>(&m_ab[0]), s - 1);
                const_cast<T&>(div) = 1 / (*mc - blas::dot(s - 1, ma, 1, &m_ab[s - 1], 1));
            }
 
            for (size_t i = 0; i != nrhs; ++i) {
                const double x2 = x[ldx * i + s - 1] =
                      div
                      * (b[ldb * i + s - 1] - blas::dot(s - 1, b + ldb * i, 1, &m_ab[s - 1], 1));
                noconst_this->solver->resolve(s - 1, 1, b + ldb * i, ldb, x + ldx * i, ldx);
                blas::axpy(s - 1, -x2, &m_ab[0], 1, x + ldx * i, 1);
            }
 
            // noconst_this->solver->resolve(s, nrhs, b, ldb, x, ldx,
            //    ma, mb, mc, ' ');
        } catch (SingularMatrixError& e) {
            if (max_null_space_vector != 0) {
                const size_t nbnv =
                      max_null_space_vector > 0
                                  && ssize_t(e.null_space_size) > max_null_space_vector
                            ? max_null_space_vector
                            : e.null_space_size;
                null_space.resize(s, nbnv);
                noconst_this->solver->get_null_space(s, nbnv, &null_space(0, 0), s);
            }
            throw;
        }
        const size_t nss = noconst_this->solver->get_null_space_size();
        if (max_null_space_vector != 0 && nss > 0) {
            const size_t nbnv = max_null_space_vector > 0 && ssize_t(nss) > max_null_space_vector
                                      ? max_null_space_vector
                                      : nss;
            null_space.resize(s, nbnv);
            noconst_this->solver->get_null_space(s, nbnv, &null_space(0, 0), s);
        }
    }
 
    size_t get_null_space_size() {
        if (Matrix<T, Msym_compressed_col_update_inv>::si.size() <= 1) { return 0; }
        if (this->solver == 0) { return 0; }
        Matrix* noconst_this = const_cast<Matrix<T, Msym_compressed_col_update_inv_ext>*>(this);
        return noconst_this->solver->get_null_space_size();
    }
 
    void get_null_space(Matrix<T, Mrectangle_ref> nks) {
        if (Matrix<T, Msym_compressed_col_update_inv>::si.size() <= 1) { return; }
        if (this->solver == 0) { return; }
        Matrix* noconst_this = const_cast<Matrix<T, Msym_compressed_col_update_inv_ext>*>(this);
        noconst_this->solver->get_null_space(nks.size1(), nks.size2(), nks.m, nks.size1());
    }
 
    static Matrix null;
 
private:
    size_t size_ext;
    std::vector<double> m_ab;
    double div;
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Msym_compressed_col_update_inv_ext> Matrix<T, Msym_compressed_col_update_inv_ext>::null;
 
struct Msym_compressed_col_update_sub_ref {
    typedef Mlower_packed copy;
};
 
template <typename T>
class Matrix<T, Msym_compressed_col_update_sub_ref> {
public:
    Matrix(Matrix<T, Msym_compressed_col_update_inv>& m_, const Interval& i_, const Interval& j_)
        : m(m_), i(i_), j(j_) {}
 
    std::pair<size_t, size_t> size() const { return std::pair<size_t, size_t>(i.size(), j.size()); }
 
    size_t size1() const { return i.size(); }
 
    size_t size2() const { return j.size(); }
 
    Matrix& operator+=(const Matrix<T, Mcompressed_col_constref>& m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrix do not have the same size, " << size() << " and "
                        << m_.size() << "." << THROW;
        }
 
        if (i.end <= j.start) {
            // Only the lower part of the matrix is stored and must by updated
            return *this;
        }
 
        std::vector<std::pair<size_t, T>> res;
        size_t ii = m_.si[0];
        for (size_t jj = 0; jj != m_.size2(); ++jj) {
            const size_t ii_end = m_.si[jj + 1];
            res.clear();
            res.reserve(ii_end - ii);
            for (; ii != ii_end; ++ii) {
                if (i[0] + m_.index[ii] >= j[0] + jj) {
                    res.push_back(std::pair<size_t, T>(i[0] + m_.index[ii], m_.m[ii]));
                }
            }
            m.add_colomn(j[0] + jj, res.begin(), res.end());
        }
        return *this;
    }
 
    Matrix& operator+=(
          const Matrix<T, MMProd<Mcompressed_col_st_constref, Mcompressed_col_st_constref>>& m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrix do not have the same size, " << size() << " and "
                        << m_.size() << "." << THROW;
        }
 
        if (i.end < j.start) {
            // Only the lower part of the matrix is stored and must by updated
            return *this;
        }
 
        if (!(m_.trans || m_.m1.trans || m_.m2.trans)) {
            size_t* tmp_index = TemporaryBuffer<size_t>::get(m_.m1.size1());
            const size_t end_list_flag = m_.m1.size1();
            T* tmp_value = TemporaryBuffer<T>::get(m_.m1.size1());
 
            const size_t* b_colind_begin = m_.m2.si;
            const size_t* const b_colind_end = b_colind_begin + m_.m2.size2();
            size_t end_list = end_list_flag;
            size_t col = j[0];
            while (b_colind_begin != b_colind_end) {
                const size_t* b_rowind_begin = &m_.m2.index[*b_colind_begin];
                const T* b_value_begin = &m_.m2.m[*b_colind_begin];
                const size_t* const b_rowind_end = &m_.m2.index[*++b_colind_begin];
                while (b_rowind_begin != b_rowind_end) {
                    const size_t* a_rowind_begin = m_.m1.index + m_.m1.si[*b_rowind_begin];
                    const size_t* const a_rowind_end = m_.m1.index + m_.m1.si[*b_rowind_begin + 1];
                    const T* a_value_begin = m_.m1.m + m_.m1.si[*b_rowind_begin++];
                    const T b_value_begin_v = *b_value_begin++;
                    while (a_rowind_begin != a_rowind_end) {
                        tmp_value[*a_rowind_begin] += *a_value_begin++ * b_value_begin_v;
                        if (!(std::numeric_limits<size_t>::max() - tmp_index[*a_rowind_begin])) {
                            tmp_index[*a_rowind_begin] = end_list;
                            end_list = *a_rowind_begin;
                        }
                        ++a_rowind_begin;
                    }
                }
                std::vector<std::pair<size_t, T>> res_tmp;
                while (end_list != end_list_flag) {
                    res_tmp.push_back(std::pair<size_t, T>(i[end_list], tmp_value[end_list]));
                    const size_t end_list_next = tmp_index[end_list];
                    tmp_index[end_list] = std::numeric_limits<size_t>::max();
                    tmp_value[end_list] = T(0);
                    end_list = end_list_next;
                }
                std::sort(res_tmp.begin(), res_tmp.end(), CompareFirstOfPair());
                m.add_colomn(col++, res_tmp.begin(), res_tmp.end());
            }
        } else if (!m_.trans && m_.m1.trans && m_.m2.trans) {
            std::vector<std::vector<std::pair<size_t, T>>> res_tmp(j.size());
 
            size_t* tmp_index = TemporaryBuffer<size_t>::get(m_.m2.size2());
            const size_t end_list_flag = m_.m2.size2();
            T* tmp_value = TemporaryBuffer<T>::get(m_.m2.size2());
 
            const size_t* b_colind_begin = m_.m1.si;
            const size_t* const b_colind_end = b_colind_begin + m_.m1.size1();
            size_t end_list = end_list_flag;
            size_t col = 0;
            while (b_colind_begin != b_colind_end) {
                const size_t* b_rowind_begin = &m_.m1.index[*b_colind_begin];
                const T* b_value_begin = &m_.m1.m[*b_colind_begin];
                const size_t* const b_rowind_end = &m_.m1.index[*++b_colind_begin];
                while (b_rowind_begin != b_rowind_end) {
                    const size_t* a_rowind_begin = m_.m2.index + m_.m2.si[*b_rowind_begin];
                    const size_t* const a_rowind_end = m_.m2.index + m_.m2.si[*b_rowind_begin + 1];
                    const T* a_value_begin = m_.m2.m + m_.m2.si[*b_rowind_begin++];
                    const T b_value_begin_v = *b_value_begin++;
                    while (a_rowind_begin != a_rowind_end) {
                        tmp_value[*a_rowind_begin] += *a_value_begin++ * b_value_begin_v;
                        if (!(std::numeric_limits<size_t>::max() - tmp_index[*a_rowind_begin])) {
                            tmp_index[*a_rowind_begin] = end_list;
                            end_list = *a_rowind_begin;
                        }
                        ++a_rowind_begin;
                    }
                }
                while (end_list != end_list_flag) {
                    res_tmp[end_list].push_back(std::pair<size_t, T>(i[col], tmp_value[end_list]));
                    const size_t end_list_next = tmp_index[end_list];
                    tmp_index[end_list] = std::numeric_limits<size_t>::max();
                    tmp_value[end_list] = T(0);
                    end_list = end_list_next;
                }
                ++col;
            }
            for (size_t jj = 0; jj != res_tmp.size(); ++jj) {
                if (!res_tmp[jj].empty()) {
                    std::sort(res_tmp[jj].begin(), res_tmp[jj].end(), CompareFirstOfPair());
                    m.add_colomn(j[jj], res_tmp[jj].begin(), res_tmp[jj].end());
                }
            }
        } else if (!m_.trans && !m_.m1.trans && m_.m2.trans && m_.m1.m == m_.m2.m) {
            Matrix<T, Mcompressed_col> m_m2;
            m_m2 = m_.m2;
 
            size_t* tmp_index = TemporaryBuffer<size_t>::get(m_.m1.size1());
            const size_t end_list_flag = m_.m1.size1();
            T* tmp_value = TemporaryBuffer<T>::get(m_.m1.size1());
 
            T scale = m_.scale * m_.m1.scale;
            const size_t* b_colind_begin = &m_m2.si[0];
            const size_t* const b_colind_end = b_colind_begin + m_m2.size2();
            size_t end_list = end_list_flag;
            size_t col = 0;
            while (b_colind_begin != b_colind_end) {
                const size_t* b_rowind_begin = &m_m2.index[*b_colind_begin];
                const T* b_value_begin = &m_m2.m[*b_colind_begin];
                const size_t* const b_rowind_end = &m_m2.index[*++b_colind_begin];
                while (b_rowind_begin != b_rowind_end) {
                    const size_t* a_rowind_begin = m_.m1.index + m_.m1.si[*b_rowind_begin];
                    const size_t* const a_rowind_end =
                          m_.m1.index + m_.m1.si[*b_rowind_begin++ + 1];
                    while (a_rowind_begin != a_rowind_end && *a_rowind_begin < col) {
                        ++a_rowind_begin;
                    }
                    const T* a_value_begin = m_.m1.m + (a_rowind_begin - m_.m1.index);
                    const T b_value_begin_v = *b_value_begin++;
                    while (a_rowind_begin != a_rowind_end) {
                        tmp_value[*a_rowind_begin] += *a_value_begin++ * b_value_begin_v;
                        if (!(std::numeric_limits<size_t>::max() - tmp_index[*a_rowind_begin])) {
                            tmp_index[*a_rowind_begin] = end_list;
                            end_list = *a_rowind_begin;
                        }
                        ++a_rowind_begin;
                    }
                }
                std::vector<std::pair<size_t, T>> res_tmp;
                while (end_list != end_list_flag) {
                    res_tmp.push_back(
                          std::pair<size_t, T>(i[end_list], scale * tmp_value[end_list]));
                    const size_t end_list_next = tmp_index[end_list];
                    tmp_index[end_list] = std::numeric_limits<size_t>::max();
                    tmp_value[end_list] = T(0);
                    end_list = end_list_next;
                }
 
                std::sort(res_tmp.begin(), res_tmp.end(), CompareFirstOfPair());
                m.add_colomn(j[col++], res_tmp.begin(), res_tmp.end());
            }
        } else {
            UnimplementedError() << THROW;
        }
 
        return *this;
    }
 
private:
    Matrix<T, Msym_compressed_col_update_inv>& m;
    const Interval& i;
    const Interval& j;
    MVFRIEND;
};
 
}  // namespace b2000::b2linalg
 
#endif