оператор + в разреженной матрице: я не нашел ошибку

Привет всем, я реализую разреженную матрицу, хранящуюся в разбросанной строке Modified Compressed! конструктор работает нормально, я могу это проверить, но оператор + имеет странное поведение: если бы у меня было ненулевое значение, сумма не вычисляет правильный результат.
описан модифицированный метод сжатых разреженных строк Вот
мой минимальный рабочий код следующий:

# include <initializer_list>
# include <vector>
# include <iosfwd>
# include <string>
# include <cstdlib>
# include <cassert>
# include <iomanip>
# include <cmath>
# include <set>
# include <fstream>
# include <algorithm>
# include <exception>

template <typename data_type>
class MCSRmatrix;

template <typename T>
std::ostream& operator<<(std::ostream& os ,const MCSRmatrix<T>& m) noexcept ;template <typename T>
MCSRmatrix<T> operator+(const MCSRmatrix<T>& m1, const MCSRmatrix<T>& m2 ) ;

template <typename data_type>
class MCSRmatrix {

public:
using itype = std::size_t ;

template <typename T>
friend std::ostream& operator<<(std::ostream& os ,const MCSRmatrix<T>& m) noexcept ;

template <typename T>
friend MCSRmatrix<T> operator+(const MCSRmatrix<T>& m1, const MCSRmatrix<T>& m2 ) ;

public:

constexpr MCSRmatrix( std::initializer_list<std::initializer_list<data_type>> rows);

constexpr MCSRmatrix(const std::size_t& ) noexcept;

const data_type& operator()(const itype r , const itype c) const noexcept ;

data_type operator()(const itype r , const itype c) noexcept ;

auto constexpr printMCSR() const noexcept ;

private:

std::vector<data_type> aa_ ;    // vector of value
std::vector<itype>     ja_ ;    // pointer vector

int dim ;

std::size_t constexpr findIndex(const itype row ,  const itype col) const noexcept ;
};

template <typename T>
constexpr MCSRmatrix<T>::MCSRmatrix( std::initializer_list<std::initializer_list<T>> rows)
{
this->dim  = rows.size();
auto _rows = *(rows.begin());

aa_.resize(dim+1);
ja_.resize(dim+1);

if(dim != _rows.size())
{
throw std::runtime_error("Error in costructor! MCSR format require square matrix!");
}

itype w = 0 ;
ja_.at(w) = dim+2 ;
for(auto ii = rows.begin(), i=1; ii != rows.end() ; ++ii, i++)
{
for(auto ij = ii->begin(), j=1, elemCount = 0 ; ij != ii->end() ; ++ij, j++ )
{
if(i==j)
aa_[i-1] = *ij ;
else if( i != j && *ij != 0 )
{
ja_.push_back(j);
aa_.push_back(*ij);
elemCount++ ;
}
ja_[i] = ja_[i-1] + elemCount;
}
}
printMCSR();
}
template <typename T>
constexpr MCSRmatrix<T>::MCSRmatrix(const std::size_t& n ) noexcept
{
this->dim = n;
aa_.resize(dim+1);
ja_.resize(dim+1);

for(std::size_t i = 0; i < aa_.size()-1 ; i++)
aa_.at(i) = 1 ;

for(std::size_t i = 0; i < ja_.size() ; i++)
ja_.at(i) = aa_.size()+1 ;}template <typename T>
std::size_t constexpr MCSRmatrix<T>::findIndex(const itype row ,  const itype col) const noexcept
{
assert( row > 0 && row <= dim && col > 0 && col <= dim );
if(row == col)
{
return row-1;
}
int i = -1;

for(int i = ja_.at(row-1)-1 ; i < ja_.at(row)-1 ; i++ )
{
if( ja_.at(i) == col )
{
return i ;
}
}
return -1;
}template <typename T>
inline auto constexpr MCSRmatrix<T>::printMCSR() const noexcept
{
for(auto& x : aa_ )
std::cout << x << ' ' ;
std::cout << std::endl;

for(auto& x : ja_ )
std::cout << x << ' ' ;
std::cout << std::endl;
}

template <typename T>
const T& MCSRmatrix<T>::operator()(const itype r , const itype c) const noexcept
{
auto i = findIndex(r,c);
if( i != -1 && i < aa_.size() )
{
return aa_.at(i) ;
}
else
{
return aa_.at(dim) ;
}
}template <typename T>
T MCSRmatrix<T>::operator()(const itype r , const itype c) noexcept
{
auto i = findIndex(r,c);
if( i != -1 && i < aa_.size() )
{
return aa_.at(i) ;
}
else
{
return aa_.at(dim) ;
}
}

// non member operator

