b2matrix_compressed.H

//------------------------------------------------------------------------
// b2matrix_compressed.H --
//
//     A framework for generic linear algebraic (matrix) computations.
//
// written by Mathias Doreille
//            Thomas Blome <thomas.blome@dlr.de>
//
// Copyright (c) 2004-2012,2016,2017 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_COMPRESSED_H__
#define __B2MATRIX_COMPRESSED_H__
 
#include <cassert>
#include <map>
#include <set>
#include <vector>
 
#include "b2linear_algebra_def.H"
#include "b2sparse_solver.H"
#include "b2vector_compressed.H"
#include "b2vector_dense.H"
#include "b2vector_index.H"
#include "utils/b2dictionary.H"
#include "utils/b2util.H"
 
namespace b2000::b2linalg {
 
struct Mcompressed_col {
    using base = Mcompressed_col_ref;
    using const_base = Mcompressed_col_st_constref;
    using copy = Mcompressed_col;
};
 
template <typename T>
class Matrix<T, Mcompressed_col> {
public:
    Matrix() : s1(0) { si.assign(1, 0); }
 
    Matrix(size_t s1_, size_t s2_, size_t snn_) : 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(const Matrix<T1, Mcompressed_col>& m_, bool set_zero_same_structure = false)
        : s1(m_.s1), si(m_.si), m(m_.m.begin(), m_.m.end()), index(m_.index) {}
 
    template <typename T1>
    Matrix& operator=(const Matrix<T1, Mcompressed_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& operator=(const Matrix<T1, Msym_compressed_col_st_constref>& m_) {
        s1 = m_.s1;
        si.assign(s1 + 2, 0);
 
        const size_t* m_index = m_.index + m_.si[0];
        size_t* si1 = &si[0] + 2;
        for (size_t j = 0; j != s1; ++j) {
            const size_t* m_index_end = m_.index + m_.si[j + 1];
            si1[j] += m_index_end - m_index;
            if (m_index != m_index_end && *m_index == j) { ++m_index; }
            for (; m_index < m_index_end; ++m_index) { ++si1[*m_index]; }
        }
 
        std::partial_sum(si1, si1 + s1, si1);
 
        --si1;
        index.assign(si1[s1], 0);
        m.assign(si1[s1], T(0));
        size_t i = m_.si[0];
        for (size_t j = 0; j != s1; ++j) {
            const size_t i_end = m_.si[j + 1];
            std::copy(m_.index + i, m_.index + i_end, index.begin() + si1[j]);
            {
                const T* b = m_.m + i;
                const T* b_end = m_.m + i_end;
                T* bs = &m[si1[j]];
                for (; b != b_end; ++b, ++bs) { *bs = m_.scale * *b; }
            }
            si1[j] += i_end - i;
            if (i != i_end && m_.index[i] == j) { ++i; }
            for (; i < i_end; ++i) {
                size_t& ii = si1[m_.index[i]];
                index[ii] = j;
                m[ii] = m_.scale * m_.m[i];
                ++ii;
            }
        }
 
        si.pop_back();
        return *this;
    }
 
    Matrix(const Matrix<T, Mcompressed_col_st_constref>& m_) { *this = m_; }
 
    Matrix(const Matrix<T, Msym_compressed_col_st_constref>& m_, const Index& index_) {
        resize(m_.size1(), index_.size(), 0);
        std::vector<size_t> ind(index_);
        std::sort(ind.begin(), ind.end());
        const size_t s2 = index_.size();
        size_t ij = 0;
        for (size_t j = 0; j != m_.s1; ++j) {
            size_t ii = m_.si[j];
            const size_t ii_end = m_.si[j + 1];
            if (ii == ii_end) { continue; }
            if (ij != s2 && j == ind[ij]) {
                si[ij] += ii_end - ii;
                ++ij;
            }
            if (m_.index[ii] == j) { ++ii; }
            size_t i = 0;
            while (i != s2 && ii != ii_end) {
                if (m_.index[ii] < ind[i]) {
                    ++ii;
                } else if (m_.index[ii] > ind[i]) {
                    ++i;
                } else {
                    ++si[i];
                    ++ii;
                    ++i;
                }
            }
        }
        for (size_t i = 0; i != s2; ++i) { si[i + 1] += si[i]; }
        index.assign(si[s2], 0);
        m.assign(si[s2], 0);
        for (size_t i = s2; i > 1; --i) { si[i] = si[i - 2]; }
        si[0] = si[1] = 0;
        size_t* siptr = &si[1];
        ij = 0;
        for (size_t j = 0; j != m_.s1; ++j) {
            size_t ii = m_.si[j];
            const size_t ii_end = m_.si[j + 1];
            if (ii == ii_end) { continue; }
            if (ij != s2 && j == ind[ij]) {
                std::copy(m_.index + ii, m_.index + ii_end, &index[siptr[ij]]);
                std::copy(m_.m + ii, m_.m + ii_end, &m[siptr[ij]]);
                siptr[ij] += ii_end - ii;
                ++ij;
            }
            if (m_.index[ii] == j) { ++ii; }
            size_t i = 0;
            while (i != s2 && ii != ii_end) {
                if (m_.index[ii] < ind[i]) {
                    ++ii;
                } else if (m_.index[ii] > ind[i]) {
                    ++i;
                } else {
                    index[siptr[i]] = j;
                    m[siptr[i]] = m_.m[ii];
                    ++siptr[i];
                    ++ii;
                    ++i;
                }
            }
        }
    }
 
    void this_prod_index(Index& i_index) const {
        const size_t* colind = &si[0];
        const size_t* rowind = &index[0];
        std::set<size_t> tmp_rowind;
        const size_t* const i_end = &i_index[0] + i_index.size();
        const size_t i_max = si.size() - 1;
        for (const size_t* i = &i_index[0]; i != i_end; ++i) {
            if (*i >= i_max) { Exception() << THROW; }
            tmp_rowind.insert(rowind + colind[*i], rowind + colind[*i + 1]);
        }
        i_index.assign(tmp_rowind.begin(), tmp_rowind.end());
    }
 
    void remove_row(const Index& index_row) {
        std::vector<size_t> tmp_index(s1);
        {
            size_t ii = 0;
            size_t i_index = 0;
            for (size_t i = 0; i != s1; ++i) {
                if (i_index != index_row.size() && i == index_row[i_index]) {
                    tmp_index[i] = s1;
                    ++i_index;
                } else {
                    tmp_index[i] = ii++;
                }
            }
        }
        size_t i_o = si[0];
        si[0] = 0;
        size_t i_n = 0;
        for (size_t j = 1; j != si.size(); ++j) {
            const size_t i_oe = si[j];
            for (; i_o != i_oe; ++i_o) {
                const size_t tmp_i = tmp_index[index[i_o]];
                if (tmp_i != s1) {
                    index[i_n] = tmp_i;
                    m[i_n] = m[i_o];
                    ++i_n;
                }
            }
            si[j] = i_n;
        }
        index.resize(i_n);
        m.resize(i_n);
        s1 -= index_row.size();
    }
 
    void remove_col(const Index& index_col) {
        size_t index_col_ptr = 0;
        size_t si_dest = 0;
        size_t dest = 0;
        for (size_t j = 0; j != si.size() - 1; ++j) {
            if (index_col_ptr != index_col.size() && j == index_col[index_col_ptr]) {
                ++index_col_ptr;
            } else {
                if (dest < si[j]) {
                    std::copy(
                          index.begin() + si[j], index.begin() + si[j + 1], index.begin() + dest);
                    std::copy(m.begin() + si[j], m.begin() + si[j + 1], m.begin() + dest);
                }
                const size_t s = si[j + 1] - si[j];
                si[si_dest++] = dest;
                dest += s;
            }
        }
        si[si_dest++] = dest;
        si.resize(si_dest);
        index.resize(dest);
        m.resize(dest);
    }
 
    void remove_nonzero_in_cols(const Index& index_col) {
        size_t i_o = si[0];
        size_t i_n = i_o;
        size_t jj = 0;
        for (size_t j = 0; j != index_col.size(); ++j) {
            const size_t i_l = si[index_col[j]];
            if (i_n != i_o) {
                std::copy(index.begin() + i_o, index.begin() + i_l, index.begin() + i_n);
                std::copy(m.begin() + i_o, m.begin() + i_l, m.begin() + i_n);
                const size_t diff = i_o - i_n;
                for (; jj != index_col[j]; ++jj) { si[jj] -= diff; }
            }
            i_n += i_l - i_o;
            i_o = si[index_col[j] + 1];
        }
        if (i_n != i_o) {
            const size_t i_l = si.back();
            std::copy(index.begin() + i_o, index.begin() + i_l, index.begin() + i_n);
            std::copy(m.begin() + i_o, m.begin() + i_l, m.begin() + i_n);
            const size_t diff = i_o - i_n;
            for (; jj != si.size(); ++jj) { si[jj] -= diff; }
        }
    }
 
    bool get_nonzero_unit_row_for_column(const Index& index_col, Index& index_row) {
        index_row.clear();
        index_row.reserve(index_col.size());
        for (size_t j = 0; j != index_col.size(); ++j) {
            size_t jj = index_col[j];
            if (si[jj + 1] - si[jj] != 1 || m[si[jj]] != T(1)) { return false; }
            index_row.push_back(index[si[jj]]);
        }
        std::sort(index_row.begin(), index_row.end());
        return true;
    }
 
    Matrix& operator=(const Matrix<T, Mcompressed_col_st_constref>& m_) {
        s1 = m_.size1();
        if (!m_.trans) {
            si.assign(m_.si, m_.si + m_.s2 + 1);
            index.assign(m_.index, m_.index + m_.si[m_.s2]);
            if (m_.scale == T(1)) {
                m.assign(m_.m, m_.m + m_.si[m_.s2]);
            } else {
                m.reserve(m_.si[m_.s2]);
                for (const T* i = m_.m; i != m_.m + m_.si[m_.s2]; ++i) {
                    m.push_back(m_.scale * *i);
                }
            }
        } else {
            si.assign(m_.s1 + 1, 0);
            m.assign(m_.si[m_.s2], 0);
            index.assign(m_.si[m_.s2], 0);
 
            size_t* tmp_col = new size_t[m_.s1];
            std::fill_n(tmp_col, m_.s1, 0);
 
            const size_t* colind_begin = m_.si;
            const size_t* rowind_begin = m_.index + *colind_begin;
            const size_t* const colind_end = colind_begin + s1;
            while (colind_begin != colind_end) {
                const size_t* const rowind_end = m_.index + *++colind_begin;
                while (rowind_begin != rowind_end) { ++(tmp_col[*rowind_begin++]); }
            }
            si[0] = 0;
            std::partial_sum(tmp_col, tmp_col + m_.s1, si.begin() + 1);
 
            rowind_begin = m_.index + m_.si[0];
            const T* value_begin = m_.m;
            for (size_t col = 0; col != s1; ++col) {
                const size_t* const rowind_end = m_.index + m_.si[col + 1];
                while (rowind_begin != rowind_end) {
                    size_t i = (si[*rowind_begin++])++;
                    index[i] = col;
                    m[i] = m_.scale * *value_begin++;
                }
            }
 
            si[0] = 0;
            std::partial_sum(tmp_col, tmp_col + m_.s1, si.begin() + 1);
            delete[] tmp_col;
        }
        return *this;
    }
 
    Matrix(const Matrix<T, MMProd<Mcompressed_col_st_constref, Mcompressed_col_st_constref>>& m_)
        : s1(m_.size1()), si(m_.size2() + 1) {
        *this = m_;
    }
 
    Matrix& operator=(
          const Matrix<T, MMProd<Mcompressed_col_st_constref, Mcompressed_col_st_constref>>& m_) {
        if (m_.trans || m_.m1.trans || m_.m2.trans) { UnimplementedError() << THROW; }
 
        s1 = m_.size1();
        si.resize(m_.size2() + 1);
        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());
        si[0] = 0;
 
        T scale = m_.scale * m_.m1.scale * m_.m2.scale;
 
        {
            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 = 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>(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());
 
                for (typename std::vector<std::pair<size_t, T>>::iterator i = res_tmp.begin();
                     i != res_tmp.end(); ++i) {
                    index.push_back(i->first);
                    m.push_back(i->second * scale);
                }
                si[++col] = m.size();
            }
        }
        return *this;
    }
 
