| subroutine dtrco(t,ldt,n,rcond,z,job) |
| integer ldt,n,job |
| double precision t(ldt,n),z(n) |
| double precision rcond |
| c |
| c dtrco estimates the condition of a double precision triangular |
| c matrix. |
| c |
| c on entry |
| c |
| c t double precision(ldt,n) |
| c t contains the triangular matrix. the zero |
| c elements of the matrix are not referenced, and |
| c the corresponding elements of the array can be |
| c used to store other information. |
| c |
| c ldt integer |
| c ldt is the leading dimension of the array t. |
| c |
| c n integer |
| c n is the order of the system. |
| c |
| c job integer |
| c = 0 t is lower triangular. |
| c = nonzero t is upper triangular. |
| c |
| c on return |
| c |
| c rcond double precision |
| c an estimate of the reciprocal condition of t . |
| c for the system t*x = b , relative perturbations |
| c in t and b of size epsilon may cause |
| c relative perturbations in x of size epsilon/rcond . |
| c if rcond is so small that the logical expression |
| c 1.0 + rcond .eq. 1.0 |
| c is true, then t may be singular to working |
| c precision. in particular, rcond is zero if |
| c exact singularity is detected or the estimate |
| c underflows. |
| c |
| c z double precision(n) |
| c a work vector whose contents are usually unimportant. |
| c if t is close to a singular matrix, then z is |
| c an approximate null vector in the sense that |
| c norm(a*z) = rcond*norm(a)*norm(z) . |
| c |
| c linpack. this version dated 08/14/78 . |
| c cleve moler, university of new mexico, argonne national lab. |
| c |
| c subroutines and functions |
| c |
| c blas daxpy,dscal,dasum |
| c fortran dabs,dmax1,dsign |
| c |
| c internal variables |
| c |
| double precision w,wk,wkm,ek |
| double precision tnorm,ynorm,s,sm,dasum |
| integer i1,j,j1,j2,k,kk,l |
| logical lower |
| c |
| lower = job .eq. 0 |
| c |
| c compute 1-norm of t |
| c |
| tnorm = 0.0d0 |
| do 10 j = 1, n |
| l = j |
| if (lower) l = n + 1 - j |
| i1 = 1 |
| if (lower) i1 = j |
| tnorm = dmax1(tnorm,dasum(l,t(i1,j),1)) |
| 10 continue |
| c |
| c rcond = 1/(norm(t)*(estimate of norm(inverse(t)))) . |
| c estimate = norm(z)/norm(y) where t*z = y and trans(t)*y = e . |
| c trans(t) is the transpose of t . |
| c the components of e are chosen to cause maximum local |
| c growth in the elements of y . |
| c the vectors are frequently rescaled to avoid overflow. |
| c |
| c solve trans(t)*y = e |
| c |
| ek = 1.0d0 |
| do 20 j = 1, n |
| z(j) = 0.0d0 |
| 20 continue |
| do 100 kk = 1, n |
| k = kk |
| if (lower) k = n + 1 - kk |
| if (z(k) .ne. 0.0d0) ek = dsign(ek,-z(k)) |
| if (dabs(ek-z(k)) .le. dabs(t(k,k))) go to 30 |
| s = dabs(t(k,k))/dabs(ek-z(k)) |
| call dscal(n,s,z,1) |
| ek = s*ek |
| 30 continue |
| wk = ek - z(k) |
| wkm = -ek - z(k) |
| s = dabs(wk) |
| sm = dabs(wkm) |
| if (t(k,k) .eq. 0.0d0) go to 40 |
| wk = wk/t(k,k) |
| wkm = wkm/t(k,k) |
| go to 50 |
| 40 continue |
| wk = 1.0d0 |
| wkm = 1.0d0 |
| 50 continue |
| if (kk .eq. n) go to 90 |
| j1 = k + 1 |
| if (lower) j1 = 1 |
| j2 = n |
| if (lower) j2 = k - 1 |
| do 60 j = j1, j2 |
| sm = sm + dabs(z(j)+wkm*t(k,j)) |
| z(j) = z(j) + wk*t(k,j) |
| s = s + dabs(z(j)) |
| 60 continue |
| if (s .ge. sm) go to 80 |
| w = wkm - wk |
| wk = wkm |
| do 70 j = j1, j2 |
| z(j) = z(j) + w*t(k,j) |
| 70 continue |
| 80 continue |
| 90 continue |
| z(k) = wk |
| 100 continue |
| s = 1.0d0/dasum(n,z,1) |
| call dscal(n,s,z,1) |
| c |
| ynorm = 1.0d0 |
| c |
| c solve t*z = y |
| c |
| do 130 kk = 1, n |
| k = n + 1 - kk |
| if (lower) k = kk |
| if (dabs(z(k)) .le. dabs(t(k,k))) go to 110 |
| s = dabs(t(k,k))/dabs(z(k)) |
| call dscal(n,s,z,1) |
| ynorm = s*ynorm |
| 110 continue |
| if (t(k,k) .ne. 0.0d0) z(k) = z(k)/t(k,k) |
| if (t(k,k) .eq. 0.0d0) z(k) = 1.0d0 |
| i1 = 1 |
| if (lower) i1 = k + 1 |
| if (kk .ge. n) go to 120 |
| w = -z(k) |
| call daxpy(n-kk,w,t(i1,k),1,z(i1),1) |
| 120 continue |
| 130 continue |
| c make znorm = 1.0 |
| s = 1.0d0/dasum(n,z,1) |
| call dscal(n,s,z,1) |
| ynorm = s*ynorm |
| c |
| if (tnorm .ne. 0.0d0) rcond = ynorm/tnorm |
| if (tnorm .eq. 0.0d0) rcond = 0.0d0 |
| return |
| end |
