數據結構——矩陣壓縮與壓縮矩陣的轉置與乘法計算

　　為什么要對矩陣進行壓縮存儲呢？對於一個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;
}