    Matrix& operator=(
          const Matrix<
                T, Sum<Mcompressed_col_st_constref,
                       MMProd<Mcompressed_col_st_constref, Mcompressed_col_st_constref>>>& m_) {
        if (m_.trans || m_.m1.trans || m_.m2.m1.trans || m_.m2.m2.trans) {
            UnimplementedError() << THROW;
        }
 
        s1 = m_.size1();
        si.resize(m_.size2() + 1);
        size_t* tmp_index = TemporaryBuffer<size_t>::get(m_.m2.m1.size1());
        const size_t end_list_flag = m_.m2.m1.size1();
        T* tmp_value = TemporaryBuffer<T>::get(m_.m2.m1.size1());
        si[0] = 0;
 
        T scale_1 = m_.scale * m_.m1.scale;
        T scale_2 = m_.scale * m_.m2.scale * m_.m2.m1.scale * m_.m2.m2.scale;
 
        {
            const size_t* c_colind_begin = m_.m1.si;
            const size_t* b_colind_begin = m_.m2.m2.si;
            const size_t* const b_colind_end = b_colind_begin + m_.m2.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.m2.index[*b_colind_begin];
                const T* b_value_begin = &m_.m2.m2.m[*b_colind_begin];
                const size_t* b_rowind_end = &m_.m2.m2.index[*++b_colind_begin];
                while (b_rowind_begin != b_rowind_end) {
                    const size_t* a_rowind_begin = m_.m2.m1.index + m_.m2.m1.si[*b_rowind_begin];
                    const size_t* const a_rowind_end =
                          m_.m2.m1.index + m_.m2.m1.si[*b_rowind_begin + 1];
                    const T* a_value_begin = m_.m2.m1.m + m_.m2.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] += scale_2 * *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;
                    }
                }
 
                b_rowind_begin = &m_.m1.index[*c_colind_begin];
                b_value_begin = &m_.m1.m[*c_colind_begin];
                b_rowind_end = &m_.m1.index[*++c_colind_begin];
                while (b_rowind_begin != b_rowind_end) {
                    tmp_value[*b_rowind_begin] += scale_1 * *b_value_begin++;
                    if (!(std::numeric_limits<size_t>::max() - tmp_index[*b_rowind_begin])) {
                        tmp_index[*b_rowind_begin] = end_list;
                        end_list = *b_rowind_begin;
                    }
                    ++b_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, 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());
 
