行逻辑链接顺序表的稀疏矩阵乘法

假设以下有两个矩阵

• M

  M
  3	0	0	5
  0	-1	0	0
  2	0	0	0

• N

  N
  0	2
  1	0
  -2	4
  0	0

建立一个与N的列同宽的数组temp缓存结果矩阵中当前行的结果

1. (通过行索引)遍历M的所有行，假设当前遍历到了M的第一行[ 3, 0, 0, 5 ]
2. 清空temp数组
3. 遍历当前行的所有列（我们现在讨论的是稀疏矩阵的乘法运算，所以是所有非零元）
   1. 先选中了第一列。那么便在矩阵N中找到第一行，遍历N的这一行，与3相乘后加到数组temp的对应位置上。例如这一次操作完成后，temp == [ 0, 6 ]
   2. 选中第四列，那么便在矩阵N中找到第四行，遍历N的这一行，与5相乘后加到数组temp的对应位置上。
4. 将temp数组压缩存储到结果（稀疏）矩阵中（即得到了结果矩阵这一行），清空temp数组，返回3遍历下一行。

typedef struct
{
    unsigned int i, j;
    int e;
} Element_t;
typedef struct
{
    unsigned int nR, nC, nU;
    Element_t data[16];
    //行偏移量，即各行第一个非零元的位置表
    unsigned int rowOffset[16];
} Matrix_t;

inline void matrixMultiply(
    const Matrix_t &left,
    const Matrix_t &right
)
{
    for (
        unsigned int currentRowIndex = 0; 
        currentRowIndex < left.nR; 
        currentRowIndex++
        )//遍历左矩阵的所有行
    {
        //建立临时数组来缓存结果行
        int *temp = reinterpret_cast<int *>(
            malloc(right.nC * sizeof(int)));
        memset(temp, 0, right.nC * sizeof(int));
        //遍历当前行的所有非零元素
        for (unsigned int currentUnitIndex = left.rowOffset[currentRowIndex];
            currentUnitIndex < (currentRowIndex < (left.nR-1) ? left.rowOffset[currentRowIndex + 1] : left.nU);
            ++currentUnitIndex)
        {
            //根据当前元素所在的列，找到右边矩阵对应的行进行遍历
            for (unsigned int rightUnitIndex = right.rowOffset[left.data[currentUnitIndex].j];
                rightUnitIndex < (left.data[currentUnitIndex].j < left.nC ? right.rowOffset[left.data[currentUnitIndex].j + 1] : right.nU);
                rightUnitIndex++)
            {
                temp[right.data[rightUnitIndex].j] += left.data[currentUnitIndex].e * right.data[rightUnitIndex].e;
            }
        }

        //输出当前行，其实应该是存入结果矩阵的
        for (unsigned int i = 0; i < right.nC; i++)
        {
            std::cout << *(temp + i) << "\t" << std::flush;
        }
        std::cout << std::endl << std::flush;
    }
}

int main()
{
    Matrix_t M, N;
    M.data[0].i = 0, M.data[0].j = 0, M.data[0].e = 3;
    M.data[1].i = 0, M.data[1].j = 3, M.data[1].e = 5;
    M.data[2].i = 1, M.data[2].j = 1, M.data[2].e = -1;
    M.data[3].i = 2, M.data[3].j = 0, M.data[3].e = 2;
    M.nC = 4;
    M.nR = 3;
    M.nU = 4;
    M.rowOffset[0] = 0, M.rowOffset[1] = 2, M.rowOffset[2] = 3;

    N.data[0].i = 0, N.data[0].j = 1, N.data[0].e = 2;
    N.data[1].i = 1, N.data[1].j = 0, N.data[1].e = 1;
    N.data[2].i = 2, N.data[2].j = 0, N.data[2].e = -2;
    N.data[3].i = 2, N.data[3].j = 1, N.data[3].e = 4;
    N.nC = 2;
    N.nR = 4;
    N.nU = 4;
    N.rowOffset[0] = 0, N.rowOffset[1] = 1, N.rowOffset[2] = 2;

    matrixMultiply(M, N);
}