数据结构——矩阵压缩与压缩矩阵的转置与乘法计算

　　为什么要对矩阵进行压缩存储呢？对于一个n*m的矩阵，我们一般会想到开一个n*m的二维数组来存储，这样计算操作都很方便模拟，但当一个矩阵很大时，这样对于空间的开销与浪费是很可怕的，尤其是当矩阵变成多维时。但我们往往不会在矩阵每一个位置都存有数据，很多矩阵元素其实是0，我们需要记录的只是那些非零元素，于是我们可以记录非零元素的位置与值，这样便可以大大减少空间上的浪费。

矩阵压缩存储代码（注意要按行顺序遍历存储矩阵非零元素）：

typedef struct {
    int i, j;           //非零元行下标与列下标
    ElemType e;
}Triple;

typedef struct {
    Triple data[MAXSIZE + 1];
    int mu, nu, tu;     //行数， 列数， 非零元个数
}TSMatrix;

　　矩阵转置是矩阵运算中极为重要的操作，但矩阵压缩后我们不能再按原来的方式用n*m的遍历方式进行转置操作，那我们应该怎么进行转置呢，首先想到的是以列为参考值对元素进行遍历转置，这样访问元素的顺序正好是转制后存储的顺序，代码如下：

Status TransposeSMatrix(TSMatrix M, TSMatrix & T) {
    T.mu = M.nu;
    T.nu = M.mu;
    T.tu = M.tu;
    if(T.tu) {
        int q = 1;
        for(int col = 1; col <= M.nu; ++col) 
            for(int p = 1; p <= M.tu; ++p)
                if(M.data[p].j == col) {
                    T.data[q].i = M.data[p].j;
                    T.data[q].j = M.data[p].i;
                    T.data[q].e = M.data[p].e;
                    ++q;
                }
    }
    return OK;
}

　　我们注意到当矩阵不那么稀疏时，比如如果这是一个稠密矩阵甚至n*m的空间里每一个位置都是一个非零元素，这时这个代码的复杂度达到了O(mu*nu*nu)，比直接用数组存储的O(mu*nu)还要大一个数量级，这样节省一点空间浪费大量时间的做法显然是不合适的，我们考虑改进算法，我们注意到如果我们预先知道矩阵M每一列的第一个非零元素在T中应有的位置，那么转置中我们就可以直接把他放在那个位置，这时，我们考虑开辅助数组记录每一列第一个非零元素的位置，以及每一列非零元素的数量，但其实，每一列的第一个非零元素的位置就是前一列第一个非零元素位置加上前一列非零元素数量，我们每次维护剩余元素在剩余矩阵中每一列第一个非零元素的位置数组，就得到了下面的O(nu + tu)代码:

Status FastTransposeSMtrix(TSMatrix M, TSMatrix &T) {
    T.mu = M.nu;
    T.nu = M.mu;
    T.tu = M.tu;
    if(T.tu) {
        int *num = new int[M.nu + 1], *cpot = new int[M.tu + 1];
        for(int col = 1; col <= M.nu; ++ col)
            num[col] = 0;
        for(int t = 1; t <= M.tu; ++t) ++num[M.data[t].j];
        cpot[1] = 1;
        for(int col = 2; col <= M.nu; ++ col)
            cpot[col] = cpot[col - 1] + num[col - 1];
        for(int p = 1; p <= M.tu; ++p) {
            int col = M.data[p].j, q = cpot[col];
            T.data[q].i = M.data[p].j;
            T.data[q].j = M.data[p].i;
            T.data[q].e = M.data[p].e;
            ++cpot[col];
        }
        delete num;
        delete cpot;
    }
    return OK;
}

矩阵转置到此就结束了，写进一个cpp效果就像下面这样，感兴趣的可以自己测试下效果：

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#include <cstring>
#include <vector>
#include <string>
#include <queue>
#include <map>
#include <set>

#define FRER() freopen("in.txt", "r", stdin);
#define INF 0x3f3f3f3f

using namespace std;

//函数状态码定义
#define TRUE        1
#define FALSE       0
#define OK          1
#define ERROR       0
#define INFEASIBLE -1
//#define OVERFLOW   -2

typedef int Status;
typedef int ElemType;

#define MAXSIZE 12500

typedef struct {
    int i, j;           //非零元行下标与列下标
    ElemType e;
}Triple;

typedef struct {
    Triple data[MAXSIZE + 1];
    int mu, nu, tu;     //行数， 列数， 非零元个数
}TSMatrix;