                for (typename std::vector<std::pair<size_t, T>>::iterator i = res_tmp.begin();
                     i != res_tmp.end(); ++i) {
                    index.push_back(i->first);
                    m.push_back(i->second);
                }
                si[++col] = m.size();
            }
        }
        return *this;
    }
 
    bool is_null() const { return this == &null; }
 
    void set_zero() {
        std::fill(si.begin(), si.end(), 0);
        index.clear();
        m.clear();
    }
 
    void clear() {
        s1 = 0;
        si.clear();
        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) {
        s1 = s1_;
        si.assign(s2_ + 1, 0);
        index.resize(snn);
        m.resize(snn);
    }
 
    void resize(const Index& idx) {
        s1 = idx.size();
        if (idx.sorted()) {
            for (std::vector<size_t>::iterator i = index.begin(); i != index.end(); ++i) {
                *i = idx[*i];
            }
        } else {
            for (size_t j = 0; j != size2(); ++j) {
                size_t* i_begin = &index[si[j]];
                size_t* i_end = &index[si[j + 1]];
                T* v_begin = &m[si[j]];
                std::vector<std::pair<size_t, T>> tmp;
                tmp.reserve(i_end - i_begin);
                while (i_begin != i_end) {
                    tmp.push_back(std::pair<size_t, T>(idx[*i_begin++], *v_begin++));
                }
                std::sort(tmp.begin(), tmp.end(), CompareFirstOfPair());
                i_begin = &index[si[j]];
                v_begin = &m[si[j]];
                for (typename std::vector<std::pair<size_t, T>>::iterator i = tmp.begin();
                     i != tmp.end(); ++i) {
                    *i_begin++ = i->first;
                    *v_begin++ = i->second;
                }
            }
        }
    }
 
    size_t get_nb_nonzero() const { return si.back(); }
 
    size_t get_nb_nonzero(size_t col) const { return si[col + 1] - si[col]; }
 
    void push_back(size_t row, size_t col, T value) {
        while (si.size() <= col + 1) { si.push_back(si.back()); }
        index.push_back(row);
        m.push_back(value);
        ++si.back();
        s1 = std::max(s1, row + 1);
    }
 
    void push_back(size_t row, size_t col, const Vector<T, Vdense_constref>& v) {
        while (si.size() <= col + 1) { si.push_back(si.back()); }
        m.insert(m.end(), v.v, v.v + v.s);
        for (size_t i = 0; i != v.s; ++i) { index.push_back(row + i); }
        si.back() = m.size();
        s1 = std::max(s1, row + v.s);
    }
 
    void push_back(size_t row, size_t col, const Index& ind, const Vector<T, Vdense_constref>& v) {
        while (si.size() <= col + 1) { si.push_back(si.back()); }
        for (size_t i = 0; i != v.size(); ++i) {
            m.push_back(v[i]);
            index.push_back(row + ind[i]);
        }
        si.back() = m.size();
        s1 = std::max(s1, row + v.s);
    }
 
    void push_back(size_t row, size_t col, const Vector<T, Vcompressed_scale_constref>& v) {
        while (si.size() <= col + 1) { si.push_back(si.back()); }
        for (size_t i = 0; i != v.snn; ++i) {
            m.push_back(v.scale * v.v[i]);
            index.push_back(row + v.index[i]);
        }
        si.back() = m.size();
        s1 = std::max(s1, row + v.s);
    }
 
    void push_back(size_t i, const Matrix<T, Mcompressed_col_st_constref>& m_) {
        if (m_.trans) { UnimplementedError() << THROW; }
        if (m_.size2() == 0) { return; }
        if (si.empty()) { si.push_back(0); }
        std::size_t ii = index.size();
        index.insert(index.end(), m_.index + m_.si[0], m_.index + m_.si[m_.s2]);
        for (std::vector<size_t>::iterator iii = index.begin() + ii; iii != index.end(); ++iii) {
            *iii += i;
        }
        size_t j = m.size();
        m.insert(m.end(), m_.m + m_.si[0], m_.m + m_.si[m_.s2]);
        if (m_.scale != T(1)) {
            for (; j != m.size(); ++j) { m[j] *= m_.scale; }
        }
 
        size_t is = si.size();
        si.insert(si.end(), m_.si + 1, m_.si + m_.s2 + 1);
        size_t iv = si[is - 1] - m_.si[0];
        for (std::vector<size_t>::iterator iii = si.begin() + is; iii != si.end(); ++iii) {
            *iii += iv;
        }
        s1 = std::max(s1, i + m_.size1());
    }
 
    void push_back(
          size_t i1, const Matrix<T, Mcompressed_col_constref>& m1, size_t i2,
          const Matrix<T, Mcompressed_col_constref>& m2) {
        if (i1 + m1.size1() > i2 && m1.size2() != m2.size2()) { UnimplementedError() << THROW; }
        if (si.empty()) { si.push_back(0); }
        for (size_t j = 0; j != m1.size2(); ++j) {
            size_t ii = index.size();
            index.insert(index.end(), m1.index + m1.si[j], m1.index + m1.si[j + 1]);
            for (std::vector<size_t>::iterator i = index.begin() + ii; i != index.end(); ++i) {
                *i += i1;
            }
            ii = index.size();
            index.insert(index.end(), m2.index + m2.si[j], m2.index + m2.si[j + 1]);
            for (std::vector<size_t>::iterator i = index.begin() + ii; i != index.end(); ++i) {
                *i += i2;
            }
            m.insert(m.end(), m1.m + m1.si[j], m1.m + m1.si[j + 1]);
            m.insert(m.end(), m2.m + m2.si[j], m2.m + m2.si[j + 1]);
            si.push_back(m.size());
        }
        s1 = std::max(s1, i2 + m2.size1());
    }
 
    void set_identity(size_t i1, size_t i2) {
        s1 = i1;
        si.resize(i2 + 1);
        {
            const size_t s = std::min(i1, i2);
            index.resize(s);
            m.resize(s);
        }
        size_t i = 0;
        for (; i != index.size(); ++i) {
            si[i] = i;
            index[i] = i;
            m[i] = 1;
        }
        std::fill(si.begin() + i, si.end(), i);
    }
 
    void push_back_identity(size_t i1, size_t i2) {
        if (si.empty()) {
            index.reserve(std::min(i1, i2));
            si.reserve(i2 + 1);
            m.reserve(index.size());
            si.push_back(0);
        }
        while (si.size() < i1 + 1) { si.push_back(si.back()); }
        while (i1 != i2) {
            index.push_back(i1);
            m.push_back(1);
            si.push_back(m.size());
            ++i1;
        }
        s1 = std::max(s1, i2);
    }
 
    void get_values(size_t*& colind, size_t*& rowind, T*& value) {
        colind = &si[0];
        rowind = &index[0];
        value = &m[0];
    }
 
    size_t get_col_value(const size_t col, size_t*& rowind, T*& value) {
        rowind = &index[si[col]];
        value = &m[si[col]];
        return si[col + 1] - si[col];
    }
 
    void set_value(
          size_t s1_, size_t s2_, size_t snn, const size_t* colind, const size_t* rowind,
          const T* value) {
        s1 = s1_;
        si.resize(s2_ + 1);
        std::copy(colind, colind + s2_ + 1, si.begin());
        index.resize(snn);
        std::copy(rowind, rowind + snn, index.begin());
        m.resize(snn);
        std::copy(value, value + snn, m.begin());
    }
 
    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; }
 
    T operator()(size_t i, size_t j) const {
        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;
    }
 
    Matrix<T, Msub_constref<Mcompressed_col_constref, Index, Index>> operator()(
          const Index& i, const Index& j) const {
        return Matrix<T, Msub_constref<Mcompressed_col_constref, Index, Index>>(
              Matrix<T, Mcompressed_col_constref>(*this), i, j);
    }
 
    Matrix<T, Mcompressed_col_constref> operator()(const Interval& j) const {
        return Matrix<T, Mcompressed_col_constref>(
              s1, j.size(), si[j.end] - si[j.start], &si[j.start], &index[0], &m[0]);
    }
 
    Vector<T, Vcompressed_scale_constref> operator[](size_t i) const {
        if (i >= size2()) { Exception() << THROW; }
        return Vector<T, Vcompressed_scale_constref>(
              s1, si[i + 1] - si[i], &index[si[i]], &m[si[i]], 1);
    }
 
    Matrix& operator=(const Matrix<T, Msub_constref<Mcompressed_col_constref, Index, Index>>& m_) {
        s1 = m_.size1();
        si.reserve(m_.size2() + 1);
        si.push_back(0);
        const size_t* jindex;
        const size_t* jindex_end;
        Index index_tmp;
        if (m_.index2.is_all()) {
            index_tmp.resize(m_.m.size1());
            for (size_t i = 0; i != index_tmp.size(); ++i) { index_tmp[i] = i; }
            jindex = &index_tmp[0];
            jindex_end = jindex + index_tmp.size();
        } else {
            jindex = &m_.index2[0];
            jindex_end = jindex + m_.index2.size();
        }
        if (m_.index1.is_all()) {
            for (; jindex != jindex_end; ++jindex) {
                size_t j = m_.m.si[*jindex];
                size_t j_end = m_.m.si[*jindex + 1];
                m.insert(m.end(), m_.m.m + j, m_.m.m + j_end);
                index.insert(index.end(), m_.m.index + j, m_.m.index + j_end);
                si.push_back(index.size());
            }
        } else {
            if (m_.index1.sorted()) {
                Index dual_iindex = m_.index1.make_dual();
                for (; jindex != jindex_end; ++jindex) {
                    size_t j = m_.m.si[*jindex];
                    size_t j_end = m_.m.si[*jindex + 1];
                    for (; j != j_end; ++j) {
                        if (dual_iindex[m_.m.index[j]] != dual_iindex.size()) {
                            index.push_back(m_.m.index[j]);
                            m.push_back(m_.m.m[j]);
                        }
                    }
                    si.push_back(index.size());
                }
            } else {
                Index dual_iindex = m_.index1.make_dual();
                for (; jindex != jindex_end; ++jindex) {
                    size_t j = m_.m.si[*jindex];
                    size_t j_end = m_.m.si[*jindex + 1];
                    std::vector<std::pair<size_t, T>> tmp;
                    tmp.reserve(j_end - j);
                    for (; j != j_end; ++j) {
                        tmp.push_back(std::pair<size_t, T>(dual_iindex[m_.m.index[j]], m_.m.m[j]));
                    }
                    std::sort(tmp.begin(), tmp.end(), CompareFirstOfPair());
                    for (typename std::vector<std::pair<size_t, T>>::iterator i = tmp.begin();
                         i != tmp.end(); ++i) {
                        index.push_back(i->first);
                        m.push_back(i->second);
                    }
                    si.push_back(index.size());
                }
            }
        }
        return *this;
    }
 
    void swap(Matrix& m_) {
        std::swap(s1, m_.s1);
        si.swap(m_.si);
        m.swap(m_.m);
        index.swap(m_.index);
    }
 
    void LUFactorization(
          Matrix<T, Mcompressed_col>& trans_L, Matrix<T, Mcompressed_col>& U, Index& P, Index& Q,
          Vector<double, Vdense>& R) {
        Matrix<T, Mcompressed_col_constref>(*this).LUFactorization(trans_L, U, P, Q, R);
    }
 
    size_t LUIFactorization(
          Matrix<T, Mcompressed_col>& L, Matrix<T, Mcompressed_col>& trans_U, Index& P, Index& Q,
          const bool compute_LU, Vector<double>& w = Vector<double>::null,
          const bool rook_pivot = false, double tol_pivot = -1, const double tol_drop = 3e-13,
          const double tol_rank_abs = 3.7e-11, const double tol_rank_rel = 3.7e-11,
          const Index& col_to_remove = Index::null);
 
    void LUIFactorization(
          Matrix<T, Mcompressed_col>& L, Matrix<T, Mcompressed_col>& trans_U, Index& P, Index& Q,
          const double tol_drop = 3e-13, const double tol_pivot_rank_search = 4.0,
          const double tol_pivot_LU_fact = 10.0, const double tol_rank_abs = 3.7e-11,
          const double tol_rank_rel = 3.7e-11) {
        // rank search
        Vector<double> w;
        const size_t rank = LUIFactorization(
              L, trans_U, Index::null, Index::null, false, w, true, tol_pivot_rank_search, tol_drop,
              tol_rank_abs, tol_rank_rel);
        Index col_to_remove;
        col_to_remove.reserve(w.size() - rank);
        for (size_t i = 0; i != w.size(); ++i) {
            if (w[i] <= 0) { col_to_remove.push_back(i); }
        }
 
        // LU factorisation
        LUIFactorization(
              L, trans_U, P, Q, true, Vector<double>::null, false, tol_pivot_LU_fact, tol_drop,
              tol_rank_abs, tol_rank_rel, col_to_remove);
    }
 
    void get_dep_columns(
          Index& Q, const double tol_drop = 3e-13, const double tol_pivot = 4.0,
          const double tol_rank_abs = 3.7e-11, const double tol_rank_rel = 3.7e-11) {
        Vector<double> w;
        Matrix<T, Mcompressed_col> L;
        Matrix<T, Mcompressed_col> trans_U;
        const size_t rank = LUIFactorization(
              L, trans_U, Index::null, Index::null, false, w, true, tol_pivot, tol_drop,
              tol_rank_abs, tol_rank_rel);
        Q.clear();
        Q.reserve(w.size() - rank);
        for (size_t i = 0; i != w.size(); ++i) {
            if (w[i] <= 0) { Q.push_back(i); }
        }
    }
 
    Matrix& operator*=(const double s) {
        for (size_t i = 0; i != m.size(); ++i) { m[i] *= s; }
        return *this;
    }
 
    Matrix& operator=(
          const Matrix<
                T,
                MMProd<Minverse_constref<Mcompressed_col_constref>, Mcompressed_col_st_constref>>&
                m_) {
        if (m_.scale != T(1) || m_.m2.scale != T(1) || m_.m2.trans) {
            UnimplementedError() << THROW;
        }
        if (m_.m1.size1() != m_.m1.size2()) { Exception() << THROW; }
 
        s1 = m_.m1.m.s1;
        si.clear();
        si.push_back(0);
        index.clear();
        m.clear();
        for (size_t k = 0; k != m_.m2.size2(); ++k) {
            if (m_.m2.si[k + 1] - m_.m2.si[k] > 0) {
                std::list<std::pair<size_t, T>> tmp;
                for (size_t j = m_.m2.si[k]; j != m_.m2.si[k + 1]; ++j) {
                    tmp.push_back(std::pair<size_t, T>(m_.m2.index[j], m_.m2.m[j]));
                }
                for (typename std::list<std::pair<size_t, T>>::reverse_iterator j = tmp.rbegin();
                     j != tmp.rend(); ++j) {
                    size_t jj_piv_p = m_.m1.m.si[j->first + 1] - 1;
                    if (m_.m1.m.index[jj_piv_p] != j->first) { Exception() << THROW; }
                    T piv = (j->second /= m_.m1.m.m[jj_piv_p]);
                    if (b2000::norm(piv) < 1e-15) { continue; }
                    typename std::list<std::pair<size_t, T>>::iterator jc = tmp.begin();
                    for (size_t jj = m_.m1.m.si[j->first]; jj < jj_piv_p; ++jj) {
                        size_t ii = m_.m1.m.index[jj];
                        T v = -m_.m1.m.m[jj] * piv;
                        while (jc != tmp.end() && jc->first < ii) { ++jc; }
                        if (jc != tmp.end() && jc->first == ii) {
                            jc->second += v;
                        } else {
                            tmp.insert(jc, std::pair<size_t, T>(ii, v));
                        }
                    }
                }
                for (typename std::list<std::pair<size_t, T>>::iterator i = tmp.begin();
                     i != tmp.end(); ++i) {
                    if (b2000::norm(i->second) > 1e-15) {
                        index.push_back(i->first);
                        m.push_back(i->second);
                    }
                }
            }
            si.push_back(index.size());
        }
        return *this;
    }
 
    Matrix& operator=(
          const Matrix<
                T,
                MMProd<Mcompressed_col_st_constref, Minverse_constref<Mcompressed_col_constref>>>&
                m_) {
        if (m_.scale != T(1) || m_.m1.scale != T(1) || m_.m1.trans) {
            UnimplementedError() << THROW;
        }
        if (m_.m2.size1() != m_.m2.size2()) { Exception() << THROW; }
 
        s1 = m_.m1.s1;
        si.clear();
        si.push_back(0);
        index.clear();
        m.clear();
        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());
        for (size_t k = 0; k != m_.m2.m.size1(); ++k) {
            size_t end_list = end_list_flag;
            const size_t* rowind_begin = m_.m1.index + m_.m1.si[k];
            const size_t* rowind_end = m_.m1.index + m_.m1.si[k + 1];
            const T* value_begin = m_.m1.m + m_.m1.si[k];
            while (rowind_begin != rowind_end) {
                tmp_value[*rowind_begin] = *value_begin++;
                tmp_index[*rowind_begin] = end_list;
                end_list = *rowind_begin++;
            }
            rowind_begin = m_.m2.m.index + m_.m2.m.si[k];
            value_begin = m_.m2.m.m + m_.m2.m.si[k];
            while (*rowind_begin < k) {
                const size_t* a_rowind_begin = &index[si[*rowind_begin]];
                const T* a_value_begin = &m[si[*rowind_begin]];
                const size_t* a_rowind_end = &index[si[*rowind_begin + 1]];
                ++rowind_begin;
                while (a_rowind_begin != a_rowind_end) {
                    tmp_value[*a_rowind_begin] -= *a_value_begin++ * *value_begin;
                    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;
                }
                ++value_begin;
            }
            T pivot = T(1) / *value_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, pivot * 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());
 
            for (typename std::vector<std::pair<size_t, T>>::iterator i = res_tmp.begin();
                 i != res_tmp.end(); ++i) {
                index.push_back(i->first);
                m.push_back(i->second);
            }
            si.push_back(m.size());
        }
 
        return *this;
    }
 
    template <typename T1, typename STORAGE>
    void scale_row(const Vector<T1, STORAGE>& v) {
        for (size_t i = 0; i != m.size(); ++i) { m[i] *= v[index[i]]; }
    }
 
    template <typename T1, typename STORAGE>
    void scale_col(const Vector<T1, STORAGE>& v) {
        if (m.empty()) { return; }
        for (size_t j = 0; j != size2(); ++j) {
            T* i = &m[si[j]];
            T* i_end = &m[si[j + 1]];
            T1 tmp = v[j];
            while (i != i_end) { *i++ *= tmp; }
        }
    }
 
    void scale_col_to_norminf(Vector<double, Vdense>& v) {
        if (m.empty()) { return; }
        for (size_t j = 0; j != size2(); ++j) {
            T* i_end = &m[si[j + 1]];
            T tmp = 0;
            for (T* i = &m[si[j]]; i != i_end; ++i) {
                tmp = max_abs(tmp, *i);  // std::max(tmp, b2000::norm(*i));
            }
            if (tmp == T(0)) { continue; }
            tmp = T(1) / tmp;
            v[j] = b2000::real(tmp);
            for (T* i = &m[si[j]]; i != i_end; ++i) { *i *= tmp; }
        }
    }
 
    template <typename T1>
    void scale_invert_col(const Vector<T1, Vdense_constref>& v) {
        if (m.empty()) { return; }
        T* ii = &m[0];
        for (size_t i = 0; i != v.size(); ++i) {
            T1 tmp = T1(1) / v[i];
            for (T const* i_end = &m[si[i + 1]]; ii != i_end; ++ii) { *ii *= tmp; }
        }
    }
 
    bool is_null_value() const {
        for (size_t i = 0; i != m.size(); ++i) {
            if (m[i] != T(0)) { return false; }
        }
        return true;
    }
 
    void remove_zero(const double tol = 0) {
        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);
    }
 
    void row_permute(size_t new_s1, const Index perm) {
        if (perm.size() != s1) { Exception() << THROW; }
        for (size_t i = 0; i != index.size(); ++i) { index[i] = perm[index[i]]; }
        s1 = new_s1;
    }
 
    void set_full() {
        {
            std::vector<T> m_tmp(s1 * size2());
            size_t ptri = 0;
            size_t ptro = 0;
            for (size_t j = 1; j != si.size(); ++j, ptro += s1) {
                for (; ptri != si[j]; ++ptri) { m_tmp[ptro + index[ptri]] = m[ptri]; }
            }
            m.swap(m_tmp);
        }
        index.resize(m.size());
        si[0] = 0;
        size_t ptr = 0;
        for (size_t j = 1; j != si.size(); ++j) {
            si[j] = si[j - 1] + s1;
            for (size_t i = 0; i != s1; ++i, ++ptr) { index[ptr] = i; }
        }
    }
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        l << "column compressed matrix of size (" << m.size1() << ", " << m.size2() << ") ";
        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;
    }
 
    friend std::ostream& operator<<(std::ostream& out, const Matrix& m) {
        /*
        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)
                out << "(" << m.index[i] << ", " << j << ") = " << m.m[i] << std::endl;
 
        */
        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;
    }
 
    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, Mcompressed_col> Matrix<T, Mcompressed_col>::null;
 
