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

Annotation of /trunk/linpack/dgbsl.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (hide annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (17 years, 2 months ago) by aw0a
File size: 3972 byte(s)
Setting up web subdirectory in repository
1 aw0a 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