Status TransposeSMatrix(TSMatrix M, TSMatrix & T) {
    T.mu = M.nu;
    T.nu = M.mu;
    T.tu = M.tu;
    if(T.tu) {
        int q = 1;
        for(int col = 1; col <= M.nu; ++col) 
            for(int p = 1; p <= M.tu; ++p)
                if(M.data[p].j == col) {
                    T.data[q].i = M.data[p].j;
                    T.data[q].j = M.data[p].i;
                    T.data[q].e = M.data[p].e;
                    ++q;
                }
    }
    return OK;
}

Status FastTransposeSMtrix(TSMatrix M, TSMatrix &T) {
    T.mu = M.nu;
    T.nu = M.mu;
    T.tu = M.tu;
    if(T.tu) {
        int *num = new int[M.nu + 1], *cpot = new int[M.tu + 1];
        for(int col = 1; col <= M.nu; ++ col)
            num[col] = 0;
        for(int t = 1; t <= M.tu; ++t) ++num[M.data[t].j];
        cpot[1] = 1;
        for(int col = 2; col <= M.nu; ++ col)
            cpot[col] = cpot[col - 1] + num[col - 1];
        for(int p = 1; p <= M.tu; ++p) {
            int col = M.data[p].j, q = cpot[col];
            T.data[q].i = M.data[p].j;
            T.data[q].j = M.data[p].i;
            T.data[q].e = M.data[p].e;
            ++cpot[col];
        }
        delete num;
        delete cpot;
    }
    return OK;
}

int main()
{
    //FRER()

    TSMatrix MatrixA, MatrixB;
    cin >> MatrixA.mu >> MatrixA.nu >> MatrixA.tu;
    for(int i = 1; i <= MatrixA.tu; ++i)
        cin >> MatrixA.data[i].i >> MatrixA.data[i].j >> MatrixA.data[i].e;
    //TransposeSMatrix(MatrixA, MatrixB);
    FastTransposeSMtrix(MatrixA, MatrixB);
    for(int i = 1; i <= MatrixB.tu; ++i)
        cout << MatrixB.data[i].i << ' ' << MatrixB.data[i].j << ' ' << MatrixB.data[i].e << endl;

    return 0;
}

/*测试数据
6 7 8
1 2 12
1 3 9
3 1 -3
3 6 14
4 3 24
5 2 18
6 1 15
6 4 -7
 */

 　　接下来是压缩矩阵的乘法，类似转置，我们直接在矩阵结构体里开一个数组用来存矩阵每一行第一个非零元素的在M中位置，就像这样：

typedef struct {
    Triple data[MAXSIZE + 1];   //非零元三元组表
    int rops[MAXSIZE + 1];      //各行第一个非零元位置
    int mu, nu, tu;             //行数，列数，非零元个数
}RLSMatrix;

　　由于稀疏矩阵相乘不一定还是稀疏矩阵，所以我们要根据结果判断元素是否是非零元，模拟乘法运算如下：

Status MultSMatrix(RLSMatrix M, RLSMatrix N, RLSMatrix & Q) {
    if(M.nu != N.mu) return ERROR;
    Q.mu = M.mu;
    Q.nu = N.nu;
    Q.tu = 0;
    int *ctemp = new int[Q.nu + 1];

    if(M.mu * N.nu != 0) {
        for(int arow = 1; arow <= M.mu; ++arow) {
            int tp;
            memset(ctemp, 0, (Q.nu + 1) * sizeof(int));
            Q.rops[arow] = Q.tu + 1;
            if(arow < M.mu)
                tp = M.rops[arow + 1];
            else
                tp = M.tu + 1;
            for(int p = M.rops[arow]; p < tp; ++p) {        //对当前行的每一个非零元进行计算
                int brow = M.data[p].j;                     //找到对应元在N中的标号
                int t;
                if(brow < N.mu) t = N.rops[brow + 1];
                else    t = N.tu + 1;
                for(int q = N.rops[brow]; q < t; ++q) {
                    int ccol = N.data[q].j;
                    ctemp[ccol] += M.data[p].e * N.data[q].e;
                }
            }
            for(int ccol = 1; ccol <= Q.nu; ++ccol)         //压缩存储该行非零元
                if(ctemp[ccol]) {
                    if(++Q.tu > MAXSIZE) return ERROR;
                    Q.data[Q.tu].i = arow;
                    Q.data[Q.tu].j = ccol;
                    Q.data[Q.tu].e = ctemp[ccol];
                }
        }
    }
    return OK;
}