1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
|
subroutine solbt (m, n, a, b, c, y, ip)
integer m, n, ip(m,n)
double precision a(m,m,n), b(m,m,n), c(m,m,n), y(m,n)
clll. optimize
c-----------------------------------------------------------------------
c solution of block-tridiagonal linear system.
c coefficient matrix must have been previously processed by decbt.
c m, n, a, b, c, and ip must not have been changed since call to decbt.
c written by a. c. hindmarsh.
c input..
c m = order of each block.
c n = number of blocks in each direction of matrix.
c a,b,c = m by m by n arrays containing block lu decomposition
c of coefficient matrix from decbt.
c ip = m by n integer array of pivot information from decbt.
c y = array of length m*n containg the right-hand side vector
c (treated as an m by n array here).
c output..
c y = solution vector, of length m*n.
c
c external routines required.. dgesl (linpack) and ddot (blas).
c-----------------------------------------------------------------------
c
integer nm1, nm2, i, k, kb, km1, kp1
double precision dp, ddot
nm1 = n - 1
nm2 = n - 2
c forward solution sweep. ----------------------------------------------
call dgesl (a, m, m, ip, y, 0)
do 30 k = 2,nm1
km1 = k - 1
do 20 i = 1,m
dp = ddot (m, c(i,1,k), m, y(1,km1), 1)
y(i,k) = y(i,k) - dp
20 continue
call dgesl (a(1,1,k), m, m, ip(1,k), y(1,k), 0)
30 continue
do 50 i = 1,m
dp = ddot (m, c(i,1,n), m, y(1,nm1), 1)
1 + ddot (m, b(i,1,n), m, y(1,nm2), 1)
y(i,n) = y(i,n) - dp
50 continue
call dgesl (a(1,1,n), m, m, ip(1,n), y(1,n), 0)
c backward solution sweep. ---------------------------------------------
do 80 kb = 1,nm1
k = n - kb
kp1 = k + 1
do 70 i = 1,m
dp = ddot (m, b(i,1,k), m, y(1,kp1), 1)
y(i,k) = y(i,k) - dp
70 continue
80 continue
do 100 i = 1,m
dp = ddot (m, c(i,1,1), m, y(1,3), 1)
y(i,1) = y(i,1) - dp
100 continue
return
c----------------------- end of subroutine solbt ---------------------
end
|
