/[ascend]/trunk/linpack/dgbsl.f
ViewVC logotype

Contents of /trunk/linpack/dgbsl.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (20 years ago) by aw0a
File size: 3972 byte(s)
Setting up web subdirectory in repository
1 C dgbsl.f
2 C is freely available from netlib. It is not subject to any GNU License
3 C set by the authors of the ASCEND math programming system.
4 C $Date: 1996/04/30 18:17:11 $ $Revision: 1.1.1.1 $
5 C
6 subroutine dgbsl(abd,lda,n,ml,mu,ipvt,b,job)
7 integer lda,n,ml,mu,ipvt(1),job
8 double precision abd(lda,1),b(1)
9 c
10 c dgbsl solves the double precision band system
11 c a * x = b or trans(a) * x = b
12 c using the factors computed by dgbco or dgbfa.
13 c
14 c on entry
15 c
16 c abd double precision(lda, n)
17 c the output from dgbco or dgbfa.
18 c
19 c lda integer
20 c the leading dimension of the array abd .
21 c
22 c n integer
23 c the order of the original matrix.
24 c
25 c ml integer
26 c number of diagonals below the main diagonal.
27 c
28 c mu integer
29 c number of diagonals above the main diagonal.
30 c
31 c ipvt integer(n)
32 c the pivot vector from dgbco or dgbfa.
33 c
34 c b double precision(n)
35 c the right hand side vector.
36 c
37 c job integer
38 c = 0 to solve a*x = b ,
39 c = nonzero to solve trans(a)*x = b , where
40 c trans(a) is the transpose.
41 c
42 c on return
43 c
44 c b the solution vector x .
45 c
46 c error condition
47 c
48 c a division by zero will occur if the input factor contains a
49 c zero on the diagonal. technically this indicates singularity
50 c but it is often caused by improper arguments or improper
51 c setting of lda . it will not occur if the subroutines are
52 c called correctly and if dgbco has set rcond .gt. 0.0
53 c or dgbfa has set info .eq. 0 .
54 c
55 c to compute inverse(a) * c where c is a matrix
56 c with p columns
57 c call dgbco(abd,lda,n,ml,mu,ipvt,rcond,z)
58 c if (rcond is too small) go to ...
59 c do 10 j = 1, p
60 c call dgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0)
61 c 10 continue
62 c
63 c linpack. this version dated 08/14/78 .
64 c cleve moler, university of new mexico, argonne national lab.
65 c
66 c subroutines and functions
67 c
68 c blas daxpy,ddot
69 c fortran min0
70 c
71 c internal variables
72 c
73 double precision ddot,t
74 integer k,kb,l,la,lb,lm,m,nm1
75 c
76 m = mu + ml + 1
77 nm1 = n - 1
78 if (job .ne. 0) go to 50
79 c
80 c job = 0 , solve a * x = b
81 c first solve l*y = b
82 c
83 if (ml .eq. 0) go to 30
84 if (nm1 .lt. 1) go to 30
85 do 20 k = 1, nm1
86 lm = min0(ml,n-k)
87 l = ipvt(k)
88 t = b(l)
89 if (l .eq. k) go to 10
90 b(l) = b(k)
91 b(k) = t
92 10 continue
93 call daxpy(lm,t,abd(m+1,k),1,b(k+1),1)
94 20 continue
95 30 continue
96 c
97 c now solve u*x = y
98 c
99 do 40 kb = 1, n
100 k = n + 1 - kb
101 b(k) = b(k)/abd(m,k)
102 lm = min0(k,m) - 1
103 la = m - lm
104 lb = k - lm
105 t = -b(k)
106 call daxpy(lm,t,abd(la,k),1,b(lb),1)
107 40 continue
108 go to 100
109 50 continue
110 c
111 c job = nonzero, solve trans(a) * x = b
112 c first solve trans(u)*y = b
113 c
114 do 60 k = 1, n
115 lm = min0(k,m) - 1
116 la = m - lm
117 lb = k - lm
118 t = ddot(lm,abd(la,k),1,b(lb),1)
119 b(k) = (b(k) - t)/abd(m,k)
120 60 continue
121 c
122 c now solve trans(l)*x = y
123 c
124 if (ml .eq. 0) go to 90
125 if (nm1 .lt. 1) go to 90
126 do 80 kb = 1, nm1
127 k = n - kb
128 lm = min0(ml,n-k)
129 b(k) = b(k) + ddot(lm,abd(m+1,k),1,b(k+1),1)
130 l = ipvt(k)
131 if (l .eq. k) go to 70
132 t = b(l)
133 b(l) = b(k)
134 b(k) = t
135 70 continue
136 80 continue
137 90 continue
138 100 continue
139 return
140 end

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22