我正在实现修改后的压缩稀疏行矩阵 [reference] , 但是我有 Matrix * vector 乘法的问题，我写了这个函数但是我没有找到错误!
该类使用 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_[ja_[row]]中的ja_[row]的列索引的范围是ja[0]到ja[dim+1],所以colum索引在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;
}
但结果与 Octave 得到的结果不一样!
我已经更正了代码，现在可以编译了!它给了我结果:
0 5.291 25.536 9.68
但是使用 Octave 获得的正确结果是:
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 代码编写的完全相同的 vector
谢谢我听从了你的建议@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)
最佳答案
调试器是你的 friend !为了将来引用，这里有一个关于调试小程序的非常好的博客文章的链接:How to debug small programs .
您的代码中有几个错误。如果您在链接到的引用中创建用作示例的 4 x 4 矩阵，您将看到 ja_您计算的值都相差一个。 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++)
关于c++ - 矩阵模式。稀疏行运算符*(矩阵， vector )，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/47368287/