# Contents of /trunk/linpack/dgbsl.f

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