// Right looking LU factorisation
 
// The output L is triangular inferior unit matrix. The output
// trans_L is triangular inferior matrix without zero on the diagonal.
// At output, L.size2() = trans_U.size2() = rank of the input  matrix.
template <typename T>
size_t Matrix<T, Mcompressed_col>::LUIFactorization(
      Matrix<T, Mcompressed_col>& L, Matrix<T, Mcompressed_col>& trans_U, Index& P, Index& Q,
      const bool compute_LU, Vector<double>& w, const bool rook_pivot, double tol_pivot,
      const double drop_tol, const double tol_rank_abs, const double tol_rank_rel,
      const Index& col_to_remove) {
    if (tol_pivot < 0) { tol_pivot = rook_pivot ? 4.0 : 10.0; }
 
    const size_t s2 = size2();
    const size_t min_mn = std::min(s1, s2);
 
    std::set<size_t> empty_row;
    std::set<size_t> empty_col;
    size_t* P_empty_row = 0;
    size_t* Q_empty_col = 0;
    if (compute_LU) {
        P.resize(s1);
        P_empty_row = &P.back();
        Q.resize(s2);
        Q_empty_col = &Q.back();
    }
 
    L.s1 = s1;
    L.si.clear();
    L.si.push_back(0);
    L.index.clear();
    L.m.clear();
 
    trans_U.s1 = s2;
    trans_U.si.clear();
    trans_U.si.push_back(0);
    trans_U.index.clear();
    trans_U.m.clear();
 
    Vector<double> U_diag;
    if (!w.is_null()) {
        U_diag.resize(s2);
        U_diag.set_zero();
        w.resize(s2);
        w.set_zero();
    }
 
    std::vector<std::map<size_t, T>> col(s2);
    std::vector<std::map<size_t, T*>> row(s1);
    std::vector<double> col_max(s2);
    std::vector<double> row_max(s1);
    {
        const size_t* col_to_remove_ptr = 0;
        const size_t* col_to_remove_ptr_end = 0;
        if (!col_to_remove.is_null() && !col_to_remove.empty()) {
            col_to_remove_ptr = &col_to_remove[0];
            col_to_remove_ptr_end = &col_to_remove.back() + 1;
        }
        size_t i = si[0];
        for (size_t j = 0; j != s2; ++j) {
            if (col_to_remove_ptr && col_to_remove_ptr != col_to_remove_ptr_end
                && j == *col_to_remove_ptr) {
                ++col_to_remove_ptr;
                i = si[j + 1];
                continue;
            }
            const size_t i_end = si[j + 1];
            for (; i != i_end; ++i) {
                const size_t ii = index[i];
                const double nm = b2000::norm(m[i]);
                if (nm > drop_tol) {
                    T* ptr;
                    {
                        typename std::map<size_t, T>& col_j = col[j];
                        ptr = &(col_j.insert(col_j.end(), std::pair<size_t, T>(ii, m[i]))->second);
                    }
                    {
                        typename std::map<size_t, T*>& row_i = row[ii];
                        row_i.insert(row_i.end(), std::pair<size_t, T*>(j, ptr));
                    }
                    if (nm > col_max[j]) { col_max[j] = nm; }
                    if (nm > row_max[ii]) { row_max[ii] = nm; }
                }
            }
        }
    }
 
    // initialise the row/com degree map
    using degree_t = std::multimap<size_t, size_t>;
    using degree_iter_t = std::vector<degree_t::iterator>;
 
    degree_t col_degree;
    degree_iter_t col_degree_iter(s2);
    degree_t row_degree;
    degree_iter_t row_degree_iter(s1);
    for (size_t i = 0; i != s2; ++i) {
        if (!col[i].empty()) {
            col_degree_iter[i] = col_degree.insert(std::pair<size_t, size_t>(col[i].size(), i));
        } else {
            if (Q_empty_col) { *Q_empty_col-- = i; }
            empty_col.insert(i);
        }
    }
    for (size_t i = 0; i != s1; ++i) {
        if (!row[i].empty()) {
            row_degree_iter[i] = row_degree.insert(std::pair<size_t, size_t>(row[i].size(), i));
        } else {
            if (P_empty_row) { *P_empty_row-- = i; }
            empty_row.insert(i);
        }
    }
 
    // iteration on the LU decomposition
    size_t k;
    for (k = 0; k != min_mn; ++k) {
        size_t i_pivot = s1;
        size_t j_pivot = s2;
 
        // chose a pivot using the Markowitz strategy to reduce the
        // fill-in and the threshold partial pivoting or the threshold
        // rook pivoting strategy to stabilise the LU factorisation.
        {
            degree_t::const_iterator j = col_degree.begin();
            degree_t::const_iterator i = row_degree.begin();
            if (j == col_degree.end() || i == row_degree.end()) { break; }
            assert(col[j->second].size() == j->first);
            assert(row[i->second].size() == i->first);
            double min_tol = std::numeric_limits<double>::max();
            double min_degree = s1 * s2;
            const size_t min_degree_all = i->first * j->first;
            for (;;) {
                const size_t min_degree_all1 = i->first * j->first;
                if (j->first <= i->first) {
                    typename std::map<size_t, T>::const_iterator i1 = col[j->second].begin();
                    if (j->first == 1 && row_degree_iter[i1->first] != row_degree.end()
                        && col_degree_iter[j->second] != col_degree.end()) {
                        // The order we take it does not influence the
                        // fill-in nor the numerical stability of the
                        // factorization.
                        i_pivot = i1->first;
                        j_pivot = j->second;
                        assert(col[j_pivot].find(i_pivot) != col[j_pivot].end());
                        break;
                    }
                    const typename std::map<size_t, T>::const_iterator i1_end =
                          col[j->second].end();
                    for (; i1 != i1_end; ++i1) {
                        const size_t d = j->first * row[i1->first].size();
                        if (d <= min_degree) {
                            const T iv = T(1) / i1->second;
                            double vc = b2000::norm(col_max[j->second] * iv);
                            if (rook_pivot) {
                                vc = std::max(vc, double(b2000::norm(row_max[i1->first] * iv)));
                            }
                            if (vc < min_tol && vc < tol_pivot
                                && row_degree_iter[i1->first] != row_degree.end()
                                && col_degree_iter[j->second] != col_degree.end()) {
                                min_tol = vc;
                                min_degree = d;
                                i_pivot = i1->first;
                                j_pivot = j->second;
                                assert(col[j_pivot].find(i_pivot) != col[j_pivot].end());
                            }
                        }
                    }
                    ++j;
                } else {
                    typename std::map<size_t, T*>::const_iterator j1 = row[i->second].begin();
                    if (i->first == 1 && row_degree_iter[i->second] != row_degree.end()
                        && col_degree_iter[j1->first] != col_degree.end()) {
                        // The order we take it does not influence the
                        // fill-in nor the numerical stability of the
                        // factorization.
                        i_pivot = i->second;
                        j_pivot = j1->first;
                        assert(col[j_pivot].find(i_pivot) != col[j_pivot].end());
                        break;
                    }
                    const typename std::map<size_t, T*>::const_iterator j1_end =
                          row[i->second].end();
                    for (; j1 != j1_end; ++j1) {
                        const size_t d = i->first * col[j1->first].size();
                        if (d <= min_degree) {
                            const T iv = T(1) / *(j1->second);
                            double vc = b2000::norm(col_max[j1->first] * iv);
                            if (rook_pivot) {
                                vc = std::min(vc, double(b2000::norm(row_max[i->second] * iv)));
                            }
                            if (vc < min_tol && vc < tol_pivot
                                && row_degree_iter[i->second] != row_degree.end()
                                && col_degree_iter[j1->first] != col_degree.end()) {
                                min_tol = vc;
                                min_degree = d;
                                i_pivot = i->second;
                                j_pivot = j1->first;
                                assert(col[j_pivot].find(i_pivot) != col[j_pivot].end());
                            }
                        }
                    }
                    ++i;
                }
                if (min_tol < 1.0001 && min_degree == min_degree_all) { break; }
                if (min_tol < tol_pivot && min_degree_all1 > min_degree_all) { break; }
                if (j == col_degree.end() || i == row_degree.end()) { break; }
            }
        }
 
        if (i_pivot == s1) { Exception() << THROW; }
 
        col_degree.erase(col_degree_iter[j_pivot]);
        col_degree_iter[j_pivot] = col_degree.end();
        row_degree.erase(row_degree_iter[i_pivot]);
        row_degree_iter[i_pivot] = row_degree.end();
        if (compute_LU) {
            P[k] = i_pivot;
            Q[k] = j_pivot;
        }
 
        typename std::map<size_t, T>::iterator i = col[j_pivot].find(i_pivot);
        assert(i != col[j_pivot].end());
        T value_pivot = i->second;
 
        // computation of the pivot column
        {
            const T inv_value_pivot = T(1) / value_pivot;
            i = col[j_pivot].begin();
            const typename std::map<size_t, T>::const_iterator i_end = col[j_pivot].end();
            for (; i != i_end; ++i) {
                L.index.push_back(i->first);
                L.m.push_back(i->first == i_pivot ? T(1) : i->second * inv_value_pivot);
                row[i->first].erase(j_pivot);
            }
        }
        col[j_pivot].clear();
 
        // update the sub matrix
        trans_U.index.push_back(j_pivot);
        trans_U.m.push_back(value_pivot);
        {
            typename std::map<size_t, T*>::const_iterator j = row[i_pivot].begin();
            const typename std::map<size_t, T*>::const_iterator j_end = row[i_pivot].end();
            for (; j != j_end; ++j) {
                typename std::map<size_t, T>& col_j = col[j->first];
                typename std::map<size_t, T>::iterator ii = col_j.find(i_pivot);
                assert(ii != col_j.end());
                trans_U.index.push_back(j->first);
                trans_U.m.push_back(ii->second);
                col_j.erase(ii);
                for (size_t iii = L.si.back(); iii != L.index.size(); ++iii) {
                    if (L.index[iii] != i_pivot) {
                        const T v = -L.m[iii] * trans_U.m.back();
                        ii = col_j.find(L.index[iii]);
                        if (ii == col_j.end()) {
                            T* ptr = &(col_j.insert(std::pair<size_t, T>(L.index[iii], v))
                                             .first->second);
                            row[L.index[iii]].insert(std::pair<size_t, T*>(j->first, ptr));
                        } else {
                            ii->second += v;
                            if (b2000::norm(ii->second) <= drop_tol) {
                                col_j.erase(ii);
                                row[L.index[iii]].erase(j->first);
                            }
                        }
                    }
                }
            }
        }
        row[i_pivot].clear();
 
        // Update col degree map and maximum values.
        for (size_t jj = trans_U.si.back(); jj != trans_U.index.size(); ++jj) {
            const size_t j = trans_U.index[jj];
            if (j == j_pivot) { continue; }
            {
                typename std::map<size_t, T>::const_iterator i = col[j].begin();
                const typename std::map<size_t, T>::const_iterator i_end = col[j].end();
                double max = 0;
                for (; i != i_end; ++i) { max = std::max(max, double(b2000::norm(i->second))); }
                col_max[j] = max;
            }
            if (col_degree_iter[j] != col_degree.end()) { col_degree.erase(col_degree_iter[j]); }
            col_degree_iter[j] = col_degree.end();
            if (!col[j].empty()) {
                col_degree_iter[j] = col_degree.insert(std::pair<size_t, size_t>(col[j].size(), j));
            } else {
                if (Q_empty_col && empty_col.find(j) == empty_col.end()) { *Q_empty_col-- = j; }
                empty_col.insert(j);
            }
        }
 
        // Update row degree map and maximum values.
        for (size_t ii = L.si.back(); ii != L.index.size(); ++ii) {
            const size_t i = L.index[ii];
            if (i == i_pivot) { continue; }
            {
                typename std::map<size_t, T*>::const_iterator j = row[i].begin();
                const typename std::map<size_t, T*>::const_iterator j_end = row[i].end();
                double max = 0;
                for (; j != j_end; ++j) { max = std::max(max, double(b2000::norm(*(j->second)))); }
                row_max[i] = max;
            }
            if (row_degree_iter[i] != row_degree.end()) { row_degree.erase(row_degree_iter[i]); }
            row_degree_iter[i] = row_degree.end();
            if (!row[i].empty()) {
                row_degree_iter[i] = row_degree.insert(std::pair<size_t, size_t>(row[i].size(), i));
            } else {
                if (P_empty_row && empty_row.find(i) == empty_row.end()) { *P_empty_row-- = i; }
                empty_row.insert(i);
            }
        }
 
        // update w
        if (!w.is_null()) {
            U_diag[j_pivot] = b2000::norm(value_pivot);
            for (size_t jj = trans_U.si.back(); jj != trans_U.index.size(); ++jj) {
                const size_t j = trans_U.index[jj];
                w[j] = std::max(w[j], double(b2000::norm(trans_U.m[jj])));
            }
        }
 
        if (compute_LU) {
            L.si.push_back(L.index.size());
            trans_U.si.push_back(trans_U.index.size());
        } else {
            L.index.clear();
            L.m.clear();
            trans_U.index.clear();
            trans_U.m.clear();
        }
    }
 
    if (compute_LU) {
        {
            size_t kk = k;
            for (degree_t::const_iterator i = row_degree.begin(); i != row_degree.end();
                 ++i, ++kk) {
                P[kk] = i->second;
            }
            std::sort(P.begin() + k, P.end());
        }
 
        using p_iter =
              pair_iterator<std::vector<size_t>::iterator, typename std::vector<T>::iterator>;
        {
            Index inv_P = P.make_dual();
            size_t i = L.si[0];
            for (Index::const_iterator j = L.si.begin() + 1; j != L.si.end(); ++j) {
                const size_t i_end = *j;
                p_iter ii_begin(L.index.begin() + i, L.m.begin() + i);
                for (; i != i_end; ++i) {
                    size_t& ii = L.index[i];
                    ii = inv_P[ii];
                }
                p_iter ii_end(L.index.begin() + i, L.m.begin() + i);
                std::sort(ii_begin, ii_end, CompareFirstOfPair());
            }
        }
 
        {
            size_t kk = k;
            for (degree_t::const_iterator j = col_degree.begin(); j != col_degree.end();
                 ++j, ++kk) {
                Q[kk] = j->second;
            }
            std::sort(Q.begin() + k, Q.end());
        }
        {
            Index inv_Q = Q.make_dual();
            size_t i = trans_U.si[0];
            for (Index::const_iterator j = trans_U.si.begin() + 1; j != trans_U.si.end(); ++j) {
                const size_t i_end = *j;
                p_iter ii_begin(trans_U.index.begin() + i, trans_U.m.begin() + i);
                for (; i != i_end; ++i) {
                    size_t& ii = trans_U.index[i];
                    ii = inv_Q[ii];
                }
                p_iter ii_end(trans_U.index.begin() + i, trans_U.m.begin() + i);
                std::sort(ii_begin, ii_end, CompareFirstOfPair());
            }
        }
    }
 
    // set w[j] <= 0 if dependent column detected
    size_t rank = s2;
    if (!w.is_null()) {
        for (size_t j = 0; j != s2; ++j) {
            if (U_diag[j] < tol_rank_abs || U_diag[j] < tol_rank_rel * w[j]) {
                w[j] *= -1;
                --rank;
            }
        }
    }
    return rank;
}
 
