Math(1)---Eigen稀疏矩阵乘法

Eigen稀疏矩阵乘法

• 稀疏矩阵能够节省存储空间；
• Eigen中稀疏矩阵的存储方式：CRS, CCS;
• 稀疏矩阵缩短遍历元素的时间。

Eigen稀疏矩阵乘以向量

• 计算公式： (MatResult = Matsparse*Vec)
• 利用Eigen可以直接根据公式计算出来，但是笔者想弄楚，Eigen是怎样实现的，于是用迭代法实现计算

示例：

#include <Eigen/Sparse>
#include <vector>
#include <iostream>

using namespace Eigen;
using namespace std;

SparseMatrix<int,RowMajor> SpMat(5,5);
int main(void)
{
	MatrixXi A(5,5);
	VectorXi V(5);
	VectorXi buffer(5);
	VectorXi Ibus(5);
	buffer << 0, 0, 0, 0, 0;
	Ibus << 0, 0, 0, 0, 0;
	V << 1, 2, 3, 4, 5;
	A << 0, 3, 0, 0, 0,
		22, 0, 0, 0, 17,
		7, 5, 0, 1, 0,
		0, 0, 0, 0, 0,
		0, 0, 14, 0, 8;
	SpMat = A.sparseView();
	
	for (int k = 0; k < SpMat.outerSize(); ++k)
	{
		for (SparseMatrix<int, RowMajor>::InnerIterator it(SpMat, k); it; ++it)
		{
			//cout<<it.value()<<endl;  //Eigen 存储稀疏矩阵的值，默认是CCS，这里指定成了CRS
			//cout << it.row() << endl; //Eigen 行索引
			//cout << it.col() << endl; //Eigen 列索引
			//cout << it.index() << endl; //这里index=it.col()

			buffer(k) += it.value()*V(it.index());
		}
		Ibus(k) += buffer(k);
	}

	cout << SpMat << endl;
	cout <<	Ibus << endl;
	return 0;
}

CRS存储格式下两种实现方式

• 1.不采用Eigen库函数

#include<iostream>

using namespace std;

int main(void)
{
	int Yx[] = { 3, 22, 17, 7, 5, 1, 14, 8 }; //稀疏矩阵非零元
	int Yj[] = { 1, 0, 4, 0, 1, 3, 2, 4 };    //稀疏矩阵非零元素对应的列号
	int Yp[] = { 0, 1, 3, 6, 6, 8 };          //稀疏矩阵非零元素的偏移量
	int buffer[5] = { 0 };
	int Ibus[5] = { 0 };
	int V[] = { 1, 2, 3, 4, 5 };
	for (int r = 0; r < 5; ++r)
	{
		for (int k = Yp[r]; k < Yp[r + 1]; ++k)
		{
			buffer[r] += Yx[k] * V[Yj[k]];
		}
		Ibus[r] += buffer[r];
	}
	for (int i = 0; i < 5; ++i)
	{
		cout << Ibus[i] << endl;
	}
	return 0;
}

• 2.Eigen库函数的另一种实现方式

#include <Eigen/Dense>
#include <Eigen/Sparse>

#include <iostream>

using namespace std;
using namespace Eigen;


int main()
{
	SparseMatrix<double, RowMajor> mat(5, 5);
	
	MatrixXd A(5, 5);
	VectorXi V(5);
	VectorXi Yx(8);
	VectorXi Yj(8);
	VectorXi Yp(6);
	VectorXi buffer(5);
	VectorXi Ibus(5);
	Yx << 3, 22, 17, 7, 5, 1, 14, 8;
	Yj << 1, 0, 4, 0, 1, 3, 2, 4;
	Yp << 0, 1, 3, 6, 6, 8;
	V << 1, 2, 3, 4, 5;
	buffer << 0,0,0,0,0;
	Ibus << 0, 0, 0, 0, 0;

	//A matrix is for checking
	/*A << 0, 3, 0, 0, 0,
		22, 0, 0, 0, 17,
		7, 5, 0, 1, 0,
		0, 0, 0, 0, 0,
		0, 0, 14, 0, 8;*/
	
	for (int r = 0; r < Yp.size() - 1; ++r)
	{
		for (int k = Yp(r); k <= Yp(r + 1)-1; ++k) {
			buffer(r) += Yx(k)*V(Yj(k));
		}
		Ibus(r) += buffer(r);
	}
	cout << Ibus << endl;
	
}