|
|
Contents | Previous | Next | Subchapters |
| Syntax |
[E, lam] = BlockTriDiag(B, C, R) |
[E, lam] = BlockTriDiag(B, C, R, level) | |
| See Also | tridiag |
{e(k)}
in the system of equations
/ b(1) c(2)' 0 ... 0 \ / e(1) \ / r(1) \
| c(2) b(2) c(3)' . | | e(2) | | r(2) |
| . . | | . | = | . |
| . | | . | | . |
\ 0 ... 0 c(N) b(N) / \ e(N) / \ r(N) /
b(k) is a symmetric positive definite
n by n matrix,
each c(k) is an
n by n matrix,
each e(k)
is an n by m matrix,
and each r(k)
is an n by m matrix.
{d(k)},
generated by the algorithm,
are assumed to be invertible.
This is guarantee if the block tri-diagonal matrix above is positive definite
(see Vassilevski's article).
n * N rows and n columns.
The value of b(k) is equal to the following
submatrix of B:
B.row((k-1) * n + 1, n)
n * N rows and n columns.
The value of c(k) is equal to the following
submatrix of C:
C.row((k-1) * n + 1, n)
Note that the first n rows of C are not used.
n * N rows and m columns.
The value of r(k) is equal to the following
submatrix of R:
R.row((k-1) * n + 1, m)
n * N rows and m columns.
The value of e(k) is equal to the following
submatrix of E:
E.row((k-1) * n + 1, m)
level < 0, no tracing is done.
This is the default value for level if it is not present.
The block tri-diagonal matrix above is formed in memory and referred
to a T below.
This can be a very large matrix and is not formed unless
level > 0.
| Level | Label | Description |
level > 1 |
type(T)
|
type of the matrix T |
level > 1 |
cond(T)
|
condition number of the matrix T |
level > 1 |
lam
|
block tri-diagonal calculation of logdet(T) |
level > 1 |
logdet(T)
|
direct calculation of logdet(T) |
level > 2 |
E
| block tri-diagonal calculation of solution of equation |
level > 2 |
T \ R
| direct calculation of solution of equation |
level > 3 |
T
| the large block tri-diagonal matrix |
clear
include BlockTriDiag.oms
#
function RandomPositive(n) begin
F = rand(n, n)
return F' * F
end
#
# dimensions of the problem
N = 4
n = 2
m = 1
#
# dimension values
T = fill(0d0, n * N, n * N)
Y = fill(0d0, n * N, m)
B = fill(0d0, n * N, n)
C = fill(0d0, n * N, n)
R = fill(0d0, n * N, m)
#
# Generate the sequences q(k), a(k), u(k), b(k), and c(k)
# as per the conditions in Lemma 7
#
# q(0) and a(0)
qkm = RandomPositive(n)
akm = rand(n, n)
for k = 1 to N begin
# generate random u(k), q(k), a(k), r(k)
uk = RandomPositive(n)
qk = RandomPositive(n)
ak = rand(n, n)
rk = rand(n, m)
#
bk = uk + inv(qkm) + ak * ( qk \ ak' )
ck = qkm \ akm'
#
# Fill in block tri-diagonal system of equations
kn = (k - 1) * n + 1
T.blk(kn, kn, n, n) = bk
if k > 1 then begin
T.blk(kn, kn - n, n, n) = ck
T.blk(kn - n, kn, n, n) = ck'
end
Y.blk(kn, 1, n, m) = rk
#
# store b(k), c(k), and r(k) in form for BlockTriDiag
kn = (k - 1) * n + 1
B.row(kn, n) = bk
C.row(kn, n) = ck
R.row(kn, n) = rk
end
#
# solve the block Tri-Diagonal System
level = 3
[E, lam] = BlockTriDiag(B, C, R, level)
Reference