【C++算法】稀疏矩阵的乘法实现（行逻辑链接法）

算法背景：
①稀疏矩阵的类表示法：

template <class T>
class sparseMatrix
{
private:
	int rows;
	int cols;   //这里的rows和cols存的是矩阵理论有的行列数
	matrixTerm<T>* terms;  //这里的terms存的是实际的矩阵结构
	//每一个空间存的都是行数、列数、相应的元素值
	int size;
public:
	sparseMatrix() : rows(0), cols(0), size(0) {}
	sparseMatrix(int row, int col) : rows(row), cols(col), size(0) {}
	void multiply(const sparseMatrix<T>& m);
	sparseMatrix<T> operator=(const sparseMatrix<T>& m);
};

②三元组表的结构体表示法：

template <class T>
struct matrixTerm
{
	int row;
	int col;
	T value;
	operator T() const { return value; }
};

主要算法：行逻辑链接法。前提条件是三元组表按照行主索引法进行映射。同稀疏矩阵的转置操作一样，乘法操作也需要定义两个辅助的数组row_nums[i]和row_first_position[i]，分别表示矩阵A第i行的非零元素数目和第i行第1个非零元素所在三元组表中对应的索引。

int* row_nums = new int[rows + 1];   //各行的元素数目
int* row_first_position = new int[rows + 1];  //各行第一个非零元对应的索引

对数组row_nums的赋值如下：

for (int i = 0; i <= rows; i++)
		row_nums[i] = 0;   //初始化row_nums
for (int i = 0; i < size; i++)
		row_nums[terms[i].row]++;     //赋值row_nums

terms[i]表示矩阵A三元组表第i号索引的元素，terms[i].row表示矩阵A这行存在元素。
对数组row_first_position的赋值如下：

for (int i = 0; i <= rows; i++)   //如果第i行没有元素，则置为-1
	row_first_position[i] = -1;
row_first_position[1] = 0;       //第一行第一个非零元一定在0号索引
for (int i = 2; i <= rows; i++)   //遍历三元组表
	row_first_position[i] = row_first_position[i - 1] + row_nums[i - 1];

for循环中从第二列开始逻辑上“遍历”矩阵A，分析一下i = 2时，position[2]为第2列第一个非零元应在的索引位置，它显然等于第一行第一个非零元所应在的索引位置 + 第一行的非零元个数，这也是逻辑上显然成立的。对于之后的每一个i，也显然都有此行第一个非零元所应在的索引位置 = 上一行第一个非零元所应在的索引位置 + 上一行的非零元素个数。
因此，公式row_first_position[i] = row_first_position[i - 1] + row_nums[i - 1]是显然的。
然后，初始化累加器：

T* accumulator = new T[m.cols + 1];      //B有m.cols列，所以定义m.cols大小的累加器
for (int i = 0; i <= m.cols; i++)
	accumulator[i] = 0;          //初始化累加器,分别表示各个列

就可以进入for循环，进入矩阵的乘法了。
在此之前，先看一下累加器的使用规则：
假设有m行k列矩阵A和k行n列的矩阵B，将二者相乘放到矩阵C中，矩阵C是m行n列的。假设有A、B矩阵如下所示：

显然结果的C矩阵是5行5列的。根据矩阵的乘法公式，先列出计算C矩阵第一行的式子：


接着进行下一次计算时，把后面的一项分别存到累加器中：

对于矩阵A的某行的第一个非零元素，当遍历到它时，需要将它执行B.cols次乘法，即累加器的大小是B.cols + 1，因为不存在第0列。遍历到A11时，需要全部扫描一遍B的三元组表，把那些“矩阵A的列索引 = 矩阵B的行索引”的元素与A11相乘，存储到第B.terms[j].col号迭代器中去。然后，A11使用完毕，此时把row_first_position[i]++，表示第i行第一个元素所在三元组表的索引+1（因为三元组表是按照行主映射的，同一行相邻两个非零元的索引一定相邻）。也即，此时遍历到A12了。对于A12，也执行相同的算法，把那些“矩阵A的列索引 = 矩阵B的行索引”的元素与A12相乘，累加到第B.terms[j].col号迭代器中去。
显然，执行乘法需要两个for循环，外层循环“逻辑上”遍历A矩阵的某一行的非零元素，内层循环遍历B矩阵的三元组表，筛选那些符合条件的元素并累加到累加器中。

