-
Notifications
You must be signed in to change notification settings - Fork 315
/
dpbsl.f
83 lines (81 loc) · 2.1 KB
/
dpbsl.f
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
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
subroutine dpbsl(abd,lda,n,m,b)
integer lda,n,m
double precision abd(lda,n),b(n)
c
c dpbsl solves the double precision symmetric positive definite
c band system a*x = b
c using the factors computed by dpbco or dpbfa.
c
c on entry
c
c abd double precision(lda, n)
c the output from dpbco or dpbfa.
c
c lda integer
c the leading dimension of the array abd .
c
c n integer
c the order of the matrix a .
c
c m integer
c the number of diagonals above the main diagonal.
c
c b double precision(n)
c the right hand side vector.
c
c on return
c
c b the solution vector x .
c
c error condition
c
c a division by zero will occur if the input factor contains
c a zero on the diagonal. technically this indicates
c singularity but it is usually caused by improper subroutine
c arguments. it will not occur if the subroutines are called
c correctly and info .eq. 0 .
c
c to compute inverse(a) * c where c is a matrix
c with p columns
c call dpbco(abd,lda,n,rcond,z,info)
c if (rcond is too small .or. info .ne. 0) go to ...
c do 10 j = 1, p
c call dpbsl(abd,lda,n,c(1,j))
c 10 continue
c
c linpack. this version dated 08/14/78 .
c cleve moler, university of new mexico, argonne national lab.
c
c subroutines and functions
c
c blas daxpy,ddot
c fortran min0
c
c internal variables
c
double precision ddot,t
integer k,kb,la,lb,lm
c
c solve trans(r)*y = b
c
do 10 k = 1, n
lm = min0(k-1,m)
la = m + 1 - lm
lb = k - lm
t = ddot(lm,abd(la,k),1,b(lb),1)
b(k) = (b(k) - t)/abd(m+1,k)
10 continue
c
c solve r*x = y
c
do 20 kb = 1, n
k = n + 1 - kb
lm = min0(k-1,m)
la = m + 1 - lm
lb = k - lm
b(k) = b(k)/abd(m+1,k)
t = -b(k)
call daxpy(lm,t,abd(la,k),1,b(lb),1)
20 continue
return
end