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

Annotation of /trunk/linpack/dgesl.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (hide annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (15 years, 9 months ago) by aw0a
File size: 3339 byte(s)
Setting up web subdirectory in repository
1 aw0a 1 C dgesl.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 dgesl(a,lda,n,ipvt,b,job)
7     integer lda,n,ipvt(1),job
8     double precision a(lda,1),b(1)
9     c
10     c dgesl solves the double precision system
11     c a * x = b or trans(a) * x = b
12     c using the factors computed by dgeco or dgefa.
13     c
14     c on entry
15     c
16     c a double precision(lda, n)
17     c the output from dgeco or dgefa.
18     c
19     c lda integer
20     c the leading dimension of the array a .
21     c
22     c n integer
23     c the order of the matrix a .
24     c
25     c ipvt integer(n)
26     c the pivot vector from dgeco or dgefa.
27     c
28     c b double precision(n)
29     c the right hand side vector.
30     c
31     c job integer
32     c = 0 to solve a*x = b ,
33     c = nonzero to solve trans(a)*x = b where
34     c trans(a) is the transpose.
35     c
36     c on return
37     c
38     c b the solution vector x .
39     c
40     c error condition
41     c
42     c a division by zero will occur if the input factor contains a
43     c zero on the diagonal. technically this indicates singularity
44     c but it is often caused by improper arguments or improper
45     c setting of lda . it will not occur if the subroutines are
46     c called correctly and if dgeco has set rcond .gt. 0.0
47     c or dgefa has set info .eq. 0 .
48     c
49     c to compute inverse(a) * c where c is a matrix
50     c with p columns
51     c call dgeco(a,lda,n,ipvt,rcond,z)
52     c if (rcond is too small) go to ...
53     c do 10 j = 1, p
54     c call dgesl(a,lda,n,ipvt,c(1,j),0)
55     c 10 continue
56     c
57     c linpack. this version dated 08/14/78 .
58     c cleve moler, university of new mexico, argonne national lab.
59     c
60     c subroutines and functions
61     c
62     c blas daxpy,ddot
63     c
64     c internal variables
65     c
66     double precision ddot,t
67     integer k,kb,l,nm1
68     c
69     nm1 = n - 1
70     if (job .ne. 0) go to 50
71     c
72     c job = 0 , solve a * x = b
73     c first solve l*y = b
74     c
75     if (nm1 .lt. 1) go to 30
76     do 20 k = 1, nm1
77     l = ipvt(k)
78     t = b(l)
79     if (l .eq. k) go to 10
80     b(l) = b(k)
81     b(k) = t
82     10 continue
83     call daxpy(n-k,t,a(k+1,k),1,b(k+1),1)
84     20 continue
85     30 continue
86     c
87     c now solve u*x = y
88     c
89     do 40 kb = 1, n
90     k = n + 1 - kb
91     b(k) = b(k)/a(k,k)
92     t = -b(k)
93     call daxpy(k-1,t,a(1,k),1,b(1),1)
94     40 continue
95     go to 100
96     50 continue
97     c
98     c job = nonzero, solve trans(a) * x = b
99     c first solve trans(u)*y = b
100     c
101     do 60 k = 1, n
102     t = ddot(k-1,a(1,k),1,b(1),1)
103     b(k) = (b(k) - t)/a(k,k)
104     60 continue
105     c
106     c now solve trans(l)*x = y
107     c
108     if (nm1 .lt. 1) go to 90
109     do 80 kb = 1, nm1
110     k = n - kb
111     b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1)
112     l = ipvt(k)
113     if (l .eq. k) go to 70
114     t = b(l)
115     b(l) = b(k)
116     b(k) = t
117     70 continue
118     80 continue
119     90 continue
120     100 continue
121     return
122     end

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