行索引稀疏矩阵乘法（严蔚敏版）

原理

行索引稀疏矩阵查找某一列的元素没那么方便，所以在做矩阵乘法时（这里以M乘N=Q为例），严书的做法是：在求Q的某一行的值是，用M的一行去遍历N的每一行，对结果中同列的值进行累加，最后稀疏化存入Q的当前行中，这个过程还原成正常矩阵比较容易理解：
求Q（2,2）的第一行时，肯定是M的第一行和N的第一列逐乘再累加，然后再M的第一行和N的第二列逐乘累加
M（2,3）* N（3,2）：

直接算的话就是:
Q[1,1] = 1*1 + 2*0 + 3*1 = 1+0+3 = 4
Q[1,2] = 1*0 + 2*1 + 3*1 = 0+2+3 = 5
这回我们直接同时求两列的值，其实就是改变一下计算顺序：

这个就相当于严书代码中ctemp的作用,只不过ctemp是直接累加.

代码

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;
const int MAX_E = 1e+4+5, MAX_R = 105;

//严版特色,下标从1开始
struct Triple{
    int i,j,e;};
struct RLSMatrix{
    
    Triple data[MAX_E];
    int rpos[MAX_R];//行首索引, rpos[i]指向第i行中的首元素在data中的序号, 则rpos[i+1]-1指向第i行中最后一个非0元素在N.data中的序号.
    int rows, cols, elems;
};

int ctemp[MAX_R];//Q的第i行元素结果累加器,算完一行就压缩存储进Q.data中


RLSMatrix mul(RLSMatrix M, RLSMatrix N){
    
    RLSMatrix Q;
    Q.rows = M.rows;Q.cols = N.cols;Q.elems = 0;
    if(M.elems * N.elems != 0){
                                  // Q是非零矩阵
        for(int arow = 1; arow <= M.rows; arow++){
    
            memset(ctemp, 0, sizeof(ctemp));              //到了下一行,清零
            Q.rpos[arow] = Q.elems+1;                  //新一行的行首索引是当前data中元素个数+1,从该处继续向data中存三元组

            int tp;      //用tp指向data中该行行末的下一个位置,方便遍历
            if(arow < M.rows) tp = M.rpos[arow+1];    //如果不是最后一行, tp指向下一行行首
            else tp = M.elems+1;   //最后一行, tp直接指向data中最后一个+1

            for(int p = M.rpos[arow];p < tp; p++){
      //对当前行找到的每一个非零元
                int brow = M.data[p].j;     //找到对应元在N中的行号(M中某行第三列的元素, 肯定和N中第三行某列的元素相乘)
                int t;//和tp同样的套路
                if(brow < N.rows) t = N.rpos[brow+1];
                else t = N.elems+1;

                for(int q = N.rpos[brow]; q < t; q++){
    //这里不是直接算出M的某行乘N的某列的值
                    //而是用M的某行,去遍历N的所有行,然后对每行的应该在同列的值进行累加
                    int ccol = N.data[q].j;
                    ctemp[ccol] += M.data[p].e * N.data[q].e;//累加器
                }
            }//求得Q中第crow(=arow)行的非零元

            for(int ccol = 1; ccol <= Q.cols; ++ccol){
    //用M的一行遍历完了整个N矩阵
                if(ctemp[ccol]){
    
                    Q.data[++Q.elems] = {
    arow, ccol, ctemp[ccol]};//压缩存储
                }
            }
        }
    }
    return Q;
}

RLSMatrix makeMat(){
    
    /*
     * 输入rows, cols, 三元组个数
     * 输入三元组
     */
    RLSMatrix A;
    cin>>A.rows>>A.cols>>A.elems;
    for(int i = 1;i <= A.elems;i++){
    
        int x,y,e;
        cin>>x>>y>>e;
        A.data[i] = {
    x,y,e};
    }

    /*
     根据乘法中这段代码写的初始化rpos数组:
         int tp;      //用tp指向data中该行行末的下一个位置,方便遍历
         if(arow < M.rows) tp = M.rpos[arow+1];    //如果不是最后一行, tp指向下一行行首
         else tp = M.elems+1;   //最后一行, tp直接指向data中最后一个+1
         for(int p = M.rpos[arow];p < tp; p++) ...
        可以看出如果某行没有元素,则让tp = M.rpos[arow]则可不进行遍历,也就是M.rpos[arow] = M.rpos[arow+1]
        而当最后几行为空时,并不会执行else tp = M.elems+1;这时应该把rpos最后几行空的填成M.elems+1
     */

    int row = 1;
    for(int e = 1;e <= A.elems; e++){
    
        int arow = A.data[e].i;
        while(row <= arow){
    
            A.rpos[row++] = e;
        }
    }
    while(row <= A.rows){
    //如果最后几行没有元素,让索引指到elems+1的位置
        A.rpos[row++] = A.elems+1;
    }
    return A;
}

void printMat(RLSMatrix A){
    
    cout<<"rows:"<<A.rows<<" cols:"<<A.cols<<" elems:"<<A.elems<<endl;
    cout<<"-------------------------"<<endl;
    cout<<"data:"<<endl;
    for(int i = 1;i <= A.elems;i++) {
    
        cout<<setw(4)<<A.data[i].i<<setw(4)<<A.data[i].j<<setw(4)<<A.data[i].e<<endl;
    }
    cout<<"-------------------------"<<endl;
    cout<<"rpos: ";
    for(int j = 1;j <= A.rows; j++){
    
        cout<<A.rpos[j]<<' ';
    }
    cout<<endl<<endl;
}

int main(){
    
    ios::sync_with_stdio(false);
    RLSMatrix A = makeMat();
    printMat(A);
    RLSMatrix B = makeMat();
    printMat(B);
    RLSMatrix C = mul(A, B);
    printMat(C);
    return 0;
}
/*
 严书上的例子
 3 4 4
 1 1 3
 1 4 5
 2 2 -1
 3 1 2

 4 2 4
 1 2 2
 2 1 1
 3 1 -2
 3 2 4
 */