矩阵Mod.稀疏行运算符*(矩阵，向量)

Question

我正在实现一个修改的压缩稀疏行矩阵[reference]，但我对矩阵*向量乘法有问题，我写了函数，但我没有找到错误！

该类使用2个容器(std:：vector)来存储

• 对角线元素(aa_[0]到aa_[dim])
• 非零值偏离对角线(aa_[dim+2]到aa_[size_of_non_zero])
• 行中第一个元素的指针(ja_[0]到ja_[dim])
• 在前一个指针中，使用了以下规则：ja_[0]=dim+1；ja_[i+1]-ja[i]=第i行元素数
• 对于上述ja_[row]，存储在ja_[ja_[row]]中的列索引的范围是ja[0]到ja[dim+1]，因此列索引在ja_[dim+2]到ja_[size_of_non_zero elment]中

这里的最小代码：

# include <initializer_list>
# include <vector>
# include <iosfwd>
# include <string>
# include <cstdlib>
# include <cassert>
# include <iomanip>
# include <cmath> for(auto i=0; i< A.dim ; i++)
{
//for(auto k=A.ja_.at(i) ; k <= A.ja_.at(i+1)-1 ; k++ )
auto k=A.ja_.at(i)-1; 
do 
{    
b.at(i) += A.aa_.at(k)* x.at(A.ja_.at(k)-1);
k++ ; for(auto i=0; i< A.dim ; i++)
{
//for(auto k=A.ja_.at(i) ; k <= A.ja_.at(i+1)-1 ; k++ )
auto k=A.ja_.at(i)-1; 
do 
{    
b.at(i) += A.aa_.at(k)* x.at(A.ja_.at(k)-1);
k++ ;
}while (k < A.ja_.at(i+1)-1 ); // ;
}
return b;
}while (k < A.ja_.at(i+1)-1 ); // ;
}
return b;
# include <set>
# include <fstream>

template <typename data_type>
class MCSRmatrix {
public:
using itype = std::size_t ;
template <typename T>
friend std::vector<T> operator*(const MCSRmatrix<T>& A, const std::vector<T>& x ) noexcept ;
public:
constexpr MCSRmatrix( std::initializer_list<std::initializer_list<data_type>> rows);

private:
std::vector<data_type> aa_ ;    // vector of value 
std::vector<itype>     ja_ ;    // pointer vector 
int dim ; 
};
//constructor 
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()) for(auto i=0; i< A.dim ; i++)
{
//for(auto k=A.ja_.at(i) ; k <= A.ja_.at(i+1)-1 ; k++ )
auto k=A.ja_.at(i)-1; 
do 
{    
b.at(i) += A.aa_.at(k)* x.at(A.ja_.at(k)-1);
k++ ;
}while (k < A.ja_.at(i+1)-1 ); // ;
}
return b;
{
throw std::runtime_error("error matrix must be square");
}
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;           
}
}     
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>
std::vector<T> operator*(const MCSRmatrix<T>& A, const std::vector<T>& x ) noexcept 
{     
std::vector<T> b(A.dim); 
for(auto i=0; i < A.dim ; i++ )
b.at(i) = A.aa_.at(i)* x.at(i) ;   

for(auto i=0; i< A.dim ; i++)
{
for(auto k=A.ja_.at(i) ; k < A.ja_.at(i+1)-1 ; k++ )
{    
b.at(i) += A.aa_.at(k)* x.at(A.ja_.at(k));
}   
}
return b;
}

最后是主

# include "ModCSRmatrix.H"

using namespace std;
int main(){
std::vector<double> v1={0,1.3,4.2,0.8};
MCSRmatrix<double> m1  = {{1.01, 0 , 2.34,0}, {0, 4.07, 0,0},{3.12,0,6.08,0},{1.06,0,2.2,9.9} }; 
std::vector<double> v2 = m1*v1 ;
for(auto& x : v2)
cout << x << ' ' ;
cout << endl;
}

但是结果与八度音阶的结果不同！

我已经更正了代码，现在编译！它给了我一个结果：

0 5.291 25.536 9.68

但是使用倍频程获得的正确结果是：

9.8280 5.2910 25.5360 17.1600

奇怪的是，用Fortran编写的同一个代码也能正常工作！

MODULE MSR
IMPLICIT NONE
CONTAINS
subroutine amuxms (n, x, y, a,ja)
real*8  x(*), y(*), a(*)
integer n, ja(*)
integer i, k
do 10 i=1, n
y(i) = a(i)*x(i)
10     continue
do 100 i = 1,n
do 99 k=ja(i), ja(i+1)-1
y(i) = y(i) + a(k) *x(ja(k))
99      continue
100  continue
return
end
END MODULE
PROGRAM MSRtest
USE MSR
IMPLICIT NONE
INTEGER :: i
REAL(KIND(0.D0)), DIMENSION(4) :: y, x= (/0.,1.3,4.2,0.8/)
REAL(KIND(0.D0)), DIMENSION(9) :: AA = (/ 1.01, 4.07, 6.08, 9.9, 0., 2.34, 3.12, 1.06, 2.2/) 
INTEGER , DIMENSION(9)         :: JA = (/6, 7, 7, 8, 10, 3, 1, 1, 3/)

WRITE(6,FMT='(4F8.3)') (x(I), I=1,4)    
CALL amuxms(4,x,y,aa,ja)
WRITE(6,FMT='(4F8.3)') (y(I), I=1,4)    
END PROGRAM

在上面的代码中，aa和ja的值是由放置该成员的c++构造函数给定的

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;
}

并在构造函数的末尾调用它！现在我已经在构造函数的末尾添加了成员的行，所以如果你尝试构造函数，你会得到与fortran代码中写的完全相同的向量

谢谢，我听从了Paul H.的建议，并将运算符+重写如下：(我没有改变ja_索引，因为在我的课上我有很多或多或少没有错误的方法)

template <typename T>
std::vector<T> operator*(const MCSRmatrix<T>& A, const std::vector<T>& x ) noexcept 
{     
std::vector<T> b(A.dim); 
for(auto i=0; i < A.dim ; i++ )
b.at(i) = A.aa_.at(i)* x.at(i) ;   

for(auto i=0; i< A.dim ; i++)
{
//for(auto k=A.ja_.at(i) ; k <= A.ja_.at(i+1)-1 ; k++ )
auto k=A.ja_.at(i)-1; 
do 
{    
b.at(i) += A.aa_.at(k)* x.at(A.ja_.at(k)-1);
k++ ;
}while (k < A.ja_.at(i+1)-1 ); // ;
}
return b;
}

正如你所看到的，我用从所有ja_中减去1

• x.at(A.ja_.at(k)-1)而非x.at(A.ja_.at(k))
• 索引K的不同起始k=A.ja_.at(i)-1
• 和cicle的不同结尾(我用了do while而不是for)

Answer 1

调试器是您的朋友！为了将来参考，这里有一个链接到一篇关于调试小程序的非常好的博客文章：如何调试小程序。你的代码中有几个错误。如果您在链接到的引用中创建用作示例的4 x 4矩阵，您将看到您计算的ja_值均为1。Fortran版本之所以有效，是因为Fortran中的数组默认情况下是从1开始索引的，而不是从0开始索引的。因此在class MCSRmatrix中更改

ja_.at(w) = dim+2;

至

ja_.at(w) = dim+1;

和

ja_.push_back(j);

至

ja_.push_back(j-1);

然后在您的operator*方法中更改

for(auto k=A.ja_.at(i) ; k < A.ja_.at(i+1)-1 ; k++ )

至

for(auto k = A.ja_.at(i); k < A.ja_.at(i+1); k++)