大家好，我实现了一个存储在修改压缩稀疏行中的稀疏矩阵！构造函数工作正常，我可以验证它，但 operator+ 有一个奇怪的行为：如果我有一个非零值，则总和不会计算正确的结果。此处描述了修改后的压缩稀疏行方法， 我的最小工作代码如下：
# 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
另一个示例，如果我将这一行放入 main ：
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'