c c assorted routines from NETLIB LINPAC archive c note modified array declarations at "cgy" to allow c array bounds checking using f77 -C option c double precision function dasum(n,dx,incx) c c takes the sum of the absolute values. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dtemp integer i,incx,m,mp1,n,nincx c dasum = 0.0d0 dtemp = 0.0d0 if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx dtemp = dtemp + dabs(dx(i)) 10 continue dasum = dtemp return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,6) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dabs(dx(i)) 30 continue if( n .lt. 6 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,6 dtemp = dtemp + dabs(dx(i)) + dabs(dx(i + 1)) + dabs(dx(i + 2)) * + dabs(dx(i + 3)) + dabs(dx(i + 4)) + dabs(dx(i + 5)) 50 continue 60 dasum = dtemp return end subroutine daxpy(n,da,dx,incx,dy,incy) c c constant times a vector plus a vector. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dy(*),da integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if (da .eq. 0.0d0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dy(i) + da*dx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 dy(i) = dy(i) + da*dx(i) dy(i + 1) = dy(i + 1) + da*dx(i + 1) dy(i + 2) = dy(i + 2) + da*dx(i + 2) dy(i + 3) = dy(i + 3) + da*dx(i + 3) 50 continue return end subroutine dgeco(a,lda,n,ipvt,rcond,z) cgy integer lda,n,ipvt(1) cgy double precision a(lda,1),z(1) integer lda,n,ipvt(n) double precision a(lda,n),z(n) double precision rcond c c dgeco factors a double precision matrix by gaussian elimination c and estimates the condition of the matrix. c c if rcond is not needed, dgefa is slightly faster. c to solve a*x = b , follow dgeco by dgesl. c to compute inverse(a)*c , follow dgeco by dgesl. c to compute determinant(a) , follow dgeco by dgedi. c to compute inverse(a) , follow dgeco by dgedi. c c on entry c c a double precision(lda, n) c the matrix to be factored. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c on return c c a an upper triangular matrix and the multipliers c which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c c ipvt integer(n) c an integer vector of pivot indices. c c rcond double precision c an estimate of the reciprocal condition of a . c for the system a*x = b , relative perturbations c in a 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 a 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 a 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 linpack dgefa c blas daxpy,ddot,dscal,dasum c fortran dabs,dmax1,dsign c c internal variables c double precision ddot,ek,t,wk,wkm double precision anorm,s,dasum,sm,ynorm integer info,j,k,kb,kp1,l c c c compute 1-norm of a c anorm = 0.0d0 do 10 j = 1, n anorm = dmax1(anorm,dasum(n,a(1,j),1)) 10 continue c c factor c call dgefa(a,lda,n,ipvt,info) c c rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) . c estimate = norm(z)/norm(y) where a*z = y and trans(a)*y = e . c trans(a) is the transpose of a . the components of e are c chosen to cause maximum local growth in the elements of w where c trans(u)*w = e . the vectors are frequently rescaled to avoid c overflow. c c solve trans(u)*w = e c ek = 1.0d0 do 20 j = 1, n z(j) = 0.0d0 20 continue do 100 k = 1, n if (z(k) .ne. 0.0d0) ek = dsign(ek,-z(k)) if (dabs(ek-z(k)) .le. dabs(a(k,k))) go to 30 s = dabs(a(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 (a(k,k) .eq. 0.0d0) go to 40 wk = wk/a(k,k) wkm = wkm/a(k,k) go to 50 40 continue wk = 1.0d0 wkm = 1.0d0 50 continue kp1 = k + 1 if (kp1 .gt. n) go to 90 do 60 j = kp1, n sm = sm + dabs(z(j)+wkm*a(k,j)) z(j) = z(j) + wk*a(k,j) s = s + dabs(z(j)) 60 continue if (s .ge. sm) go to 80 t = wkm - wk wk = wkm do 70 j = kp1, n z(j) = z(j) + t*a(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 c solve trans(l)*y = w c do 120 kb = 1, n k = n + 1 - kb if (k .lt. n) z(k) = z(k) + ddot(n-k,a(k+1,k),1,z(k+1),1) if (dabs(z(k)) .le. 1.0d0) go to 110 s = 1.0d0/dabs(z(k)) call dscal(n,s,z,1) 110 continue l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t 120 continue s = 1.0d0/dasum(n,z,1) call dscal(n,s,z,1) c ynorm = 1.0d0 c c solve l*v = y c do 140 k = 1, n l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t if (k .lt. n) call daxpy(n-k,t,a(k+1,k),1,z(k+1),1) if (dabs(z(k)) .le. 1.0d0) go to 130 s = 1.0d0/dabs(z(k)) call dscal(n,s,z,1) ynorm = s*ynorm 130 continue 140 continue s = 1.0d0/dasum(n,z,1) call dscal(n,s,z,1) ynorm = s*ynorm c c solve u*z = v c do 160 kb = 1, n k = n + 1 - kb if (dabs(z(k)) .le. dabs(a(k,k))) go to 150 s = dabs(a(k,k))/dabs(z(k)) call dscal(n,s,z,1) ynorm = s*ynorm 150 continue if (a(k,k) .ne. 0.0d0) z(k) = z(k)/a(k,k) if (a(k,k) .eq. 0.0d0) z(k) = 1.0d0 t = -z(k) call daxpy(k-1,t,a(1,k),1,z(1),1) 160 continue c make znorm = 1.0 s = 1.0d0/dasum(n,z,1) call dscal(n,s,z,1) ynorm = s*ynorm c if (anorm .ne. 0.0d0) rcond = ynorm/anorm if (anorm .eq. 0.0d0) rcond = 0.0d0 return end subroutine dgesl(a,lda,n,ipvt,b,job) cgy integer lda,n,ipvt(1),job cgy double precision a(lda,1),b(1) integer lda,n,ipvt(n),job double precision a(lda,n),b(n) c c dgesl solves the double precision system c a * x = b or trans(a) * x = b c using the factors computed by dgeco or dgefa. c c on entry c c a double precision(lda, n) c the output from dgeco or dgefa. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c ipvt integer(n) c the pivot vector from dgeco or dgefa. c c b double precision(n) c the right hand side vector. c c job integer c = 0 to solve a*x = b , c = nonzero to solve trans(a)*x = b where c trans(a) is the transpose. c c on return c c b the solution vector x . c c error condition c c a division by zero will occur if the input factor contains a c zero on the diagonal. technically this indicates singularity c but it is often caused by improper arguments or improper c setting of lda . it will not occur if the subroutines are c called correctly and if dgeco has set rcond .gt. 0.0 c or dgefa has set info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call dgeco(a,lda,n,ipvt,rcond,z) c if (rcond is too small) go to ... c do 10 j = 1, p c call dgesl(a,lda,n,ipvt,c(1,j),0) c 10 continue 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,ddot c c internal variables c double precision ddot,t integer k,kb,l,nm1 c nm1 = n - 1 if (job .ne. 0) go to 50 c c job = 0 , solve a * x = b c first solve l*y = b c if (nm1 .lt. 1) go to 30 do 20 k = 1, nm1 l = ipvt(k) t = b(l) if (l .eq. k) go to 10 b(l) = b(k) b(k) = t 10 continue call daxpy(n-k,t,a(k+1,k),1,b(k+1),1) 20 continue 30 continue c c now solve u*x = y c do 40 kb = 1, n k = n + 1 - kb b(k) = b(k)/a(k,k) t = -b(k) call daxpy(k-1,t,a(1,k),1,b(1),1) 40 continue go to 100 50 continue c c job = nonzero, solve trans(a) * x = b c first solve trans(u)*y = b c do 60 k = 1, n t = ddot(k-1,a(1,k),1,b(1),1) b(k) = (b(k) - t)/a(k,k) 60 continue c c now solve trans(l)*x = y c if (nm1 .lt. 1) go to 90 do 80 kb = 1, nm1 k = n - kb b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1) l = ipvt(k) if (l .eq. k) go to 70 t = b(l) b(l) = b(k) b(k) = t 70 continue 80 continue 90 continue 100 continue return end subroutine dscal(n,da,dx,incx) c c scales a vector by a constant. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c double precision da,dx(*) integer i,incx,m,mp1,n,nincx c if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx dx(i) = da*dx(i) 10 continue return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dx(i) = da*dx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 dx(i) = da*dx(i) dx(i + 1) = da*dx(i + 1) dx(i + 2) = da*dx(i + 2) dx(i + 3) = da*dx(i + 3) dx(i + 4) = da*dx(i + 4) 50 continue return end *DECK DSDOT DOUBLE PRECISION FUNCTION DSDOT (N, SX, INCX, SY, INCY) C***BEGIN PROLOGUE DSDOT C***PURPOSE Compute the inner product of two vectors with extended C precision accumulation and result. C***LIBRARY SLATEC (BLAS) C***CATEGORY D1A4 C***TYPE DOUBLE PRECISION (DSDOT-D, DCDOT-C) C***KEYWORDS BLAS, COMPLEX VECTORS, DOT PRODUCT, INNER PRODUCT, C LINEAR ALGEBRA, VECTOR C***AUTHOR Lawson, C. L., (JPL) C Hanson, R. J., (SNLA) C Kincaid, D. R., (U. of Texas) C Krogh, F. T., (JPL) C***DESCRIPTION C C B L A S Subprogram C Description of Parameters C C --Input-- C N number of elements in input vector(s) C SX single precision vector with N elements C INCX storage spacing between elements of SX C SY single precision vector with N elements C INCY storage spacing between elements of SY C C --Output-- C DSDOT double precision dot product (zero if N.LE.0) C C Returns D.P. dot product accumulated in D.P., for S.P. SX and SY C DSDOT = sum for I = 0 to N-1 of SX(LX+I*INCX) * SY(LY+I*INCY), C where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is C defined in a similar way using INCY. C C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T. C Krogh, Basic linear algebra subprograms for Fortran C usage, Algorithm No. 539, Transactions on Mathematical C Software 5, 3 (September 1979), pp. 308-323. C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 791001 DATE WRITTEN C 890831 Modified array declarations. (WRB) C 890831 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 920310 Corrected definition of LX in DESCRIPTION. (WRB) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE DSDOT REAL SX(*),SY(*) C***FIRST EXECUTABLE STATEMENT DSDOT DSDOT = 0.0D0 IF (N .LE. 0) RETURN IF (INCX.EQ.INCY .AND. INCX.GT.0) GO TO 20 C C Code for unequal or nonpositive increments. C KX = 1 KY = 1 IF (INCX .LT. 0) KX = 1+(1-N)*INCX IF (INCY .LT. 0) KY = 1+(1-N)*INCY DO 10 I = 1,N DSDOT = DSDOT + DBLE(SX(KX))*DBLE(SY(KY)) KX = KX + INCX KY = KY + INCY 10 CONTINUE RETURN C C Code for equal, positive, non-unit increments. C 20 NS = N*INCX DO 30 I = 1,NS,INCX DSDOT = DSDOT + DBLE(SX(I))*DBLE(SY(I)) 30 CONTINUE RETURN END integer function isamax(n,sx,incx) c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c real sx(*),smax integer i,incx,ix,n c isamax = 0 if( n.lt.1 .or. incx.le.0 ) return isamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 smax = abs(sx(1)) ix = ix + incx do 10 i = 2,n if(abs(sx(ix)).le.smax) go to 5 isamax = i smax = abs(sx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 smax = abs(sx(1)) do 30 i = 2,n if(abs(sx(i)).le.smax) go to 30 isamax = i smax = abs(sx(i)) 30 continue return end real function sasum(n,sx,incx) c c takes the sum of the absolute values. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c real sx(*),stemp integer i,incx,m,mp1,n,nincx c sasum = 0.0e0 stemp = 0.0e0 if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx stemp = stemp + abs(sx(i)) 10 continue sasum = stemp return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,6) if( m .eq. 0 ) go to 40 do 30 i = 1,m stemp = stemp + abs(sx(i)) 30 continue if( n .lt. 6 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,6 stemp = stemp + abs(sx(i)) + abs(sx(i + 1)) + abs(sx(i + 2)) * + abs(sx(i + 3)) + abs(sx(i + 4)) + abs(sx(i + 5)) 50 continue 60 sasum = stemp return end subroutine saxpy(n,sa,sx,incx,sy,incy) c c constant times a vector plus a vector. c uses unrolled loop for increments equal to one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c real sx(*),sy(*),sa integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if (sa .eq. 0.0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n sy(iy) = sy(iy) + sa*sx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m sy(i) = sy(i) + sa*sx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 sy(i) = sy(i) + sa*sx(i) sy(i + 1) = sy(i + 1) + sa*sx(i + 1) sy(i + 2) = sy(i + 2) + sa*sx(i + 2) sy(i + 3) = sy(i + 3) + sa*sx(i + 3) 50 continue return end real function sdot(n,sx,incx,sy,incy) c c forms the dot product of two vectors. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c real sx(*),sy(*),stemp integer i,incx,incy,ix,iy,m,mp1,n c stemp = 0.0e0 sdot = 0.0e0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n stemp = stemp + sx(ix)*sy(iy) ix = ix + incx iy = iy + incy 10 continue sdot = stemp return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m stemp = stemp + sx(i)*sy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 stemp = stemp + sx(i)*sy(i) + sx(i + 1)*sy(i + 1) + * sx(i + 2)*sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4) 50 continue 60 sdot = stemp return end subroutine sgefa(a,lda,n,ipvt,info) cgy integer lda,n,ipvt(1),info cgy real a(lda,1) integer lda,n,ipvt(n),info real a(lda,n) c c sgefa factors a real matrix by gaussian elimination. c c sgefa is usually called by sgeco, but it can be called c directly with a saving in time if rcond is not needed. c (time for sgeco) = (1 + 9/n)*(time for sgefa) . c c on entry c c a real(lda, n) c the matrix to be factored. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c on return c c a an upper triangular matrix and the multipliers c which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c c ipvt integer(n) c an integer vector of pivot indices. c c info integer c = 0 normal value. c = k if u(k,k) .eq. 0.0 . this is not an error c condition for this subroutine, but it does c indicate that sgesl or sgedi will divide by zero c if called. use rcond in sgeco for a reliable c indication of singularity. 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 saxpy,sscal,isamax c c internal variables c real t integer isamax,j,k,kp1,l,nm1 c c c gaussian elimination with partial pivoting c info = 0 nm1 = n - 1 if (nm1 .lt. 1) go to 70 do 60 k = 1, nm1 kp1 = k + 1 c c find l = pivot index c l = isamax(n-k+1,a(k,k),1) + k - 1 ipvt(k) = l c c zero pivot implies this column already triangularized c if (a(l,k) .eq. 0.0e0) go to 40 c c interchange if necessary c if (l .eq. k) go to 10 t = a(l,k) a(l,k) = a(k,k) a(k,k) = t 10 continue c c compute multipliers c t = -1.0e0/a(k,k) call sscal(n-k,t,a(k+1,k),1) c c row elimination with column indexing c do 30 j = kp1, n t = a(l,j) if (l .eq. k) go to 20 a(l,j) = a(k,j) a(k,j) = t 20 continue call saxpy(n-k,t,a(k+1,k),1,a(k+1,j),1) 30 continue go to 50 40 continue info = k 50 continue 60 continue 70 continue ipvt(n) = n if (a(n,n) .eq. 0.0e0) info = n return end subroutine sgesl(a,lda,n,ipvt,b,job) cgy integer lda,n,ipvt(1),job cgy real a(lda,1),b(1) integer lda,n,ipvt(n),job real a(lda,n),b(n) c c sgesl solves the real system c a * x = b or trans(a) * x = b c using the factors computed by sgeco or sgefa. c c on entry c c a real(lda, n) c the output from sgeco or sgefa. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c ipvt integer(n) c the pivot vector from sgeco or sgefa. c c b real(n) c the right hand side vector. c c job integer c = 0 to solve a*x = b , c = nonzero to solve trans(a)*x = b where c trans(a) is the transpose. c c on return c c b the solution vector x . c c error condition c c a division by zero will occur if the input factor contains a c zero on the diagonal. technically this indicates singularity c but it is often caused by improper arguments or improper c setting of lda . it will not occur if the subroutines are c called correctly and if sgeco has set rcond .gt. 0.0 c or sgefa has set info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call sgeco(a,lda,n,ipvt,rcond,z) c if (rcond is too small) go to ... c do 10 j = 1, p c call sgesl(a,lda,n,ipvt,c(1,j),0) c 10 continue 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 saxpy,sdot c c internal variables c real sdot,t integer k,kb,l,nm1 c nm1 = n - 1 if (job .ne. 0) go to 50 c c job = 0 , solve a * x = b c first solve l*y = b c if (nm1 .lt. 1) go to 30 do 20 k = 1, nm1 l = ipvt(k) t = b(l) if (l .eq. k) go to 10 b(l) = b(k) b(k) = t 10 continue call saxpy(n-k,t,a(k+1,k),1,b(k+1),1) 20 continue 30 continue c c now solve u*x = y c do 40 kb = 1, n k = n + 1 - kb b(k) = b(k)/a(k,k) t = -b(k) call saxpy(k-1,t,a(1,k),1,b(1),1) 40 continue go to 100 50 continue c c job = nonzero, solve trans(a) * x = b c first solve trans(u)*y = b c do 60 k = 1, n t = sdot(k-1,a(1,k),1,b(1),1) b(k) = (b(k) - t)/a(k,k) 60 continue c c now solve trans(l)*x = y c if (nm1 .lt. 1) go to 90 do 80 kb = 1, nm1 k = n - kb b(k) = b(k) + sdot(n-k,a(k+1,k),1,b(k+1),1) l = ipvt(k) if (l .eq. k) go to 70 t = b(l) b(l) = b(k) b(k) = t 70 continue 80 continue 90 continue 100 continue return end subroutine sscal(n,sa,sx,incx) c c scales a vector by a constant. c uses unrolled loops for increment equal to 1. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c real sa,sx(*) integer i,incx,m,mp1,n,nincx c if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx sx(i) = sa*sx(i) 10 continue return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m sx(i) = sa*sx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 sx(i) = sa*sx(i) sx(i + 1) = sa*sx(i + 1) sx(i + 2) = sa*sx(i + 2) sx(i + 3) = sa*sx(i + 3) sx(i + 4) = sa*sx(i + 4) 50 continue return end subroutine sgbfa(abd,lda,n,ml,mu,ipvt,info) cgy integer lda,n,ml,mu,ipvt(1),info cgy real abd(lda,1) integer lda,n,ml,mu,ipvt(n),info real abd(lda,n) c c sgbfa factors a real band matrix by elimination. c c sgbfa is usually called by sgbco, but it can be called c directly with a saving in time if rcond is not needed. c c on entry c c abd real(lda, n) c contains the matrix in band storage. the columns c of the matrix are stored in the columns of abd and c the diagonals of the matrix are stored in rows c ml+1 through 2*ml+mu+1 of abd . c see the comments below for details. c c lda integer c the leading dimension of the array abd . c lda must be .ge. 2*ml + mu + 1 . c c n integer c the order of the original matrix. c c ml integer c number of diagonals below the main diagonal. c 0 .le. ml .lt. n . c c mu integer c number of diagonals above the main diagonal. c 0 .le. mu .lt. n . c more efficient if ml .le. mu . c on return c c abd an upper triangular matrix in band storage and c the multipliers which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c c ipvt integer(n) c an integer vector of pivot indices. c c info integer c = 0 normal value. c = k if u(k,k) .eq. 0.0 . this is not an error c condition for this subroutine, but it does c indicate that sgbsl will divide by zero if c called. use rcond in sgbco for a reliable c indication of singularity. c c band storage c c if a is a band matrix, the following program segment c will set up the input. c c ml = (band width below the diagonal) c mu = (band width above the diagonal) c m = ml + mu + 1 c do 20 j = 1, n c i1 = max0(1, j-mu) c i2 = min0(n, j+ml) c do 10 i = i1, i2 c k = i - j + m c abd(k,j) = a(i,j) c 10 continue c 20 continue c c this uses rows ml+1 through 2*ml+mu+1 of abd . c in addition, the first ml rows in abd are used for c elements generated during the triangularization. c the total number of rows needed in abd is 2*ml+mu+1 . c the ml+mu by ml+mu upper left triangle and the c ml by ml lower right triangle are not referenced. 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 saxpy,sscal,isamax c fortran max0,min0 c c internal variables c real t integer i,isamax,i0,j,ju,jz,j0,j1,k,kp1,l,lm,m,mm,nm1 c c m = ml + mu + 1 info = 0 c c zero initial fill-in columns c j0 = mu + 2 j1 = min0(n,m) - 1 if (j1 .lt. j0) go to 30 do 20 jz = j0, j1 i0 = m + 1 - jz do 10 i = i0, ml abd(i,jz) = 0.0e0 10 continue 20 continue 30 continue jz = j1 ju = 0 c c gaussian elimination with partial pivoting c nm1 = n - 1 if (nm1 .lt. 1) go to 130 do 120 k = 1, nm1 kp1 = k + 1 c c zero next fill-in column c jz = jz + 1 if (jz .gt. n) go to 50 if (ml .lt. 1) go to 50 do 40 i = 1, ml abd(i,jz) = 0.0e0 40 continue 50 continue c c find l = pivot index c lm = min0(ml,n-k) l = isamax(lm+1,abd(m,k),1) + m - 1 ipvt(k) = l + k - m c c zero pivot implies this column already triangularized c if (abd(l,k) .eq. 0.0e0) go to 100 c c interchange if necessary c if (l .eq. m) go to 60 t = abd(l,k) abd(l,k) = abd(m,k) abd(m,k) = t 60 continue c c compute multipliers c t = -1.0e0/abd(m,k) call sscal(lm,t,abd(m+1,k),1) c c row elimination with column indexing c ju = min0(max0(ju,mu+ipvt(k)),n) mm = m if (ju .lt. kp1) go to 90 do 80 j = kp1, ju l = l - 1 mm = mm - 1 t = abd(l,j) if (l .eq. mm) go to 70 abd(l,j) = abd(mm,j) abd(mm,j) = t 70 continue call saxpy(lm,t,abd(m+1,k),1,abd(mm+1,j),1) 80 continue 90 continue go to 110 100 continue info = k 110 continue 120 continue 130 continue ipvt(n) = n if (abd(m,n) .eq. 0.0e0) info = n return end subroutine sgbsl(abd,lda,n,ml,mu,ipvt,b,job) cgy integer lda,n,ml,mu,ipvt(1),job cgy real abd(lda,1),b(1) integer lda,n,ml,mu,ipvt(n),job real abd(lda,n),b(n) c c sgbsl solves the real band system c a * x = b or trans(a) * x = b c using the factors computed by sgbco or sgbfa. c c on entry c c abd real(lda, n) c the output from sgbco or sgbfa. c c lda integer c the leading dimension of the array abd . c c n integer c the order of the original matrix. c c ml integer c number of diagonals below the main diagonal. c c mu integer c number of diagonals above the main diagonal. c c ipvt integer(n) c the pivot vector from sgbco or sgbfa. c c b real(n) c the right hand side vector. c c job integer c = 0 to solve a*x = b , c = nonzero to solve trans(a)*x = b , where c trans(a) is the transpose. c c on return c c b the solution vector x . c c error condition c c a division by zero will occur if the input factor contains a c zero on the diagonal. technically this indicates singularity c but it is often caused by improper arguments or improper c setting of lda . it will not occur if the subroutines are c called correctly and if sgbco has set rcond .gt. 0.0 c or sgbfa has set info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call sgbco(abd,lda,n,ml,mu,ipvt,rcond,z) c if (rcond is too small) go to ... c do 10 j = 1, p c call sgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0) c 10 continue 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 saxpy,sdot c fortran min0 c c internal variables c real sdot,t integer k,kb,l,la,lb,lm,m,nm1 c m = mu + ml + 1 nm1 = n - 1 if (job .ne. 0) go to 50 c c job = 0 , solve a * x = b c first solve l*y = b c if (ml .eq. 0) go to 30 if (nm1 .lt. 1) go to 30 do 20 k = 1, nm1 lm = min0(ml,n-k) l = ipvt(k) t = b(l) if (l .eq. k) go to 10 b(l) = b(k) b(k) = t 10 continue call saxpy(lm,t,abd(m+1,k),1,b(k+1),1) 20 continue 30 continue c c now solve u*x = y c do 40 kb = 1, n k = n + 1 - kb b(k) = b(k)/abd(m,k) lm = min0(k,m) - 1 la = m - lm lb = k - lm t = -b(k) call saxpy(lm,t,abd(la,k),1,b(lb),1) 40 continue go to 100 50 continue c c job = nonzero, solve trans(a) * x = b c first solve trans(u)*y = b c do 60 k = 1, n lm = min0(k,m) - 1 la = m - lm lb = k - lm t = sdot(lm,abd(la,k),1,b(lb),1) b(k) = (b(k) - t)/abd(m,k) 60 continue c c now solve trans(l)*x = y c if (ml .eq. 0) go to 90 if (nm1 .lt. 1) go to 90 do 80 kb = 1, nm1 k = n - kb lm = min0(ml,n-k) b(k) = b(k) + sdot(lm,abd(m+1,k),1,b(k+1),1) l = ipvt(k) if (l .eq. k) go to 70 t = b(l) b(l) = b(k) b(k) = t 70 continue 80 continue 90 continue 100 continue return end subroutine scopy(n,sx,incx,sy,incy) c c copies a vector, x, to a vector, y. c uses unrolled loops for increments equal to 1. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c real sx(*),sy(*) integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n sy(iy) = sx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,7) if( m .eq. 0 ) go to 40 do 30 i = 1,m sy(i) = sx(i) 30 continue if( n .lt. 7 ) return 40 mp1 = m + 1 do 50 i = mp1,n,7 sy(i) = sx(i) sy(i + 1) = sx(i + 1) sy(i + 2) = sx(i + 2) sy(i + 3) = sx(i + 3) sy(i + 4) = sx(i + 4) sy(i + 5) = sx(i + 5) sy(i + 6) = sx(i + 6) 50 continue return end double precision function ddot(n,dx,incx,dy,incy) c c forms the dot product of two vectors. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dy(*),dtemp integer i,incx,incy,ix,iy,m,mp1,n c ddot = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dtemp + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy 10 continue ddot = dtemp return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dx(i)*dy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 50 continue 60 ddot = dtemp return end subroutine dgbfa(abd,lda,n,ml,mu,ipvt,info) cgy integer lda,n,ml,mu,ipvt(1),info cgy double precision abd(lda,1) integer lda,n,ml,mu,ipvt(n),info double precision abd(lda,n) c c dgbfa factors a double precision band matrix by elimination. c c dgbfa is usually called by dgbco, but it can be called c directly with a saving in time if rcond is not needed. c c on entry c c abd double precision(lda, n) c contains the matrix in band storage. the columns c of the matrix are stored in the columns of abd and c the diagonals of the matrix are stored in rows c ml+1 through 2*ml+mu+1 of abd . c see the comments below for details. c c lda integer c the leading dimension of the array abd . c lda must be .ge. 2*ml + mu + 1 . c c n integer c the order of the original matrix. c c ml integer c number of diagonals below the main diagonal. c 0 .le. ml .lt. n . c c mu integer c number of diagonals above the main diagonal. c 0 .le. mu .lt. n . c more efficient if ml .le. mu . c on return c c abd an upper triangular matrix in band storage and c the multipliers which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c c ipvt integer(n) c an integer vector of pivot indices. c c info integer c = 0 normal value. c = k if u(k,k) .eq. 0.0 . this is not an error c condition for this subroutine, but it does c indicate that dgbsl will divide by zero if c called. use rcond in dgbco for a reliable c indication of singularity. c c band storage c c if a is a band matrix, the following program segment c will set up the input. c c ml = (band width below the diagonal) c mu = (band width above the diagonal) c m = ml + mu + 1 c do 20 j = 1, n c i1 = max0(1, j-mu) c i2 = min0(n, j+ml) c do 10 i = i1, i2 c k = i - j + m c abd(k,j) = a(i,j) c 10 continue c 20 continue c c this uses rows ml+1 through 2*ml+mu+1 of abd . c in addition, the first ml rows in abd are used for c elements generated during the triangularization. c the total number of rows needed in abd is 2*ml+mu+1 . c the ml+mu by ml+mu upper left triangle and the c ml by ml lower right triangle are not referenced. 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,idamax c fortran max0,min0 c c internal variables c double precision t integer i,idamax,i0,j,ju,jz,j0,j1,k,kp1,l,lm,m,mm,nm1 c c m = ml + mu + 1 info = 0 c c zero initial fill-in columns c j0 = mu + 2 j1 = min0(n,m) - 1 if (j1 .lt. j0) go to 30 do 20 jz = j0, j1 i0 = m + 1 - jz do 10 i = i0, ml abd(i,jz) = 0.0d0 10 continue 20 continue 30 continue jz = j1 ju = 0 c c gaussian elimination with partial pivoting c nm1 = n - 1 if (nm1 .lt. 1) go to 130 do 120 k = 1, nm1 kp1 = k + 1 c c zero next fill-in column c jz = jz + 1 if (jz .gt. n) go to 50 if (ml .lt. 1) go to 50 do 40 i = 1, ml abd(i,jz) = 0.0d0 40 continue 50 continue c c find l = pivot index c lm = min0(ml,n-k) l = idamax(lm+1,abd(m,k),1) + m - 1 ipvt(k) = l + k - m c c zero pivot implies this column already triangularized c if (abd(l,k) .eq. 0.0d0) go to 100 c c interchange if necessary c if (l .eq. m) go to 60 t = abd(l,k) abd(l,k) = abd(m,k) abd(m,k) = t 60 continue c c compute multipliers c t = -1.0d0/abd(m,k) call dscal(lm,t,abd(m+1,k),1) c c row elimination with column indexing c ju = min0(max0(ju,mu+ipvt(k)),n) mm = m if (ju .lt. kp1) go to 90 do 80 j = kp1, ju l = l - 1 mm = mm - 1 t = abd(l,j) if (l .eq. mm) go to 70 abd(l,j) = abd(mm,j) abd(mm,j) = t 70 continue call daxpy(lm,t,abd(m+1,k),1,abd(mm+1,j),1) 80 continue 90 continue go to 110 100 continue info = k 110 continue 120 continue 130 continue ipvt(n) = n if (abd(m,n) .eq. 0.0d0) info = n return end subroutine dgefa(a,lda,n,ipvt,info) cgy integer lda,n,ipvt(1),info cgy double precision a(lda,1) integer lda,n,ipvt(n),info double precision a(lda,n) c c dgefa factors a double precision matrix by gaussian elimination. c c dgefa is usually called by dgeco, but it can be called c directly with a saving in time if rcond is not needed. c (time for dgeco) = (1 + 9/n)*(time for dgefa) . c c on entry c c a double precision(lda, n) c the matrix to be factored. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c on return c c a an upper triangular matrix and the multipliers c which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c c ipvt integer(n) c an integer vector of pivot indices. c c info integer c = 0 normal value. c = k if u(k,k) .eq. 0.0 . this is not an error c condition for this subroutine, but it does c indicate that dgesl or dgedi will divide by zero c if called. use rcond in dgeco for a reliable c indication of singularity. 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,idamax c c internal variables c double precision t integer idamax,j,k,kp1,l,nm1 c c c gaussian elimination with partial pivoting c info = 0 nm1 = n - 1 if (nm1 .lt. 1) go to 70 do 60 k = 1, nm1 kp1 = k + 1 c c find l = pivot index c l = idamax(n-k+1,a(k,k),1) + k - 1 ipvt(k) = l c c zero pivot implies this column already triangularized c if (a(l,k) .eq. 0.0d0) go to 40 c c interchange if necessary c if (l .eq. k) go to 10 t = a(l,k) a(l,k) = a(k,k) a(k,k) = t 10 continue c c compute multipliers c t = -1.0d0/a(k,k) call dscal(n-k,t,a(k+1,k),1) c c row elimination with column indexing c do 30 j = kp1, n t = a(l,j) if (l .eq. k) go to 20 a(l,j) = a(k,j) a(k,j) = t 20 continue call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1) 30 continue go to 50 40 continue info = k 50 continue 60 continue 70 continue ipvt(n) = n if (a(n,n) .eq. 0.0d0) info = n return end integer function idamax(n,dx,incx) c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dmax integer i,incx,ix,n c idamax = 0 if( n.lt.1 .or. incx.le.0 ) return idamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 dmax = dabs(dx(1)) ix = ix + incx do 10 i = 2,n if(dabs(dx(ix)).le.dmax) go to 5 idamax = i dmax = dabs(dx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 dmax = dabs(dx(1)) do 30 i = 2,n if(dabs(dx(i)).le.dmax) go to 30 idamax = i dmax = dabs(dx(i)) 30 continue return end subroutine dgbsl(abd,lda,n,ml,mu,ipvt,b,job) cgy integer lda,n,ml,mu,ipvt(1),job cgy double precision abd(lda,1),b(1) integer lda,n,ml,mu,ipvt(n),job double precision abd(lda,n),b(n) c c dgbsl solves the double precision band system c a * x = b or trans(a) * x = b c using the factors computed by dgbco or dgbfa. c c on entry c c abd double precision(lda, n) c the output from dgbco or dgbfa. c c lda integer c the leading dimension of the array abd . c c n integer c the order of the original matrix. c c ml integer c number of diagonals below the main diagonal. c c mu integer c number of diagonals above the main diagonal. c c ipvt integer(n) c the pivot vector from dgbco or dgbfa. c c b double precision(n) c the right hand side vector. c c job integer c = 0 to solve a*x = b , c = nonzero to solve trans(a)*x = b , where c trans(a) is the transpose. c c on return c c b the solution vector x . c c error condition c c a division by zero will occur if the input factor contains a c zero on the diagonal. technically this indicates singularity c but it is often caused by improper arguments or improper c setting of lda . it will not occur if the subroutines are c called correctly and if dgbco has set rcond .gt. 0.0 c or dgbfa has set info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call dgbco(abd,lda,n,ml,mu,ipvt,rcond,z) c if (rcond is too small) go to ... c do 10 j = 1, p c call dgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0) c 10 continue 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,ddot c fortran min0 c c internal variables c double precision ddot,t integer k,kb,l,la,lb,lm,m,nm1 c m = mu + ml + 1 nm1 = n - 1 if (job .ne. 0) go to 50 c c job = 0 , solve a * x = b c first solve l*y = b c if (ml .eq. 0) go to 30 if (nm1 .lt. 1) go to 30 do 20 k = 1, nm1 lm = min0(ml,n-k) l = ipvt(k) t = b(l) if (l .eq. k) go to 10 b(l) = b(k) b(k) = t 10 continue call daxpy(lm,t,abd(m+1,k),1,b(k+1),1) 20 continue 30 continue c c now solve u*x = y c do 40 kb = 1, n k = n + 1 - kb b(k) = b(k)/abd(m,k) lm = min0(k,m) - 1 la = m - lm lb = k - lm t = -b(k) call daxpy(lm,t,abd(la,k),1,b(lb),1) 40 continue go to 100 50 continue c c job = nonzero, solve trans(a) * x = b c first solve trans(u)*y = b c do 60 k = 1, n lm = min0(k,m) - 1 la = m - lm lb = k - lm t = ddot(lm,abd(la,k),1,b(lb),1) b(k) = (b(k) - t)/abd(m,k) 60 continue c c now solve trans(l)*x = y c if (ml .eq. 0) go to 90 if (nm1 .lt. 1) go to 90 do 80 kb = 1, nm1 k = n - kb lm = min0(ml,n-k) b(k) = b(k) + ddot(lm,abd(m+1,k),1,b(k+1),1) l = ipvt(k) if (l .eq. k) go to 70 t = b(l) b(l) = b(k) b(k) = t 70 continue 80 continue 90 continue 100 continue return end double precision function d1mach(i) integer i c c double-precision machine constants c d1mach( 1) = b**(emin-1), the smallest positive magnitude. c d1mach( 2) = b**emax*(1 - b**(-t)), the largest magnitude. c d1mach( 3) = b**(-t), the smallest relative spacing. c d1mach( 4) = b**(1-t), the largest relative spacing. c d1mach( 5) = log10(b) c integer small(2) integer large(2) integer right(2) integer diver(2) integer log10(2) integer sc, cray1(38), j common /d9mach/ cray1 save small, large, right, diver, log10, sc double precision dmach(5) double precision b, x equivalence (dmach(1),small(1)) equivalence (dmach(2),large(1)) equivalence (dmach(3),right(1)) equivalence (dmach(4),diver(1)) equivalence (dmach(5),log10(1)) c this version adapts automatically to most current machines. c r1mach can handle auto-double compiling, but this version of c d1mach does not, because we do not have quad constants for c many machines yet. c to compile on older machines, add a c in column 1 c on the next line data sc/0/ c and remove the c from column 1 in one of the sections below. c constants for even older machines can be obtained by c mail netlib@research.bell-labs.com c send old1mach from blas c please send corrections to dmg or ehg@bell-labs.com. c c machine constants for the honeywell dps 8/70 series. c data small(1),small(2) / o402400000000, o000000000000 / c data large(1),large(2) / o376777777777, o777777777777 / c data right(1),right(2) / o604400000000, o000000000000 / c data diver(1),diver(2) / o606400000000, o000000000000 / c data log10(1),log10(2) / o776464202324, o117571775714 /, sc/987/ c c machine constants for pdp-11 fortrans supporting c 32-bit integers. c data small(1),small(2) / 8388608, 0 / c data large(1),large(2) / 2147483647, -1 / c data right(1),right(2) / 612368384, 0 / c data diver(1),diver(2) / 620756992, 0 / c data log10(1),log10(2) / 1067065498, -2063872008 /, sc/987/ c c machine constants for the univac 1100 series. c data small(1),small(2) / o000040000000, o000000000000 / c data large(1),large(2) / o377777777777, o777777777777 / c data right(1),right(2) / o170540000000, o000000000000 / c data diver(1),diver(2) / o170640000000, o000000000000 / c data log10(1),log10(2) / o177746420232, o411757177572 /, sc/987/ c c on first call, if no data uncommented, test machine types. if (sc .ne. 987) then dmach(1) = 1.d13 if ( small(1) .eq. 1117925532 * .and. small(2) .eq. -448790528) then * *** ieee big endian *** small(1) = 1048576 small(2) = 0 large(1) = 2146435071 large(2) = -1 right(1) = 1017118720 right(2) = 0 diver(1) = 1018167296 diver(2) = 0 log10(1) = 1070810131 log10(2) = 1352628735 else if ( small(2) .eq. 1117925532 * .and. small(1) .eq. -448790528) then * *** ieee little endian *** small(2) = 1048576 small(1) = 0 large(2) = 2146435071 large(1) = -1 right(2) = 1017118720 right(1) = 0 diver(2) = 1018167296 diver(1) = 0 log10(2) = 1070810131 log10(1) = 1352628735 else if ( small(1) .eq. -2065213935 * .and. small(2) .eq. 10752) then * *** vax with d_floating *** small(1) = 128 small(2) = 0 large(1) = -32769 large(2) = -1 right(1) = 9344 right(2) = 0 diver(1) = 9472 diver(2) = 0 log10(1) = 546979738 log10(2) = -805796613 else if ( small(1) .eq. 1267827943 * .and. small(2) .eq. 704643072) then * *** ibm mainframe *** small(1) = 1048576 small(2) = 0 large(1) = 2147483647 large(2) = -1 right(1) = 856686592 right(2) = 0 diver(1) = 873463808 diver(2) = 0 log10(1) = 1091781651 log10(2) = 1352628735 else if ( small(1) .eq. 1120022684 * .and. small(2) .eq. -448790528) then * *** convex c-1 *** small(1) = 1048576 small(2) = 0 large(1) = 2147483647 large(2) = -1 right(1) = 1019215872 right(2) = 0 diver(1) = 1020264448 diver(2) = 0 log10(1) = 1072907283 log10(2) = 1352628735 else if ( small(1) .eq. 815547074 * .and. small(2) .eq. 58688) then * *** vax g-floating *** small(1) = 16 small(2) = 0 large(1) = -32769 large(2) = -1 right(1) = 15552 right(2) = 0 diver(1) = 15568 diver(2) = 0 log10(1) = 1142112243 log10(2) = 2046775455 else c c------------------------------------------ c --- Modified by Gary Wilson / Ramboll c Fri Oct 29 16:52:14 PDT 2010 c Use F90 Intrinsic functions to get c the value of d1mach(4) c------------------------------------------ c x = 1.0D0 b = RADIX(X) dmach(4) = b**(1-DIGITS(x)) if( i .NE. 4 ) then write(iout,'(//,A)') 'ERROR in LINPAC:' write(iout,'(/,2A)') ' This version only supports ', & 'D1MACH for an argument of 4.' write(iout,'(2A)') ' Other values returned from ', & 'D1MACH are not accurate.' call camxerr() endif c c------------------------------------------ c CENVIRON dmach(2) = 1.d27 + 1 CENVIRON dmach(3) = 1.d27 CENVIRON large(2) = large(2) - right(2) CENVIRON if (large(2) .eq. 64 .and. small(2) .eq. 0) then CENVIRON cray1(1) = 67291416 CENVIRON do 10 j = 1, 20 CENVIRON cray1(j+1) = cray1(j) + cray1(j) CENVIRON 10 continue CENVIRON cray1(22) = cray1(21) + 321322 CENVIRON do 20 j = 22, 37 CENVIRON cray1(j+1) = cray1(j) + cray1(j) CENVIRON 20 continue CENVIRON if (cray1(38) .eq. small(1)) then CENVIRON* *** cray *** CENVIRON call i1mcry(small(1), j, 8285, 8388608, 0) CENVIRON small(2) = 0 CENVIRON call i1mcry(large(1), j, 24574, 16777215, 16777215) CENVIRON call i1mcry(large(2), j, 0, 16777215, 16777214) CENVIRON call i1mcry(right(1), j, 16291, 8388608, 0) CENVIRON right(2) = 0 CENVIRON call i1mcry(diver(1), j, 16292, 8388608, 0) CENVIRON diver(2) = 0 CENVIRON call i1mcry(log10(1), j, 16383, 10100890, 8715215) CENVIRON call i1mcry(log10(2), j, 0, 16226447, 9001388) CENVIRON else CENVIRON write(*,9000) CENVIRON stop 779 CENVIRON end if CENVIRON else CENVIRON write(*,9000) CENVIRON stop 779 CENVIRON end if end if sc = 987 end if * sanity check if (dmach(4) .ge. 1.0d0) stop 778 if (i .lt. 1 .or. i .gt. 5) then write(*,*) 'd1mach(i): i =',i,' is out of bounds.' stop end if d1mach = dmach(i) return 9000 format(/' adjust d1mach by uncommenting data statements'/ *' appropriate for your machine.') * /* standard c source for d1mach -- remove the * in column 1 */ *#include *#include *#include *double d1mach_(long *i) *{ * switch(*i){ * case 1: return dbl_min; * case 2: return dbl_max; * case 3: return dbl_epsilon/flt_radix; * case 4: return dbl_epsilon; * case 5: return log10((double)flt_radix); * } * fprintf(stderr, "invalid argument: d1mach(%ld)\n", *i); * exit(1); return 0; /* some compilers demand return values */ *} end subroutine i1mcry(a, a1, b, c, d) **** special computation for old cray machines **** integer a, a1, b, c, d a1 = 16777216*b + c a = 16777216*a1 + d end real function r1mach(i) integer i c c single-precision machine constants c r1mach(1) = b**(emin-1), the smallest positive magnitude. c r1mach(2) = b**emax*(1 - b**(-t)), the largest magnitude. c r1mach(3) = b**(-t), the smallest relative spacing. c r1mach(4) = b**(1-t), the largest relative spacing. c r1mach(5) = log10(b) c integer small(2) integer large(2) integer right(2) integer diver(2) integer log10(2) c needs to be (2) for autodouble, harris slash 6, ... integer sc save small, large, right, diver, log10, sc real rmach(5) equivalence (rmach(1),small(1)) equivalence (rmach(2),large(1)) equivalence (rmach(3),right(1)) equivalence (rmach(4),diver(1)) equivalence (rmach(5),log10(1)) integer j, k, l, t3e(3) data t3e(1) / 9777664 / data t3e(2) / 5323660 / data t3e(3) / 46980 / c this version adapts automatically to most current machines, c including auto-double compilers. c to compile on older machines, add a c in column 1 c on the next line data sc/0/ c and remove the c from column 1 in one of the sections below. c constants for even older machines can be obtained by c mail netlib@research.bell-labs.com c send old1mach from blas c please send corrections to dmg or ehg@bell-labs.com. c c machine constants for the honeywell dps 8/70 series. c data rmach(1) / o402400000000 / c data rmach(2) / o376777777777 / c data rmach(3) / o714400000000 / c data rmach(4) / o716400000000 / c data rmach(5) / o776464202324 /, sc/987/ c c machine constants for pdp-11 fortrans supporting c 32-bit integers (expressed in integer and octal). c data small(1) / 8388608 / c data large(1) / 2147483647 / c data right(1) / 880803840 / c data diver(1) / 889192448 / c data log10(1) / 1067065499 /, sc/987/ c data rmach(1) / o00040000000 / c data rmach(2) / o17777777777 / c data rmach(3) / o06440000000 / c data rmach(4) / o06500000000 / c data rmach(5) / o07746420233 /, sc/987/ c c machine constants for the univac 1100 series. c data rmach(1) / o000400000000 / c data rmach(2) / o377777777777 / c data rmach(3) / o146400000000 / c data rmach(4) / o147400000000 / c data rmach(5) / o177464202324 /, sc/987/ c if (sc .ne. 987) then * *** check for autodouble *** small(2) = 0 rmach(1) = 1e13 if (small(2) .ne. 0) then * *** autodoubled *** if ( small(1) .eq. 1117925532 * .and. small(2) .eq. -448790528) then * *** ieee big endian *** small(1) = 1048576 small(2) = 0 large(1) = 2146435071 large(2) = -1 right(1) = 1017118720 right(2) = 0 diver(1) = 1018167296 diver(2) = 0 log10(1) = 1070810131 log10(2) = 1352628735 else if ( small(2) .eq. 1117925532 * .and. small(1) .eq. -448790528) then * *** ieee little endian *** small(2) = 1048576 small(1) = 0 large(2) = 2146435071 large(1) = -1 right(2) = 1017118720 right(1) = 0 diver(2) = 1018167296 diver(1) = 0 log10(2) = 1070810131 log10(1) = 1352628735 else if ( small(1) .eq. -2065213935 * .and. small(2) .eq. 10752) then * *** vax with d_floating *** small(1) = 128 small(2) = 0 large(1) = -32769 large(2) = -1 right(1) = 9344 right(2) = 0 diver(1) = 9472 diver(2) = 0 log10(1) = 546979738 log10(2) = -805796613 else if ( small(1) .eq. 1267827943 * .and. small(2) .eq. 704643072) then * *** ibm mainframe *** small(1) = 1048576 small(2) = 0 large(1) = 2147483647 large(2) = -1 right(1) = 856686592 right(2) = 0 diver(1) = 873463808 diver(2) = 0 log10(1) = 1091781651 log10(2) = 1352628735 else write(*,9010) stop 777 end if else rmach(1) = 1234567. if (small(1) .eq. 1234613304) then * *** ieee *** small(1) = 8388608 large(1) = 2139095039 right(1) = 864026624 diver(1) = 872415232 log10(1) = 1050288283 else if (small(1) .eq. -1271379306) then * *** vax *** small(1) = 128 large(1) = -32769 right(1) = 13440 diver(1) = 13568 log10(1) = 547045274 else if (small(1) .eq. 1175639687) then * *** ibm mainframe *** small(1) = 1048576 large(1) = 2147483647 right(1) = 990904320 diver(1) = 1007681536 log10(1) = 1091781651 else if (small(1) .eq. 1251390520) then * *** convex c-1 *** small(1) = 8388608 large(1) = 2147483647 right(1) = 880803840 diver(1) = 889192448 log10(1) = 1067065499 else do 10 l = 1, 3 j = small(1) / 10000000 k = small(1) - 10000000*j if (k .ne. t3e(l)) go to 20 small(1) = j 10 continue * *** cray t3e *** call i1mcra(small, k, 16, 0, 0) call i1mcra(large, k, 32751, 16777215, 16777215) call i1mcra(right, k, 15520, 0, 0) call i1mcra(diver, k, 15536, 0, 0) call i1mcra(log10, k, 16339, 4461392, 10451455) go to 30 20 call i1mcra(j, k, 16405, 9876536, 0) if (small(1) .ne. j) then write(*,9020) stop 777 end if * *** cray 1, xmp, 2, and 3 *** call i1mcra(small(1), k, 8195, 8388608, 1) call i1mcra(large(1), k, 24574, 16777215, 16777214) call i1mcra(right(1), k, 16338, 8388608, 0) call i1mcra(diver(1), k, 16339, 8388608, 0) call i1mcra(log10(1), k, 16383, 10100890, 8715216) end if end if 30 sc = 987 end if * sanity check if (rmach(4) .ge. 1.0) stop 776 if (i .lt. 1 .or. i .gt. 5) then write(*,*) 'r1mach(i): i =',i,' is out of bounds.' stop end if r1mach = rmach(i) return 9010 format(/' adjust autodoubled r1mach by getting data'/ *' appropriate for your machine from d1mach.') 9020 format(/' adjust r1mach by uncommenting data statements'/ *' appropriate for your machine.') * /* c source for r1mach -- remove the * in column 1 */ *#include *#include *#include *float r1mach_(long *i) *{ * switch(*i){ * case 1: return flt_min; * case 2: return flt_max; * case 3: return flt_epsilon/flt_radix; * case 4: return flt_epsilon; * case 5: return log10((double)flt_radix); * } * fprintf(stderr, "invalid argument: r1mach(%ld)\n", *i); * exit(1); return 0; /* else complaint of missing return value */ *} end subroutine i1mcra(a, a1, b, c, d) **** special computation for cray machines **** integer a, a1, b, c, d a1 = 16777216*b + c a = 16777216*a1 + d end