for (int k = 1; k <= row_nums[i]; k++) //逻辑上遍历A的第i行
{
	matrixTerm<T>& row_num = terms[row_first_position[i]];
	for (int j = 0; j < m.size; j++)  //对于第i行的各个元素，都遍历一遍B矩阵
	{
		matrixTerm<T>& col_num = m.terms[j];
		if (row_num.row == i && row_num.col == col_num.row)  //确保还在当前行
			accumulator[col_num.col] += row_num.value * col_num.value;
	}
	row_first_position[i]++;
}

为了方便，用row_num表示当前第i行使用的元素，col_num表示当前第j列使用的元素。外层循环的判断用到了row_nums数组，第i行有x个非零元素，就执行x次外层循环。if语句判断row_num还在第i行，且行元素和列元素满足累加条件。内层循环完之后，把row_first_position[i]++，表示三元组表的索引加1

for (int k = 1; k <= m.cols; k++)
{
	if (accumulator[k] != 0)   //矩阵C的第i行第k列有值
	{
		result.size++;
		result.terms[cursor_C].row = i;
		result.terms[cursor_C].col = k;
		result.terms[cursor_C].value = accumulator[k];
		cursor_C++;
	}
}

紧接着，把第i行1~m.cols号累加器中的非零值按照行主映射存到结果矩阵中去。存完后，清空累加器，继续进入下一次循环。
这个算法需要嵌套3层for循环：

注意它们之间的联系！
完整函数代码附下：

template <class T>
void sparseMatrix<T>::multiply(const sparseMatrix<T>& m)
{
	sparseMatrix<T> result;  //结果矩阵
	//使结果矩阵的行数等于A矩阵的行数，列数等于B矩阵的列数
	result.rows = rows;
	result.cols = m.cols;
	int cursor_C = 0;    //使用cursor_C来标记结果矩阵的三元组表
	result.terms = new matrixTerm<T>[rows * m.cols];  //按照最大情况来初始化
	result.size = 0;

	//用一个数组来表示A每一行第一个元素在三元组表中对应的索引
	int* row_nums = new int[rows + 1];   //各行的元素数目
	int* row_first_position = new int[rows + 1];  //各行第一个非零元对应的索引

	for (int i = 0; i <= rows; i++)
		row_nums[i] = 0;   //初始化row_nums
	for (int i = 0; i < size; i++)
		row_nums[terms[i].row]++;     //赋值row_nums

	for (int i = 0; i <= rows; i++)   //如果第i行没有元素，则置为-1
		row_first_position[i] = -1;
	row_first_position[1] = 0;       //第一行第一个非零元一定在0号索引
	for (int i = 2; i <= rows; i++)   //遍历三元组表
		row_first_position[i] = row_first_position[i - 1] + row_nums[i - 1];
	
	T* accumulator = new T[m.cols + 1];      //B有m.cols列，所以定义m.cols大小的累加器
	for (int i = 0; i <= m.cols; i++)
		accumulator[i] = 0;          //初始化累加器,分别表示各个列

	for (int i = 1; i <= rows; i++)    //逻辑上遍历矩阵A的每一行
	{
		if (row_first_position[i] == -1)   //矩阵的第i行为空
		{
			//do nothing
		}
		else   //矩阵的第i行非空，此时要遍历B的三元表，并在累加器中存值
		{
			for (int k = 1; k <= row_nums[i]; k++) //逻辑上遍历A的第i行
			{
				matrixTerm<T>& row_num = terms[row_first_position[i]];
				for (int j = 0; j < m.size; j++)  //对于第i行的各个元素，都遍历一遍B矩阵
				{
					matrixTerm<T>& col_num = m.terms[j];
					if (row_num.row == i && row_num.col == col_num.row)  //确保还在当前行
						accumulator[col_num.col] += row_num.value * col_num.value;
				}
				row_first_position[i]++;
			}
			//如此便得到了C矩阵第i行的各个元素的值，分别存在累加器中
			for (int k = 1; k <= m.cols; k++)
			{
				if (accumulator[k] != 0)   //矩阵C的第i行第k列有值
				{
					result.size++;
					result.terms[cursor_C].row = i;
					result.terms[cursor_C].col = k;
					result.terms[cursor_C].value = accumulator[k];
					cursor_C++;
				}
			}
			//如此便运算完了第i行并成功存值，接下来要归零累加器
			for (int k = 1; k <= m.cols; k++)
				accumulator[k] = 0;          //初始化累加器,分别表示各个列
		}
	}
	*this = result;
}