Extract and create sparse band and diagonal matrices
B = spdiags(A)
[B,d] = spdiags(A)
B = spdiags(A,d)
A = spdiags(B,d,A)
A = spdiags(B,d,m,n)
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 |
The spdiags
function deals with three matrices,
in various combinations, as both input and output.
A | An |
B | A |
d | A vector of length |
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.
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.
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.
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)
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
.
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
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)
.
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