template <typename T>
std::ostream& operator<<(std::ostream& os ,const MCSRmatrix<T>& m) noexcept
{
for(std::size_t i=1 ; i <= m.dim ; i++ )
{
for(std::size_t j=1 ; j <= m.dim ; j++)
os << std::setw(8) << m(i,j) << ' ' ;
os << std::endl;
}
return os;
}// perform sum by 2 Mod Comp Sparse Row matrices
template <typename T>
MCSRmatrix<T> operator+(const MCSRmatrix<T>& m1, const MCSRmatrix<T>& m2 )
{
if(m1.dim != m2.dim)
{
throw std::runtime_error("Matrixs dimension does not match! Error in operator +");
}
else
{
MCSRmatrix<T> res(m1.dim);for(auto i=0; i < res.dim ; i++)
res.aa_.at(i) = m1.aa_.at(i)  + m2.aa_.at(i) ;

res.ja_.at(0) = res.dim+2;

std::set<unsigned int> ctrl;int n1=0, n2=0, j1=0 , j2 =0;
for(auto i=0 ; i < res.dim  ; i++)
{

ctrl.clear();

n1 = m1.ja_.at(i+1)- m1.ja_.at(i) ;
n2 = m2.ja_.at(i+1)- m2.ja_.at(i) ;j1=0 , j2=0 ;

auto sum1 = 0. , sum2 = 0. , sum=0.;for(auto j = 0; j < std::max(n1,n2) ; j++ )
{
if(n1 && n2)
{
if(m1.ja_.at(m1.ja_.at(i)-1 + j1 ) == m2.ja_.at(m2.ja_.at(i)-1 + j2))
{
ctrl.insert(m1.ja_.at(m1.ja_.at(i)-1 + j1 ));

sum = m1.aa_.at(m1.ja_.at(i)-1 + j1 ) + m2.aa_.at(m2.ja_.at(i)-1 + j2) ;}
else if(m1.ja_.at(m1.dim+1 + j1 ) != m2.ja_.at(m2.dim+1 + j2))
{
ctrl.insert(m1.ja_.at(m1.ja_.at(i)-1 + j1 ));
ctrl.insert(m2.ja_.at(m1.ja_.at(i)-1 + j2 ));

sum1 = m1.aa_.at(m1.ja_.at(i)-1 + j1);
sum2 = m2.aa_.at(m2.ja_.at(i)-1 + j1);
}
}
else if(n1)
{
ctrl.insert(m1.ja_.at(m1.ja_.at(i)-1 + j1 ));
sum1 = m1.aa_.at(m1.ja_.at(i)-1 + j1);
}
else if(n2)
{
ctrl.insert(m2.ja_.at(m2.ja_.at(i)-1 + j2 ));
sum2 = m2.aa_.at(m2.ja_.at(i)-1 + j2);
}
if(sum1)
{
res.aa_.push_back(sum1);
res.ja_.push_back(m1.ja_.at(m1.ja_.at(i)-1 + j1)) ;
}
if(sum2)
{
res.aa_.push_back(sum2);
res.ja_.push_back(m2.ja_.at(m2.ja_.at(i)-1 + j2 ));
}
if(sum)
{

res.aa_.push_back(sum);
res.ja_.push_back(m1.ja_.at(m1.ja_.at(i)-1 + j1 ));
}

if(j1 < n1) j1++ ;
if(j2 < n2) j2++ ;
}res.ja_.at(i+1) = res.ja_.at(i) + ctrl.size() ;

}
return res ;
}
}

это основная программа:

# include "ModifiedCSRmatrix.H"
using namespace std ;int main() {

MCSRmatrix<int> m01 =  {{0,1,0,0},{0,3,0,0},{0,0,0,3},{0,0,0,1}};
MCSRmatrix<int> m02 =  {{1,1,0,14},{0,22,0,3},{0,0,33,34},{0,0,1,44}};

cout << endl;
MCSRmatrix<int> m03 = m01+m02 ;

cout << m03 << endl;
cout << endl;
m03.printMCSR();return 0;
}

именно эта матрица 2 дает мне правильный результат! то есть :

1        2        0       17
0       25        0        3
0        0       33       37
0        0        1       45

но если я изменю значение элемента (3,2) в m01 матрица ..
и так матрица стала:

 MCSRmatrix<int> m01 =  {{0,1,0,0},{0,3,0,0},{0,2,0,3},{0,0,0,1}};

код дает мне этот неверный результат:

  1        2        0        0
0       25        0       14
0        2       33        0
0        0        0       45

но смотрит на вектор aa_ (вектор значения, в котором первый matrix.dim являются элементом диагонали, а последующим элементом является ненулевой ненулевой элемент, который выглядит так, что матрица отличается от этого результата)
они есть :

aa_ = 1 25 33 45 0 2 2 14 2 3 3 1 1
ja_ = 6 8 9 10 11 2 2 4 2 4 4 3 3

еще один пример, если у меня была эта строка в основной:

 MCSRmatrix<double> m3  = {{1.01, 0 , 0,0}, {0, 4.07, 0,0},{0,0,6.08,0},{1.06,0,2.2,9.9} };

MCSRmatrix<double> m4  = {{ 1.21, 0 , 1.06,0 }, {0, 3.31, 1.06,0},{0,1.06,-6.01,0},{4.12,0,1.06,-8.3}};

MCSRmatrix<double> m5 = m3+m4 ;

cout << m5 ;

дай мне правильный результат:

    2.22        0     1.06        0
0     7.38     1.06        0
0     1.06     0.07        0
5.18        0     3.26      1.6

если я, например, немного изменю m1 и m2:

 MCSRmatrix<double> m6  = {{1.01, 0 , 0,0}, {0, 4.07, 0,0},{0,0,6.08,0},{1.06,0,2.2,9.9} };
MCSRmatrix<double> m7  = {{ 1.21, 0 , 0,2.15 }, {0, 3.31, 1.03,0},{0,1.06,-6.01,1.01},{4.12,1.1,1.06,-8.3}};
cout << m6+m7 ;

эта программа выхода с terminate called after throwing an instance of 'std::out_of_range'