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

Contents of /trunk/linpack/dgesl.f

Parent Directory Parent Directory | Revision Log Revision Log


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