spdiags (MATLAB Functions)

Extract and create sparse band and diagonal matrices

Syntax

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

Description

The spdiags function generalizes the function diag. Four different operations, distinguished by the number of input arguments, are possible.

B = spdiags(A) extracts all nonzero diagonals from the m-by-n matrix A. B is a min(m,n)-by-p matrix whose columns are the p nonzero diagonals of A.

[B,d] = spdiags(A) returns a vector d of length p, whose integer components specify the diagonals in A.

B = spdiags(A,d) extracts the diagonals specified by d.

A = spdiags(B,d,A) replaces the diagonals specified by d with the columns of B. The output is sparse.

A = spdiags(B,d,m,n) creates an m-by-n sparse matrix by taking the columns of B and placing them along the diagonals specified by d.

Note In this syntax, if a column of B is longer than the diagonal it is replacing, and m >= n, spdiags takes elements of super-diagonals from the lower part of the column of B, and elements of sub-diagonals from the upper part of the column of B. However, if m < n, then super-diagonals are from the upper part of the column of B, and sub-diagonals from the lower part. (See Example 5, below).

Arguments

The spdiags function deals with three matrices, in various combinations, as both input and output.

A
An m-by-n matrix, usually (but not necessarily) sparse, with its nonzero or specified elements located on p diagonals.

B
A min(m,n)-by-p matrix, usually (but not necessarily) full, whose columns are the diagonals of A.

d
A vector of length p whose integer components specify the diagonals in A.

`A`	An `m`-by-`n` matrix, usually (but not necessarily) sparse, with its nonzero or specified elements located on `p` diagonals.
`B`	A `min(m,n)`-by-`p` matrix, usually (but not necessarily) full, whose columns are the diagonals of `A`.
`d`	A vector of length `p` whose integer components specify the diagonals in `A`.

Roughly, A, B, and d are related by

for k = 1:p
    B(:,k) = diag(A,d(k))
end

Some elements of B, corresponding to positions outside of A, are not defined by these loops. They are not referenced when B is input and are set to zero when B is output.

How the Diagonals of A are Listed in the Vector d

An m-by-n matrix A has m+n-1diagonals. These are specified in the vector d using indices from -m+1 to n-1. For example, if A is 5-by-6, it has 10 diagonals, which are specified in the vector d using the indices -4, -3 , ... 4, 5. The following diagram illustrates this for a vector of all ones.

Examples

Example 1

For the following matrix,

A=[0 5 0 10 0 0;...
0 0 6 0 11 0;...
3 0 0 7 0 12;...
1 4 0 0 8 0;...
0 2 5 0 0 9]

A =

     0     5     0    10     0     0
     0     0     6     0    11     0
     3     0     0     7     0    12
     1     4     0     0     8     0
     0     2     5     0     0     9

the command

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

returns

B =

     0     0     5    10
     0     0     6    11
     0     3     7    12
     1     4     8     0
     2     5     9     0


d =


    -3
    -2
     1
     3

The columns of the first output B contain the nonzero diagonals of A. The second output d lists the indices of the nonzero diagonals of A, as shown in the following diagram. See How the Diagonals of A are Listed in the Vector d.

Note that the longest nonzero diagonal in A is contained in column 3 of B. The other nonzero diagonals of A have extra zeros added to their corresponding columns in B, to give all columns of B the same length. For the nonzero diagonals below the main diagonal of A, extra zeros are added at the tops of columns. For the nonzero diagonals above the main diagonal of A, extra zeros are added at the bottoms of columns. This is illustrated by the following diagram.

Example 2

This example generates a sparse tridiagonal representation of the classic second difference operator on n points.

e = ones(n,1);
A = spdiags([e -2*e e], -1:1, n, n)

Turn it into Wilkinson's test matrix (see gallery):

A = spdiags(abs(-(n-1)/2:(n-1)/2)',0,A)

Finally, recover the three diagonals:

```
B = spdiags(A)
```

Example 3

The second example is not square.

A = [11    0   13    0
      0   22    0   24
      0    0   33    0
     41    0    0   44
      0   52    0    0
      0    0   63    0
      0    0    0   74]

Here m = 7, n = 4, and p = 3.

The statement [B,d] = spdiags(A) produces d = [-3 0 2]' and

B = [41   11    0
     52   22    0
     63   33   13
     74   44   24]

Conversely, with the above B and d, the expression spdiags(B,d,7,4) reproduces the original A.

Example 4

This example shows how spdiags creates the diagonals when the columns of B are longer than the diagonals they are replacing.

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

B =

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

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

ans =

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

Example 5A

This example illustrates the use of the syntax A = spdiags(B,d,m,n), under three conditions:

m is equal to n
m is greater than n
m is less than n

The command used in this example is

```
A = full(spdiags(B, [-2 0 2], m, n))
```

where B is the 5-by-3 matrix shown below. The resulting matrix A has dimensions m-by-n, and has nonzero diagonals at [-2 0 2] (a sub-diagonal at -2, the main diagonal, and a super-diagonal at 2).

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

The first and third columns of matrix B are used to create the sub- and super-diagonals of A respectively. In all three cases though, these two outer columns of B are longer than the resulting diagonals of A. Because of this, only a part of the columns is used in A.

When m == n or m > n, spdiags takes elements of the super-diagonal in A from the lower part of the corresponding column of B, and elements of the sub-diagonal in A from the upper part of the corresponding column of B.

When m < n, spdiags does the opposite, taking elements of the super-diagonal in A from the upper part of the corresponding column of B, and elements of the sub-diagonal in A from the lower part of the corresponding column of B.

Part 1 -- m is equal to n

A = full(spdiags(B, [-2 0 2], 5, 5))
 Matrix B                      Matrix A

1    6   11                 6    0   13    0    0
2    7   12                 0    7    0   14    0
3    8   13  == spdiags =>  1    0    8    0   15
4    9   14                 0    2    0    9    0
5   10   15                 0    0    3    0   10

A(3,1), A(4,2), and A(5,3) are taken from the upper part of B(:,1).
A(1,3), A(2,4), and A(3,5) are taken from the lower part of B(:,3).

Part 2 -- m is greater than n

A = full(spdiags(B, [-2 0 2], 5, 4))
 Matrix B                      Matrix A

1    6   11                 6    0   13    0
2    7   12                 0    7    0   14
3    8   13  == spdiags =>  1    0    8    0
4    9   14                 0    2    0    9
5   10   15                 0    0    3    0

Same as in Part A.

Part 3 -- m is less than n

A = full(spdiags(B, [-2 0 2], 4, 5))
 Matrix B                      Matrix A

1    6   11                 6    0   11    0    0
2    7   12                 0    7    0   12    0
3    8   13  == spdiags =>  3    0    8    0   13
4    9   14                 0    4    0    9    0
5   10   15

A(3,1) and A(4,2) are taken from the lower part of B(:,1).
A(1,3), A(2,4), and A(3,5) are taken from the upper part of B(:,3).

Example 5B

Extract the diagonals from the first part of this example back into a column format using the command

```
B = spdiags(A)
```

You can see that in each case the original columns are restored (minus those elements that had overflowed the super- and sub-diagonals of matrix A).

Part 1

      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

Part 2

    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

Part 3

      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

See Also

diag

spconvert speye