struct Mcompressed_col_ref {
    using base = Mcompressed_col_ref;
    using const_base = Mcompressed_col_st_constref;
    using copy = Mcompressed_col;
};
 
template <typename T>
class Matrix<T, Mcompressed_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, Mcompressed_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; }
 
    bool is_null_value() const {
        for (size_t i = si[0]; i != si[s2]; ++i) {
            if (m[i] != 0) { return false; }
        }
        return true;
    }
 
    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 {
        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;
    }
 
    template <typename T1>
    void scale_row(const Vector<T1, Vdense_constref>& v) {
        for (size_t i = si[0]; i != si[s2]; ++i) { m[i] *= v[index[i]]; }
    }
 
    template <typename T1>
    void scale_col(const Vector<T1, Vdense_constref>& v) {
        if (si[s2] == 0) { return; }
        for (size_t j = 0; j != size2(); ++j) {
            T* i = m + si[j];
            T* i_end = m + si[j + 1];
            T1 tmp = v[j];
            while (i != i_end) { *i *= tmp; }
        }
    }
 
    template <typename T1>
    void scale_invert_col(const Vector<T1, Vdense_constref>& v) {
        if (si[s2] == 0) { return; }
        T* ii = m;
        for (size_t i = 0; i != v.size(); ++i) {
            T1 tmp = T1(1) / v[i];
            for (T const* i_end = m + index[i + 1]; ii != i_end; ++ii) { *ii *= tmp; }
        }
    }
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        l << "column compressed matrix of size (" << m.size1() << ", " << m.size2() << ") ";
        l.write(m.s2 + 1, m.si, 1, "colind");
        l.write(m.si[m.s2], m.index + m.si[0], 1, "rowind");
        l.write(m.si[m.s2], m.m + m.si[0], 1, "value");
        return l;
    }
 
    friend std::ostream& operator<<(std::ostream& out, const Matrix& m) {
        size_t i = m.si[0];
        for (size_t j = 0; j != m.s2; ++j) {
            for (size_t i_end = m.si[j + 1]; i != i_end; ++i) {
                out << "(" << m.index[i] << ", " << j << ") = " << m.m[i] << std::endl;
            }
        }
        return out;
    }
 
    static Matrix null;
 
