【发布时间】：2017-11-18 16:24:20
【问题描述】：
我正在实现一个修改后的压缩稀疏行矩阵[reference]， 但是我对 Matrix * 向量乘法有疑问，我编写了函数但我没有找到错误！
该类使用 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]，所以列索引在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
但使用 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 代码编写的完全相同的向量
感谢我听从了@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）
【问题讨论】：
标签： c++ matrix sparse-matrix