10.8.3 spdiags稀疏對角指令

前文討論的對角元素處理指令已可隨意就斜角方向置入任意元素。實際上在解題時，尤其進行數值分析過程中，所形成之矩陣雖然很大，但百分之八、九十的空間大部份都是為零。對於資料處理方面是一種相當大的資源浪費。為解決是項問題，必須使用稀鬆矩陣(Sparse matrix)處理，以節省時間與空間。spdiags指令則是針對對角元素部份之存取運作，與前面討論之diag之功能可相輔相成。

稀疏對角處埋指令格式如下：

 B = spdiags(A)
 [B,d] = spdiags(A)
 B = spdiags(A,d)
 A = spdiags(B,d,A)
 A = spdiags(B,d,m,n)

spdiags(A)指令以不同的方式處理三個矩陣，可以作為輸入與輸出。其組成係利用不同之輸入參數完成各項功能。其中矩陣Ａ之大小可為mxn，只要陳述非零元素之位址及數值即可，Ｂ則為一輸出項，旨在描述對角線項之陣列，其大小為min(m,n) x p；d則表示矩陣Ａ中非零元素之之對角序列值。

矩陣Ａ中若大小為m x n，則應有m+n-1條對角線。其序列以d參數表示，其範圍為由-m+1至n-1。例如，Ａ之大小為5x6，則應有10條對角線，此10條線分別以d=-4,-3,-2,-1,0,1,2,3,4,5表示。d=0時表示正好在最長的對角線上。其位置可以由下圖表示。


例：

>> e=ones(4,1);A=diag(1:5)+diag(e,1)+diag(e,-1)

A =

    1     1     0     0     0
    1     2     1     0     0
    0     1     3     1     0
    0     0     1     4     1
    0     0     0     1     5

>> [B,d]=spdiags(A)

B =

    1     1     0
    1     2     1
    1     3     1
    1     4     1
    0     5     1


d =

   -1
    0
    1

由結果可知，d指示k值所在，而Ｂ值則是非零對角線上之值，其起點係以左邊項為基準，不足的部份補零。就本例而言第一行為 k=-1，故由基線往右下延伸，不足部分填零；第三行為 k=1，由左邊基線開始，有部份不足，須先補零。這些值可以利用下列指令取得Ｂ中之值列，例如：由矩陣Ａ中取k=1之對角線項值：

>> C=spdiags(A,1)

C =

    0
    1
    1
    1
    1

若要在現有之矩陣增加一對角線項，則可利用 spdiags(B,d,A) 這個指令函數，其中B為填入之值，必須以行的型式，所給的值若超過矩陣Ａ之範圍，則多餘的項會被截斷。d則指示置於那一條對角線上：

>> AA=spdiags([3 4 5]',-3,A)

AA =

  (1,1)        1
  (2,1)        1
  (4,1)        3
  (1,2)        1
  (2,2)        2
  (3,2)        1
  (5,2)        4
  (2,3)        1
  (3,3)        3
  (4,3)        1
  (3,4)        1
  (4,4)        4
  (5,4)        1
  (4,5)        1
  (5,5)        5

>> full(AA)

ans =

    1     1     0     0     0
    1     2     1     0     0
    0     1     3     1     0
    3     0     1     4     1
    0     4     0     1     5

上述指令之應用是在現有之矩陣填入對角線資料，下面的指令spdiags(B,d,m,n)則是創造一個新的稀疏矩陣，其大小為m-by-n。 注意陣列Ｂ與d之長度應相符合：

>> B = repmat((1:5)',[1 6])

B =

    1     1     1     1     1     1
    2     2     2     2     2     2
    3     3     3     3     3     3
    4     4     4     4     4     4
    5     5     5     5     5     5

>> d = [-4 -2 -1 0 3 4];
>> A = spdiags(B,d,5,5);
>> full(A)

ans =

    1     0     0     4     5
    1     2     0     0     5
    1     2     3     0     0
    0     2     3     4     0
    1     0     3     4     5

由所得之Ａ矩陣中可知，其大小受m x n之限制，而左下方之排例由B之高項而下，超過者將被截斷。右上方則截斷原需補足的部份，或其尾項由最後項填起。上述之說明實際上僅討論到m=n的情況，實際上當m與n不同時，有些情況因為截斷位置之不同，結果也不一樣。以下面為例：

>> B=reshape(1:15,5,3)

B =

    1     6    11
    2     7    12
    3     8    13
    4     9    14
    5    10    15

>> d=[-2 0 2]

d =

   -2     0     2

>> A=full(spdiags(B,d,5,5))

A =

    6     0    13     0     0
    0     7     0    14     0
    1     0     8     0    15
    0     2     0     9     0
    0     0     3     0    10

當m=n=5時，其填入之順序與前例相同，左下依右邊基準填入，再行截斷；右上則依底部填起，上部截斷。

>> A=full(spdiags(B,d,5,4))

A =

    6     0    13     0
    0     7     0    14
    1     0     8     0
    0     2     0     9
    0     0     3     0

當m＞n時，其左下依右邊基準填入，再行截斷，情況與前例相同；但右上因為截斷的關係，可依k=0處之截斷位置(此處為9倒數第二行)，則依底部填起，此時底部也要截斷一行，再依序往上，截斷超出者。

>> A=full(spdiags(B,d,4,5))

A =

    6     0    11     0     0
    0     7     0    12     0
    3     0     8     0    13
    0     4     0     9     0

當m＜n時，填入之順序相反過來，其左下改依底部往上填入，但依k=0處之截斷位置(此處為9倒數第二行)，則依底部填起，此時底部也要截斷一行，再依序往上，截斷超出者。右上邊則改由上面之基準填入，再依長度進行截斷。此種情況完全與前面之m=n及m＞n之情況相反，須特別注意。

在使用spdiags(A)，取出其對角線值置入Ｂ矩陣時，也有類似情況，當m＜n時，其取出之順序也會相反，例如：

當m=n時：

     Matrix A                        Matrix B

6    0   13    0    0                 1   6   0
0    7    0   14    0                 2   7   0
1    0    8    0   15  == spdiags =>  3   8  13
0    2    0    9    0                 0   9  14
0    0    3    0   10                 0  10  15

當m＞n時：

   Matrix A                     Matrix B

6    0   13    0                 1   6   0
0    7    0   14                 2   7   0
1    0    8    0  == spdiags =>  3   8  13
0    2    0    9                 0   9  14
0    0    3    0

當m＜n時：

     Matrix A                        Matrix B

6    0   11    0    0                 0   6  11
0    7    0   12    0                 0   7  12
3    0    8    0   13  == spdiags =>  3   8  13
0    4    0    9    0                 4   9   0