| 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 |
Parent Directory
| 
Revision Log