private:
    size_t s1;
    size_t s2;
    size_t* si;
    size_t* index;
    T* m;
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mcompressed_col_ref> Matrix<T, Mcompressed_col_ref>::null;
 
struct Mcompressed_col_constref {
    using base = Mcompressed_col_st_constref;
    using const_base = Mcompressed_col_st_constref;
    using copy = Mcompressed_col;
};
 
template <typename T>
class Matrix<T, Mcompressed_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, Mcompressed_col_ref>& m_)
        : s1(m_.s1), s2(m_.s2), si(m_.si), index(m_.index), m(m_.m) {}
 
    Matrix(const Matrix<T, Mcompressed_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; }
 
    bool is_null_value() const {
        for (size_t i = si[0]; i != si[s2]; ++i) {
            if (m[i] != T(0)) { return false; }
        }
        return true;
    }
 
    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 {
        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 << "column compressed matrix of size (" << m.size1() << ", " << m.size2() << ") ";
        l.write(m.s2 + 1, m.si, 1, "colind");
        l.write(m.si[m.s2], m.index + m.si[0], 1, "rowind");
        l.write(m.si[m.s2], m.m + m.si[0], 1, "value");
        return l;
    }
 
    friend std::ostream& operator<<(std::ostream& out, const Matrix& m) {
        size_t i = m.si[0];
        for (size_t j = 0; j != m.s2; ++j) {
            for (size_t i_end = m.si[j + 1]; i != i_end; ++i) {
                out << "(" << m.index[i] << ", " << j << ") = " << m.m[i] << std::endl;
            }
        }
        return out;
    }
 
    void LUFactorization(
          Matrix<T, Mcompressed_col>& trans_L, Matrix<T, Mcompressed_col>& U, Index& P, Index& Q,
          Vector<double, Vdense>& R);
 
    static Matrix null;
 
