# Contents of /trunk/linpack/dgesl.f

Revision 1 - (show annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (18 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