稀疏矩阵的运算
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了稀疏矩阵的运算相关的知识，希望对你有一定的参考价值。
参考技术A 内容假设两个稀疏矩阵A和B，他们均为m行n列，要求表写求矩阵的加法即：C=A+B的算法（C矩阵存储A与B相加的结果）
分析
利用一维数组来存储，一维数组顺序存放非零元素的行号、列号和数值，行号-1表示结束，然后进行矩阵加法运算时依次扫描矩阵A和B的行列值，并以行优先。当行列相同的时候，将第三个元素的值相加和以及行列号三个元素存入结果数组C中；不相同时，将A或B的三个元素直接存入结果数组中。
代码
// fanchen.cpp : 定义控制台应用程序的入口点。
//
#include "stdafx.h"
#include<iostream>
using namespace std;
#define LINE 10
struct Node
//行号
int row;
//列号
int line;
//数据
int data;
;
//初始化数组
int init(Node array[])
int row,line,data;
int index = 0;
while(cin>>row)
if(row == -1)
break;
cin>>line>>data;
array[index].data=data;
array[index].line = line;
array[index].row = row;
index++;
return index;
//打印数组
void printArray(Node array[],int len)
for(int i = 0;i < len;i++)
cout<<array[i].row<<" "<<array[i].line<<" "<<array[i].data<<endl;
int calc(Node a[],int lena,Node b[],int lenb,Node c[])
//计数单位
int index = 0;
int i = 0,j = 0,k;
while(i < lena && j < lenb)
int tempa = a[i].row * 10 + a[i].line;
int tempb = b[j].row * 10 + b[j].line;
if(tempa < tempb)
c[index++] = a[i++];
else if(tempa == tempb)
int l = a[i].line;
int r = a[i].row;
int data = a[i].data + b[j].data;
c[index].line = l;
c[index].row = r;
c[index].data = data;
index++;
i++;
j++;
else
c[index++] = b[j++];
while(i < lena)
c[index++] = a[i++];
while(j < lenb)
c[index++] = b[j++];
return index;
int _tmain(int argc, _TCHAR* argv[])
Node a[60],b[60],c[60];
//a的长度
int lena = init(a);
int lenb = init(b);
int lenc = calc(a,lena,b,lenb,c);
cout<<"a::::"<<endl;
printArray(a,lena);
cout<<endl;
cout<<"b::::"<<endl;
printArray(b,lenb);
cout<<endl;
printArray(c,lenc);
system("pause");
return 0;
参考技术B
M AT L A B中对满矩阵的运算和函数同样可用在稀疏矩阵中.结果是稀疏矩阵还是满矩阵,
这取决于运算符或者函数及下列的操作数:当函数用一个矩阵作为输入参数,输出参数为一个标量或者一个给定大小的向量时,输出参数的格式总是返回一个满阵形式,如命令s i z e.
当函数用一个标量或者一个向量作为输入参数,输出参数为一个矩阵时,输出参数的格式也总是返回一个满矩阵,如命令e y e.还有一些特殊的命令可以得到稀疏矩阵,如命令s p e y e.
对于单参数的其他函数来说,通常返回的结果和参数的形式是一样的,如d i a g.
对于双参数的运算或者函数来说,如果两个参数的形式一样,那么也返回同样形式的结果.在两个参数形式不一样的情况下,除非运算的需要,均以满矩阵的形式给出结果.
两个矩阵的组和[A B],如果A或B中至少有一个是满矩阵,则得到的结果就是满矩阵.
表达式右边的冒号是要求一个参数的运算符,遵守这些运算规则.
表达式左边的冒号不改变矩阵的形式. 假设有:
这是一个5×5的单位满矩阵和相应的稀疏矩阵.
(a) C = 5*B,结果为:
这是一个稀疏矩阵.
(b) D = A + B,给出的结果为:
这是一个满矩阵.
(c) x = B \\ h,结果为:
这是一个满向量.
有许多命令可以对非零元素进行操作.
命令集8 9矩阵的非零元素
n n z ( A )求矩阵A中非零元素的个数.它既可求满矩阵也可求稀疏矩阵.
s p y ( A )画出稀疏矩阵A中非零元素的分布.也可用在满矩阵中,在
这种情况下,只给出非零元素的分布.
s p y ( A , c s t r , s i z e )用指定的颜色c s t r(见表1 3 - 1 )和在s i z e规定的范围内画出稀疏
矩阵A中非零元素的分布.
n o n z e r o s ( A )按照列的顺序找出矩阵A中非零的元素.
s p o n e s ( A )把矩阵A中的非零元素全换为1.
s p a l l o c ( m , n ,产生一个m×n阶只有n z m a x个非零元素的稀疏矩阵.这样可以
n z m a x )有效地减少存储空间和提高运算速度.
n z m a x ( A )给出为矩阵A中非零元素分配的内存数.不一定和n n z ( A )得
到的数相同;参见s p a r s e或者s p a l l o c.
i s s p a r s e ( A )如果矩阵A是稀疏矩阵,则返回1;否则返回0.
s p f u n ( f c n , A )用A中所有非零元素对函数f c n求值,如果函数不是对稀疏矩
阵定义的,同样也可以求值.
s p r a n k( A )求稀疏矩阵A的结构秩.对于所有的矩阵来说,都有
s p r a n k ( A）≥r a n k ( A ). 用下面的命令定义稀疏矩阵:
创建一个大矩阵:
Big=kron(A, A)
这个矩阵B i g是什么样子呢?
K r o n e c k e r张量积给出一个大矩阵,它的元素是矩阵A的元素之间可能的乘积.因为参量都是稀疏矩阵,所以得到的矩阵也是一个稀疏矩阵.可以用命令 w h o s和i s s p a r s e来确认一下.
查看矩阵B i g的结构图,可输入s p y ( B i g ),结构图如右图所示. 从图中可以看出B i g是一个块双对角矩阵. MATLAB中有四个基本稀疏矩阵,它们是单位矩阵,随机矩阵,对称随机矩阵和对角矩阵.
命令集9 0单位稀疏矩阵
s p e y e ( n )生成n×n的单位稀疏矩阵.
s p e y e ( m , n )生成m×n的单位稀疏矩阵.
命令speye(A) 得到的结果和s p a r s e ( e y e ( A ) )是一样的,但是没有涉及到满阵的存储.
命令集9 1随机稀疏矩阵
s p r a n d ( A )生成与A有相同结构的随机稀疏矩阵,且元素服从均匀分布.
s p r a n d ( m , n , d e n s )生成一个m×n的服从均匀分布的随机稀疏矩阵,有d e n s×m×
n个非零元素,0≤d e n s≤1.参数d e n s是非零元素的分布密度.
s p r a n d ( m , n , d e n s ,生成一个近似的条件数为1 /rc,大小为m×n的随机稀疏矩
r c )阵.如果rc=rc是一个长度为l≤l ( m i n (m, n) )的向量,那么
矩阵将rci作为它l个奇异值的第一个,其他的奇异值为0.
s p r a n d n ( A )生成与A有相同结构的随机稀疏矩阵,且元素服从正态分布.
s p r a n d n ( m , n , d e n s ,生成一个m×n的服从正态分布的随机稀疏矩阵,和sprand
r c )一样.
s p r a n d s y m ( S )生成一个随机对称稀疏矩阵.它的下三角及主对角线部分与S的结构相同,矩阵元素服从正态分布.
s p r a n d s y m ( n , d e n s )生成一个m×n的随机对称稀疏矩阵.矩阵元素服从正态分布,分布密度为d e n s.
s p r a n d s y m ( n , d e n s ,生成一个近似条件数为1 /rc的随机对称稀疏矩阵.元素以0r c )对称分布,但不是正态分布.如果rc=rc是一个向量,则矩阵有特征值rci.也就是说,如果rc是一个正向量,则矩阵是正定矩阵.
s p r a n d s y m ( n , d e n s ,生成一个正定矩阵.如果k= 1,则矩阵是由一正定对称矩r c , k )阵经随机J a c o b i旋转得到的,其条件数正好等于1 /rc;如果k= 2,则矩阵为外积的换位和,其条件数近似等于1 /rc.
s p r a n d s y m ( S , d e n s ,生成一个与矩阵S结构相同的稀疏矩阵,近似条件数为1 /rc.
r c , 3 )参数d e n s被忽略,但是这个参数在这个位置以便函数能确认最后两个参数的正确与否. (a) 假设有矩阵A:
输入R a n d o m = s p r a n d n ( A ),可得到随机稀疏矩阵:
矩阵中随机数的位置和矩阵A中非零元素的位置相同.
(b) 对于( a )中的矩阵A,输入:
B = s p r a n d s y m ( A )
结果为:
这是一个用矩阵A的下三角及主对角线部分创建的对称矩阵,在非零元素的位置用随机数作为元素值.
用命令s p d i a g s可以取出对角线元素,并创建带状对角矩阵.假设矩阵A的大小为m×n,
在p个对角线上有非零元素.B的大小为m i n (m×n)×p,它的列是矩阵A的对角线.向量d的长度为p,其整型分量给定了A的对角元:
di0 主对角线上的对角线
命令集9 2对角稀疏矩阵
[ B , d ] = s p d i a g s ( A )求出A中所有的对角元,对角元保存在矩阵B中,它们的下标保存在向量d中.
s p d i a g s ( A , d )生成一个矩阵,这个矩阵包含有矩阵A中向量d规定的对角元.
s p d i a g s ( B , d , A )生成矩阵A,用矩阵B中的列替换d定义的对角元.
A = s p d i a g s ( B , d , m , n )用保存在由d定义的B中的对角元创建稀疏矩阵A.
例11 . 4给出了如何使用s p d i a g s命令来解普通微分方程组. 在许多实际应用中要保留稀疏矩阵的结构,但是在计算过程中的中间结果会减弱它的稀疏性,如L U分解.这就会导致增加浮点运算次数和存储空间.为了避免这种情况发生,在第稀疏矩阵1 2 9
M AT L A B中用命令对矩阵进行重新安排.这些命令都列在下面的命令集9 3中.通过h e l p命令
可以得到每个命令更多的帮助信息,也可见h e l p d e s k.
命令集9 3矩阵变换
c o l m m d ( A )返回一个变换向量,使得矩阵A列的秩为最小.
s y m m m d ( A )返回使对称矩阵秩为最小的变换.
s y m r c m ( A )矩阵A的C u t h i l l - M c K e e逆变换.矩阵A的非零元素在主对角线附近.
c o l p e r m ( A )返回一个矩阵A的列变换的向量.列按非零元素升序排列.有时这是L U因式分解前有用的变换:lu(A(:, j)).如果A是一个对称矩阵,对行和列进行排序,这有利于C h o l e s k y分解:chol(A(j, j)).
r a n d p e r m ( n )给出正数1,2,. . .,n的随机排列,可以用来创建随机变换矩阵.
d m p e r m ( A )对矩阵A进行D u l m a g e - M e n d e l s o h n分解,输入help dmperm可得更多信息. 创建一个秩为4的变换矩阵,可输入:
一旦运行p e r m = r a n d p e r m ( 4 ),就会得到:
给出的变换矩阵为:
如果矩阵A为:
输入命令:
运行结果为:
有两个不完全因式分解命令,它们是用来在解大线性方程组前进行预处理的.用h e l p d e s k命令可得更多信息.命令集9 4不完全因式分解c h o l i n c ( A , o p t )进行不完全C h o l e s k y分解,变量o p t取下列值之一:
d r o p t o l指定不完全分解的舍入误差,0给出完全分解.
m i c h o l如果m i c h o l = 1,则从对角线上抽取出被去掉的元素.
r d i a g用s q r t ( d r o p t o l*n o r m ( X ( : , j ) ) )代替上三角分
解因子中的零元素,j为零元素所在的列.
[ L , U , P ]=返回矩阵X的不完全分解得到的三个矩阵L,U和P,变量o p t取
l u i n c ( X , o p t )下列值之一:
d r o p t o l指定分解的舍入误差.
m i l u改变分解以便从上三角角分解因子中抽取被去掉的列元素.
u d i a g用d r o p t o l值代替上三角角分解因子中的对角线上的零元素.
t h r e s h中心极限.
解稀疏线性方程组既可用左除运算符解,也可用一些特殊命令来解.
命令集9 5稀疏矩阵和线性方程组
s p p a r m s ( k e y s t r , o p )设置稀疏矩阵算法的参数,用help spparms可得详细信息.
s p a u g m e n t ( A , c )根据[ c*l A; A\' 0 ]创建稀疏矩阵,这是二次线性方程组的最
小二乘问题.参见7 . 7节.
s y m b f a c t ( A )给出稀疏矩阵的C h o l e s k y和L U因式分解的符号分解因子.
用help symbfact可得详细信息.
稀疏矩阵的范数计算和普通满矩阵的范数计算有一个重要的区别.稀疏矩阵的欧几里德范数不能直接求得.如果稀疏矩阵是一个小矩阵,则用n o r m ( f u l l ( A ) )来计算它的范数;但是对于大矩阵来说,这样计算是不可能的.然而M AT L A B可以计算出欧几里德范数的近似值,在计算条件数时也是一样.
命令集9 6稀疏矩阵的近似欧几里德范数和条件数
n o r m e s t ( A )计算A的近似欧几里德范数,相对误差为1 0-6.
n o r m e s t ( A , t o l )计算A的近似欧几里德范数,设置相对误差t o l,而不用缺省时的1 0-6.
[ n r m , n i t ] =计算近似n r m范数,还给出计算范数迭代的次数n i t.
n o r m e s t ( A )
c o n d e s t ( A )求矩阵A条件数的1 -范数中的下界估计值.
[ c , v ]=求矩阵A的1 -范数中条件数的下界估计值c和向量v,使得
c o n d e s t ( A , t r )| |Av| | = ( | |A| | | |v| | ) / c.如果给定t r,则给出计算的过程.t r= 1,
给出每步过程;t r=-1,给出商c / r c o n d ( A ). 假设给出:
用n o r m A p p r o x = n o r m e s t ( S p r s )计算出:
用t h e N o r m = n o r m ( f u l l ( S p r s ) )得:
为了找到它们之间的差别,计算d i f f e r e n c e = t h e N o r m - n o r m A p p r o x,得:
在许多应用中,n o r m e s t计算得到的近似值是一个很好的近似欧几里德范数,它的计算步数要比n o r m要少得多;可参见7 . 6节.
用e t r e e命令来找到稀疏对称矩阵的消元树,用向量f来描述消元树,还可用e t r e e p l o t命令画出来.元素fi是矩阵的上三角C h o l e s k y分解因子中i行上第1非零元素的列下标.如果有非零元素,则fi= 0.消元树可以这样来建立:
节点i是fi的孩子,或者如果fi= 0,则节点i是树的根节点.
命令集9 7矩阵的消元树
e t r e e ( A )求A的消元树向量f,这个命令有可选参数;输入h e l p
e t r e e获取帮助.
e t r e e p l o t ( A )画出向量f定义的消元树图形.
t r e e p l o t ( p , c , d )画出指针向量p的树图形,参数c和d分别指定节点的颜色和分支数.e t r e e p l o t可以调用这个命令.
t r e e l a y o u t显示树的结构,t r e e p l o t可以调用这个命令. 假设有对称稀疏矩阵B:
运行命令b t r e e = e t r e e ( B ),结果为:
开始的数字2不难理解,它是矩阵的第1列上第1个非零元素的行数,它决定了在C h o l e s k y分解因子的第1行第2列处有一个非零元素.当缩减第1列的元素时就得到第2列的数字5.B在缩减后,在( 5 , 2 )位置的元素是非零的,这样消元树向量中第2个元素的值为5.
s p y ( c h o l ( B ) )给出了C h o l e s k y分解因子的结构图,如图9 - 2所示:
图9-2 Cholesky分解结构图
图9-3 矩阵B的消元树
这个向量消元树可以这样来建立:上三角中只有一行有非零元素,节点8,因此这就是树
唯一的根.节点1是节点2的孩子,节点2和3是节点5的孩子,而节点5是节点6的孩子.节点4和6是节点7的孩子,而节点7又是节点8的孩子,即根的孩子.
命令e t r e e p l o t ( B )给出了树的结构图,如图9 - 3所示.
消元树的形状取决于列和行序,它可以用来分析消元过程.
用g p l o t命令可以画出坐标和矩阵元素间的联系图形.必须在n×2的矩阵中给出n个坐标,
矩阵的每一行作为一个点.这样就创建出点点之间连接的n×n矩阵,如果点4连接到点8,则(4, 8)的值为1.由于是一个大矩阵,而且非零元素较少,所以它应该被建成稀疏矩阵.
这个图可以说明网络问题,如传递问题.它还包含有线性方程组中未知量之间的相关信息.
命令集9 8网络图形
g p l o t ( A , K )如果矩阵A的a(i, j)不为0,则将点ki连接到点kj.K是一个n×
2的坐标矩阵,A是一个n×n的关联矩阵.
g p l o t ( A , K , s t r )用字符串s t r给定的颜色和线型画出的同上图形.字符串s t r的
取值参见表1 3 - 1.
[ X , A ] = u n m e s h ( E )求网格边界矩阵E的L a p l a c e矩阵A和网格点的坐标矩阵X.
例七
假设有下面的坐标矩阵K和关联矩阵A:
矩阵A在稀疏化后,用命令g p l o t ( A , K )画出图9 - 4,给出了点(0, 1)和点(4, 1)之间所有可能的路径.
矩阵模组。稀疏行运算符*（矩阵，向量）
【中文标题】矩阵模组。稀疏行运算符*（矩阵，向量）【英文标题】：Matrix Mod. sparse row Operator*(matrix , vector) 【发布时间】：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）
【问题讨论】：
【参考方案1】：调试器是你的朋友！为了将来参考，这里有一个关于调试小程序的非常好的博客文章的链接：How to debug small programs。
您的代码中有几个错误。如果您在链接到的参考中创建用作示例的 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++)
【讨论】：
是的，谢谢，我已经解决了使用索引从所有 ja_ 中减去 1以上是关于稀疏矩阵的运算的主要内容，如果未能解决你的问题，请参考以下文章