private:
    static void UMFPACK_clean_incomplete_factorization(
          Matrix<T, Mcompressed_col>& trans_L, Matrix<T, Mcompressed_col>& U, const size_t nb_udiag,
          Index& P, Index& Q, Vector<double, Vdense>& R);
 
    size_t s1;
    size_t s2;
    const size_t* si;
    const size_t* index;
    const T* m;
    MVFRIEND;
};
 
template <>
void Matrix<double, Mcompressed_col_constref>::LUFactorization(
      Matrix<double, Mcompressed_col>& trans_L, Matrix<double, Mcompressed_col>& U, Index& P,
      Index& Q, Vector<double, Vdense>& R);
 
template <>
void Matrix<b2000::csda<double>, Mcompressed_col_constref>::LUFactorization(
      Matrix<b2000::csda<double>, Mcompressed_col>& trans_L,
      Matrix<b2000::csda<double>, Mcompressed_col>& U, Index& P, Index& Q,
      Vector<double, Vdense>& R);
 
template <>
void Matrix<std::complex<double>, Mcompressed_col_constref>::LUFactorization(
      Matrix<std::complex<double>, Mcompressed_col>& trans_L,
      Matrix<std::complex<double>, Mcompressed_col>& U, Index& P, Index& Q,
      Vector<double, Vdense>& R);
 
template <typename T>
Matrix<T, Mcompressed_col_constref> Matrix<T, Mcompressed_col_constref>::null;
 
struct Mcompressed_col_st_constref {
    using base = Mcompressed_col_st_constref;
    using const_base = Mcompressed_col_st_constref;
    using copy = Mcompressed_col;
};
 
template <typename T>
class Matrix<T, Mcompressed_col_st_constref> {
public:
    Matrix() : s1(0), s2(0), si(0), index(0), m(0), scale(0), trans(0) {}
 
    Matrix(
          size_t s1_, size_t s2_, size_t snn_, const size_t* si_, const size_t* index_, const T* m_,
          T scale_, bool trans_)
        : s1(s1_), s2(s2_), si(si_), index(index_), m(m_), scale(scale_), trans(trans_) {}
 
    Matrix(const Matrix& m_)
        : s1(m_.s1),
          s2(m_.s2),
          si(m_.si),
          index(m_.index),
          m(m_.m),
          scale(m_.scale),
          trans(m_.trans) {}
 
    Matrix(const Matrix<T, Mcompressed_col_ref>& m_)
        : s1(m_.s1), s2(m_.s2), si(m_.si), index(m_.index), m(m_.m), scale(1), trans(false) {}
 
    Matrix(const Matrix<T, Mcompressed_col_constref>& m_)
        : s1(m_.s1), s2(m_.s2), si(m_.si), index(m_.index), m(m_.m), scale(1), trans(false) {}
 
    Matrix(const Matrix<T, Mcompressed_col>& m_)
        : s1(m_.s1),
          s2(m_.si.size() - 1),
          si(&m_.si[0]),
          index(&m_.index[0]),
          m(&m_.m[0]),
          scale(1),
          trans(false) {}
 
    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;
        }
    }
 
    T operator()(size_t i, size_t j) const {
        if (trans) { 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;
    }
 
    Matrix& transpose() {
        trans = trans ? false : true;
        return *this;
    }
 
    static Matrix null;
 
private:
    bool operator==(const Matrix& m_) const {
        return s1 == m_.s1 && s2 == m_.s2 && si == m_.si && index == m_.index && m == m_.m
               && scale == m_.scale && trans == m_.trans;
    }
 
    size_t s1;
    size_t s2;
    const size_t* si;
    const size_t* index;
    const T* m;
    T scale;
    char trans;
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mcompressed_col_st_constref> Matrix<T, Mcompressed_col_st_constref>::null;
 
struct Mcompressed_col_update_inv {
    using base = Mcompressed_col_ref;
    using const_base = Mcompressed_col_st_constref;
    using copy = Mcompressed_col_update_inv;
    using inverse = Mcompressed_col_update_inv;
};
 
template <typename T>
class Matrix<T, Mcompressed_col_update_inv> {
public:
    Matrix(
          size_t s1_ = 0,
          SparseMatrixConnectivityType connectivity_ = sparse_matrix_connectivity_unknown,
          const Dictionary& dictionary_ = Dictionary::get_empty())
        : s1(s1_), value(s1_), solver(0), connectivity(connectivity_), dictionary(&dictionary_) {}
 
    Matrix(const Matrix& m_)
        : s1(m_.s1),
          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);
            s1 = 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()) {
        s1 = s1_;
        if (si.empty()) {
            value.resize(s2);
        } else {
            UnimplementedError() << THROW;
        }
        if (connectivity_ != sparse_matrix_connectivity_unknown) { connectivity = connectivity_; }
        if (&dictionary_ != &Dictionary::get_empty()) { dictionary = &dictionary_; }
    }
 
    void set_same_structure(const Matrix& m_) {
        s1 = m_.s1;
        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() {
        std::fill(si.begin(), si.end(), 0);
        index.resize(0);
        m.resize(0);
        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>(s1, s);
    }
 
    size_t size1() const { return s1; }
 
    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));
        }
    }
 
    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();
    }
 
    T operator()(size_t i, size_t j) const {
        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 (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, 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()]; }
        }
        static T zero = 0;
        return zero;
    }
 
    Matrix& operator+=(
          const Matrix<T, Mstructured_constref<Mrectangle_increment_st_constref>>& m_) {
        if (size() != m_.size()) {
            Exception() << "The two matrix do not have the same size, " << size() << " and "
                        << m_.size() << THROW;
        }
        for (size_t j = 0; j != m_.m.s2; ++j) {
            std::vector<std::pair<size_t, T>> tmp;
            tmp.reserve(m_.m.s1);
            for (size_t i = 0; i != m_.m.s1; ++i) {
                tmp.push_back(std::pair<size_t, T>(m_.index[i], m_.m(i, j)));
            }
            std::sort(tmp.begin(), tmp.end(), CompareFirstOfPair());
            if (m_.index2) {
                add_colomn(m_.index2[j], tmp.begin(), tmp.end());
            } else {
                add_colomn(m_.index[j], tmp.begin(), tmp.end());
            }
        }
        return *this;
    }
 
    Matrix& operator+=(const Matrix<
                       T, MMProd<
                                MMProd<
                                      Mcompressed_col_st_constref,
                                      Mstructured_constref<Mrectangle_increment_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_.m1.m1.trans == true || m_.m1.m2.m.s1 != m_.m1.m2.m.s2 || 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.s1;
        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, Mrectangle_increment_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 = 0; 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(), output_col.end());
            }
        } else {
            for (size_t i = 0; i != output_size; ++i) {
                tmp_rowind[output_dof_numbering[i].first] = 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(), 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, Mcompressed_col_st_constref>,
                                Mcompressed_col_st_constref>>& m_) {
        UnimplementedError() << THROW;
        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) { UnimplementedError() << THROW; }
 
        Matrix<T, Mcompressed_col> m_m2_tmp;
        if (m_.m2.trans) { m_m2_tmp = m_.m2; }
 
        Matrix<T, Mcompressed_col_st_constref> m_m2(
              m_.m2.trans ? Matrix<T, Mcompressed_col_st_constref>(m_m2_tmp) : 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 * m_m2.scale;
        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 = 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>(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();
        blas::axpy(m.size(), -1, &m_.m[0], 1, &m[0], 1);
        return *this;
    }
 
    Matrix& operator+=(const Matrix& m_) {
        value_to_ccarray();
        const_cast<Matrix&>(m_).value_to_ccarray();
        blas::axpy(m.size(), 1, &m_.m[0], 1, &m[0], 1);
        return *this;
    }
 
    Matrix& operator+=(const Matrix<T, Mcompressed_col_st_constref>& m_) {
        if (m_.size1() != m_.size2()) { UnimplementedError() << THROW; }
        value_to_ccarray();
        if (m_.trans) { UnimplementedError() << THROW; }
        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<Mcompressed_col_st_constref, Mcompressed_col_st_constref>>& m_) {
        if (m_.size1() != m_.size2()) { UnimplementedError() << THROW; }
        if (m_.m1.trans || m_.m2.trans || !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;
        }
        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]]; }
        }
    }
 
    void remove_empty_column(Index& index) {
        if (!si.empty()) { UnimplementedError() << THROW; }
        std::vector<std::vector<std::pair<size_t, T>>> value_tmp;
        index.clear();
        for (size_t i = 0; i != value.size(); ++i) {
            if (!value[i].empty()) {
                value_tmp.push_back(std::vector<std::pair<size_t, T>>());
                value_tmp.back().swap(value[i]);
                index.push_back(i);
            }
        }
        value.swap(value_tmp);
        value_to_ccarray();
    }
 
    void remove_empty_column(Index& index, const double tol) {
        if (!si.empty()) { UnimplementedError() << THROW; }
        std::vector<std::vector<std::pair<size_t, T>>> value_tmp;
        index.clear();
        for (size_t i = 0; i != value.size(); ++i) {
            for (typename std::vector<std::pair<size_t, T>>::const_iterator j = value[i].begin();
                 j != value[i].end(); ++j) {
                if (b2000::norm(j->second) > tol) {
                    value_tmp.push_back(std::vector<std::pair<size_t, T>>());
                    value_tmp.back().swap(value[i]);
                    index.push_back(i);
                    break;
                }
            }
        }
        value.swap(value_tmp);
        value_to_ccarray();
    }
 
    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);
    }
 
    Vector<T, Vindex1_constref> get_diagonal() {
        set_diag_index();
        return Vector<T, Vindex1_constref>(
              diag_index.size(), &m[0], diag_index.back(), &diag_index[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, Mcompressed_col_update_sub_ref> operator()(const Interval& i, const Interval& j) {
        return Matrix<T, Mcompressed_col_update_sub_ref>(*this, i, j);
    }
 
    operator const Matrix<T, Mcompressed_col_st_constref>() const {
        const_cast<Matrix*>(this)->value_to_ccarray();
        return Matrix<T, Mcompressed_col_st_constref>(
              s1, si.size() - 1, m.size(), &si[0], &index[0], &m[0], 1, false);
    }
 
    friend logging::Logger& operator<<(logging::Logger& l, const Matrix& m) {
        const_cast<Matrix&>(m).value_to_ccarray();
        l << "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;
    }
 
    size_t get_null_space_size() {
        if (si.size() <= 1) { return 0; }
        if (solver == 0) { return 0; }
        Matrix* noconst_this = const_cast<Matrix<T, Mcompressed_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; }
        Matrix* noconst_this = const_cast<Matrix<T, Mcompressed_col_update_inv>*>(this);
        noconst_this->value_to_ccarray();
        if (solver == 0) {
            noconst_this->solver = LU_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->get_null_space(nks.size1(), nks.size2(), nks.m, nks.size1());
    }
 
    static Matrix null;
 
private:
    void set_diag_index() {
        value_to_ccarray();
        if (diag_index.empty()) {
            diag_index.resize(si.size() - 1);
            size_t i = si[0];
            for (size_t j = 0; j != diag_index.size(); ++j) {
                for (const size_t i_end = si[j + 1]; i != i_end; ++i) {
                    if (index[i] == j) { diag_index[j] = i; }
                }
            }
        }
    }
 
    bool value_to_ccarray() {
        if (value.empty()) {
            if (si.empty()) { si.push_back(0); }
            return false;
        }
        if (si.empty()) {
            diag_index.clear();
            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(); }
            index.reserve(nnz);
            m.reserve(nnz);
            for (begin = value.begin(); begin != end; ++begin) {
                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, Mcompressed_col_update_inv>*>(this);
        noconst_this->value_to_ccarray();
        if (solver == 0) {
            noconst_this->solver = LU_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);
    }
 
    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);
                }
            }
        }
    }
 
    size_t s1;
    std::vector<size_t> si;     
    std::vector<T> m;           
    std::vector<size_t> index;  
    std::vector<size_t> diag_index;
    std::vector<std::vector<std::pair<size_t, T>>>
          value;  
    LU_sparse_solver<T>* solver;
    SparseMatrixConnectivityType connectivity;
    const Dictionary* dictionary;
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mcompressed_col_update_inv> Matrix<T, Mcompressed_col_update_inv>::null;
 
struct Mcompressed_col_update_inv_ext {
    using copy = Mcompressed_col_update_inv;
    using inverse = Mcompressed_col_update_inv;
};
 
template <typename T>
class Matrix<T, Mcompressed_col_update_inv_ext> : public Matrix<T, Mcompressed_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, Mcompressed_col_update_inv>(s1_, connectivity_, dictionary_),
          size_ext(size_ext_),
          solver(0) {}
 
    Matrix(const Matrix& m_) : size_ext(m_.size_ext), solver(0) {}
 
    ~Matrix() { delete solver; }
 
    std::pair<size_t, size_t> size() const {
        return std::pair<size_t, size_t>(
              Matrix<T, Mcompressed_col_update_inv>::size1() + size_ext,
              Matrix<T, Mcompressed_col_update_inv>::size2() + size_ext);
    }
 
    size_t size1() const { return Matrix<T, Mcompressed_col_update_inv>::size1() + size_ext; }
 
    size_t size2() const { return Matrix<T, Mcompressed_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, Mcompressed_col_update_inv>::resize(s, connectivity_, dictionary_);
        size_ext = s_ext;
    }
 
    void set_zero() {
        Matrix<T, Mcompressed_col_update_inv>::set_zero();
        delete solver;
        solver = 0;
    }
 
    void set_zero_same_struct() {
        Matrix<T, Mcompressed_col_update_inv>::set_zero_same_struct();
        if (solver) { solver->update_value(); }
    }
 
    bool is_null() const { return this == &null; }
 
    T operator()(size_t i, size_t j) const {
        if (i < Matrix<T, Mcompressed_col_update_inv>::size1()
            && j < Matrix<T, Mcompressed_col_update_inv>::size2()) {
            return Matrix<T, Mcompressed_col_update_inv>::operator()(i, j);
        }
        return 0;
    }
 
    Matrix<T, Mcompressed_col_update_sub_ref> operator()(const Interval& i, const Interval& j) {
        return Matrix<T, Mcompressed_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, Mcompressed_col_update_inv_ext>*>(this);
        if (noconst_this->value_to_ccarray()) {
            delete noconst_this->solver;
            noconst_this->solver = 0;
        }
        if (solver == 0) {
            noconst_this->solver = LU_extension_sparse_solver<T>::new_default(
                  Matrix<T, Mcompressed_col_update_inv>::connectivity,
                  *Matrix<T, Mcompressed_col_update_inv>::dictionary);
            noconst_this->solver->init(
                  Matrix<T, Mcompressed_col_update_inv>::si.size() - 1,
                  Matrix<T, Mcompressed_col_update_inv>::index.size(),
                  &Matrix<T, Mcompressed_col_update_inv>::si[0],
                  &Matrix<T, Mcompressed_col_update_inv>::index[0],
                  &Matrix<T, Mcompressed_col_update_inv>::m[0], size_ext,
                  Matrix<T, Mcompressed_col_update_inv>::connectivity,
                  *Matrix<T, Mcompressed_col_update_inv>::dictionary);
        }
        noconst_this->solver->resolve(s, nrhs, b, ldb, x, ldx, ma, mb, mc, left_or_right);
    }
 
    static Matrix null;
 
private:
    size_t size_ext;
    LU_extension_sparse_solver<T>* solver;
    MVFRIEND;
};
 
template <typename T>
Matrix<T, Mcompressed_col_update_inv_ext> Matrix<T, Mcompressed_col_update_inv_ext>::null;
 
struct Mcompressed_col_update_sub_ref {
    using copy = Mrectangle;
};
 
template <typename T>
class Matrix<T, Mcompressed_col_update_sub_ref> {
public:
    Matrix(Matrix<T, Mcompressed_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;
        }
        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) {
                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 (!(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 + j[0]].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) {
            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, Mcompressed_col_update_inv>& m;
    const Interval& i;
    const Interval& j;
    MVFRIEND;
};
 
}  // namespace b2000::b2linalg
 
#endif