c c # the following routines are taken from netlib c ---------------------------------------------- c subroutine dgefa(a,lda,n,ipvt,info) integer lda,n,ipvt(1),info double precision a(lda,1) 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 subroutine dgbfa(abd,lda,n,ml,mu,ipvt,info) integer lda,n,ml,mu,ipvt(1),info double precision abd(lda,1) 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 dgesl(a,lda,n,ipvt,b,job) integer lda,n,ipvt(1),job double precision a(lda,1),b(1) 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 dgbsl(abd,lda,n,ml,mu,ipvt,b,job) integer lda,n,ml,mu,ipvt(1),job double precision abd(lda,1),b(1) 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 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 double precision da,dx(1) integer i,incx,m,mp1,n,nincx c if(n.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 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 double precision dx(1),dy(1),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 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 double precision dx(1),dy(1),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 DOUBLE PRECISION FUNCTION D1MACH(I) C C DOUBLE-PRECISION MACHINE CONSTANTS C C D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. C C D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C C D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING. C C D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING. C C D1MACH( 5) = LOG10(B) C C TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT, C THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY C REMOVING THE C FROM COLUMN 1. C ON RARE MACHINES A STATIC STATEMENT MAY NEED TO BE ADDED. C (BUT PROBABLY MORE SYSTEMS PROHIBIT IT THAN REQUIRE IT.) C C FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), ONE OF THE FIRST C TWO SETS OF CONSTANTS BELOW SHOULD BE APPROPRIATE. C C WHERE POSSIBLE, DECIMAL, OCTAL OR HEXADECIMAL CONSTANTS ARE USED C TO SPECIFY THE CONSTANTS EXACTLY. SOMETIMES THIS REQUIRES USING C EQUIVALENT INTEGER ARRAYS. IF YOUR COMPILER USES HALF-WORD C INTEGERS BY DEFAULT (SOMETIMES CALLED INTEGER*2), YOU MAY NEED TO C CHANGE INTEGER TO INTEGER*4 OR OTHERWISE INSTRUCT YOUR COMPILER C TO USE FULL-WORD INTEGERS IN THE NEXT 5 DECLARATIONS. C INTEGER SMALL(4) INTEGER LARGE(4) INTEGER RIGHT(4) INTEGER DIVER(4) INTEGER LOG10(4) INTEGER SC C DOUBLE PRECISION DMACH(5) C 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 C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T C 3B SERIES AND MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T C PC 7300), IN WHICH THE MOST SIGNIFICANT BYTE IS STORED FIRST. C C DATA SMALL(1),SMALL(2) / 1048576, 0 / C DATA LARGE(1),LARGE(2) / 2146435071, -1 / C DATA RIGHT(1),RIGHT(2) / 1017118720, 0 / C DATA DIVER(1),DIVER(2) / 1018167296, 0 / C DATA LOG10(1),LOG10(2) / 1070810131, 1352628735 /, SC/987/ C C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES AND 8087-BASED C MICROS, SUCH AS THE IBM PC AND AT&T 6300, IN WHICH THE LEAST C SIGNIFICANT BYTE IS STORED FIRST. C C DATA SMALL(1),SMALL(2) / 0, 1048576 / C DATA LARGE(1),LARGE(2) / -1, 2146435071 / C DATA RIGHT(1),RIGHT(2) / 0, 1017118720 / C DATA DIVER(1),DIVER(2) / 0, 1018167296 / C DATA LOG10(1),LOG10(2) / 1352628735, 1070810131 /, SC/987/ C C MACHINE CONSTANTS FOR AMDAHL MACHINES. C C DATA SMALL(1),SMALL(2) / 1048576, 0 / C DATA LARGE(1),LARGE(2) / 2147483647, -1 / C DATA RIGHT(1),RIGHT(2) / 856686592, 0 / C DATA DIVER(1),DIVER(2) / 873463808, 0 / C DATA LOG10(1),LOG10(2) / 1091781651, 1352628735 /, SC/987/ C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA SMALL(1) / ZC00800000 / C DATA SMALL(2) / Z000000000 / C C DATA LARGE(1) / ZDFFFFFFFF / C DATA LARGE(2) / ZFFFFFFFFF / C C DATA RIGHT(1) / ZCC5800000 / C DATA RIGHT(2) / Z000000000 / C C DATA DIVER(1) / ZCC6800000 / C DATA DIVER(2) / Z000000000 / C C DATA LOG10(1) / ZD00E730E7 / C DATA LOG10(2) / ZC77800DC0 /, SC/987/ C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O0000000000000000 / C C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O0007777777777777 / C C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 /, SC/987/ C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O7770000000000000 / C C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O7777777777777777 / C C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 /, SC/987/ C C MACHINE CONSTANTS FOR FTN4 ON THE CDC 6000/7000 SERIES. C C DATA SMALL(1) / 00564000000000000000B / C DATA SMALL(2) / 00000000000000000000B / C C DATA LARGE(1) / 37757777777777777777B / C DATA LARGE(2) / 37157777777777777774B / C C DATA RIGHT(1) / 15624000000000000000B / C DATA RIGHT(2) / 00000000000000000000B / C C DATA DIVER(1) / 15634000000000000000B / C DATA DIVER(2) / 00000000000000000000B / C C DATA LOG10(1) / 17164642023241175717B / C DATA LOG10(2) / 16367571421742254654B /, SC/987/ C C MACHINE CONSTANTS FOR FTN5 ON THE CDC 6000/7000 SERIES. C C DATA SMALL(1) / O"00564000000000000000" / C DATA SMALL(2) / O"00000000000000000000" / C C DATA LARGE(1) / O"37757777777777777777" / C DATA LARGE(2) / O"37157777777777777774" / C C DATA RIGHT(1) / O"15624000000000000000" / C DATA RIGHT(2) / O"00000000000000000000" / C C DATA DIVER(1) / O"15634000000000000000" / C DATA DIVER(2) / O"00000000000000000000" / C C DATA LOG10(1) / O"17164642023241175717" / C DATA LOG10(2) / O"16367571421742254654" /, SC/987/ C C MACHINE CONSTANTS FOR CONVEX C-1 C C DATA SMALL(1),SMALL(2) / '00100000'X, '00000000'X / C DATA LARGE(1),LARGE(2) / '7FFFFFFF'X, 'FFFFFFFF'X / C DATA RIGHT(1),RIGHT(2) / '3CC00000'X, '00000000'X / C DATA DIVER(1),DIVER(2) / '3CD00000'X, '00000000'X / C DATA LOG10(1),LOG10(2) / '3FF34413'X, '509F79FF'X /, SC/987/ C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C C DATA SMALL(1) / 201354000000000000000B / C DATA SMALL(2) / 000000000000000000000B / C C DATA LARGE(1) / 577767777777777777777B / C DATA LARGE(2) / 000007777777777777776B / C C DATA RIGHT(1) / 376434000000000000000B / C DATA RIGHT(2) / 000000000000000000000B / C C DATA DIVER(1) / 376444000000000000000B / C DATA DIVER(2) / 000000000000000000000B / C C DATA LOG10(1) / 377774642023241175717B / C DATA LOG10(2) / 000007571421742254654B /, SC/987/ C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE - C STATIC DMACH(5) C C DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/ C DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/ C DATA LOG10/40423K,42023K,50237K,74776K/, SC/987/ C C MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7 C C DATA SMALL(1),SMALL(2) / '20000000, '00000201 / C DATA LARGE(1),LARGE(2) / '37777777, '37777577 / C DATA RIGHT(1),RIGHT(2) / '20000000, '00000333 / C DATA DIVER(1),DIVER(2) / '20000000, '00000334 / C DATA LOG10(1),LOG10(2) / '23210115, '10237777 /, SC/987/ C C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. C 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 THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86. C C DATA SMALL(1),SMALL(2) / Z00100000, Z00000000 / C DATA LARGE(1),LARGE(2) / Z7FFFFFFF, ZFFFFFFFF / C DATA RIGHT(1),RIGHT(2) / Z33100000, Z00000000 / C DATA DIVER(1),DIVER(2) / Z34100000, Z00000000 / C DATA LOG10(1),LOG10(2) / Z41134413, Z509F79FF /, SC/987/ C C MACHINE CONSTANTS FOR THE INTERDATA 8/32 C WITH THE UNIX SYSTEM FORTRAN 77 COMPILER. C C FOR THE INTERDATA FORTRAN VII COMPILER REPLACE C THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S. C C DATA SMALL(1),SMALL(2) / Z'00100000', Z'00000000' / C DATA LARGE(1),LARGE(2) / Z'7EFFFFFF', Z'FFFFFFFF' / C DATA RIGHT(1),RIGHT(2) / Z'33100000', Z'00000000' / C DATA DIVER(1),DIVER(2) / Z'34100000', Z'00000000' / C DATA LOG10(1),LOG10(2) / Z'41134413', Z'509F79FF' /, SC/987/ C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C DATA SMALL(1),SMALL(2) / "033400000000, "000000000000 / C DATA LARGE(1),LARGE(2) / "377777777777, "344777777777 / C DATA RIGHT(1),RIGHT(2) / "113400000000, "000000000000 / C DATA DIVER(1),DIVER(2) / "114400000000, "000000000000 / C DATA LOG10(1),LOG10(2) / "177464202324, "144117571776 /, SC/987/ C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C DATA SMALL(1),SMALL(2) / "000400000000, "000000000000 / C DATA LARGE(1),LARGE(2) / "377777777777, "377777777777 / C DATA RIGHT(1),RIGHT(2) / "103400000000, "000000000000 / C DATA DIVER(1),DIVER(2) / "104400000000, "000000000000 / C DATA LOG10(1),LOG10(2) / "177464202324, "047674776746 /, SC/987/ C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C 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 DATA SMALL(1),SMALL(2) / O00040000000, O00000000000 / C DATA LARGE(1),LARGE(2) / O17777777777, O37777777777 / C DATA RIGHT(1),RIGHT(2) / O04440000000, O00000000000 / C DATA DIVER(1),DIVER(2) / O04500000000, O00000000000 / C DATA LOG10(1),LOG10(2) / O07746420232, O20476747770 /, SC/987/ C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1),SMALL(2) / 128, 0 / C DATA SMALL(3),SMALL(4) / 0, 0 / C C DATA LARGE(1),LARGE(2) / 32767, -1 / C DATA LARGE(3),LARGE(4) / -1, -1 / C C DATA RIGHT(1),RIGHT(2) / 9344, 0 / C DATA RIGHT(3),RIGHT(4) / 0, 0 / C C DATA DIVER(1),DIVER(2) / 9472, 0 / C DATA DIVER(3),DIVER(4) / 0, 0 / C C DATA LOG10(1),LOG10(2) / 16282, 8346 / C DATA LOG10(3),LOG10(4) / -31493, -12296 /, SC/987/ C C DATA SMALL(1),SMALL(2) / O000200, O000000 / C DATA SMALL(3),SMALL(4) / O000000, O000000 / C C DATA LARGE(1),LARGE(2) / O077777, O177777 / C DATA LARGE(3),LARGE(4) / O177777, O177777 / C C DATA RIGHT(1),RIGHT(2) / O022200, O000000 / C DATA RIGHT(3),RIGHT(4) / O000000, O000000 / C C DATA DIVER(1),DIVER(2) / O022400, O000000 / C DATA DIVER(3),DIVER(4) / O000000, O000000 / C C DATA LOG10(1),LOG10(2) / O037632, O020232 / C DATA LOG10(3),LOG10(4) / O102373, O147770 /, SC/987/ C C MACHINE CONSTANTS FOR THE PRIME 50 SERIES SYSTEMS C WITH 32-BIT INTEGERS AND 64V MODE INSTRUCTIONS, C SUPPLIED BY IGOR BRAY. C C DATA SMALL(1),SMALL(2) / :10000000000, :00000100001 / C DATA LARGE(1),LARGE(2) / :17777777777, :37777677775 / C DATA RIGHT(1),RIGHT(2) / :10000000000, :00000000122 / C DATA DIVER(1),DIVER(2) / :10000000000, :00000000123 / C DATA LOG10(1),LOG10(2) / :11504046501, :07674600177 /, SC/987/ C C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000 C C DATA SMALL(1),SMALL(2) / $00000000, $00100000 / C DATA LARGE(1),LARGE(2) / $FFFFFFFF, $7FEFFFFF / C DATA RIGHT(1),RIGHT(2) / $00000000, $3CA00000 / C DATA DIVER(1),DIVER(2) / $00000000, $3CB00000 / C DATA LOG10(1),LOG10(2) / $509F79FF, $3FD34413 /, SC/987/ C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. C 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 MACHINE CONSTANTS FOR THE VAX UNIX F77 COMPILER C DATA SMALL(1),SMALL(2) / 128, 0 / DATA LARGE(1),LARGE(2) / -32769, -1 / DATA RIGHT(1),RIGHT(2) / 9344, 0 / DATA DIVER(1),DIVER(2) / 9472, 0 / DATA LOG10(1),LOG10(2) / 546979738, -805796613 /, SC/987/ C C MACHINE CONSTANTS FOR THE VAX-11 WITH C FORTRAN IV-PLUS COMPILER C C DATA SMALL(1),SMALL(2) / Z00000080, Z00000000 / C DATA LARGE(1),LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1),RIGHT(2) / Z00002480, Z00000000 / C DATA DIVER(1),DIVER(2) / Z00002500, Z00000000 / C DATA LOG10(1),LOG10(2) / Z209A3F9A, ZCFF884FB /, SC/987/ C C MACHINE CONSTANTS FOR VAX/VMS VERSION 2.2 C C DATA SMALL(1),SMALL(2) / '80'X, '0'X / C DATA LARGE(1),LARGE(2) / 'FFFF7FFF'X, 'FFFFFFFF'X / C DATA RIGHT(1),RIGHT(2) / '2480'X, '0'X / C DATA DIVER(1),DIVER(2) / '2500'X, '0'X / C DATA LOG10(1),LOG10(2) / '209A3F9A'X, 'CFF884FB'X /, SC/987/ C C *** ISSUE STOP 779 IF ALL DATA STATEMENTS ARE COMMENTED... IF (SC .NE. 987) STOP 779 IF (I .LT. 1 .OR. I .GT. 5) GOTO 999 D1MACH = DMACH(I) RETURN 999 WRITE(*,1999) I 1999 FORMAT(' D1MACH - I OUT OF BOUNDS',I10) STOP 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 double precision dx(1),dmax integer i,incx,ix,n c idamax = 0 if( n .lt. 1 ) 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 double precision function datanh (x) c june 1977 edition. w. fullerton, c3, los alamos scientific lab. double precision x, atnhcs(27), dxrel, sqeps, y, dcsevl, 1 d1mach, dlog, dsqrt external d1mach, dcsevl, dlog, dsqrt, initds c c series for atnh on the interval 0. to 2.50000e-01 c with weighted error 6.86e-32 c log weighted error 31.16 c significant figures required 30.00 c decimal places required 31.88 c data atnhcs( 1) / +.9439510239 3195492308 4289221863 3 d-1 / data atnhcs( 2) / +.4919843705 5786159472 0003457666 8 d-1 / data atnhcs( 3) / +.2102593522 4554327634 7932733175 2 d-2 / data atnhcs( 4) / +.1073554449 7761165846 4073104527 6 d-3 / data atnhcs( 5) / +.5978267249 2930314786 4278751787 2 d-5 / data atnhcs( 6) / +.3505062030 8891348459 6683488620 0 d-6 / data atnhcs( 7) / +.2126374343 7653403508 9621931443 1 d-7 / data atnhcs( 8) / +.1321694535 7155271921 2980172305 5 d-8 / data atnhcs( 9) / +.8365875501 1780703646 2360405295 9 d-10 / data atnhcs( 10) / +.5370503749 3110021638 8143458777 2 d-11 / data atnhcs( 11) / +.3486659470 1571079229 7124578429 0 d-12 / data atnhcs( 12) / +.2284549509 6034330155 2402411972 2 d-13 / data atnhcs( 13) / +.1508407105 9447930448 7422906755 8 d-14 / data atnhcs( 14) / +.1002418816 8041091261 3699572283 7 d-15 / data atnhcs( 15) / +.6698674738 1650695397 1552688298 6 d-17 / data atnhcs( 16) / +.4497954546 4949310830 8332762453 3 d-18 / data atnhcs( 17) / +.3032954474 2794535416 8236714666 6 d-19 / data atnhcs( 18) / +.2052702064 1909368264 6386141866 6 d-20 / data atnhcs( 19) / +.1393848977 0538377131 9301461333 3 d-21 / data atnhcs( 20) / +.9492580637 2245769719 5895466666 6 d-23 / data atnhcs( 21) / +.6481915448 2423076049 8244266666 6 d-24 / data atnhcs( 22) / +.4436730205 7236152726 3232000000 0 d-25 / data atnhcs( 23) / +.3043465618 5431616389 1200000000 0 d-26 / data atnhcs( 24) / +.2091881298 7923934740 4799999999 9 d-27 / data atnhcs( 25) / +.1440445411 2340505613 6533333333 3 d-28 / data atnhcs( 26) / +.9935374683 1416404650 6666666666 6 d-30 / data atnhcs( 27) / +.6863462444 3582600533 3333333333 3 d-31 / c data nterms, dxrel, sqeps / 0, 2*0.d0 / c if (nterms.ne.0) go to 10 nterms = initds (atnhcs, 27, 0.1*sngl(d1mach(3)) ) dxrel = dsqrt (d1mach(4)) sqeps = dsqrt (3.0d0*d1mach(3)) c 10 y = dabs(x) if (y.ge.1.d0) call seteru (20hdatanh dabs(x) ge 1, 20, 2, 2) c if (1.d0-y.lt.dxrel) call seteru ( 1 59hdatanh answer lt half precision because dabs(x) too near 1, 2 59, 1, 1) c datanh = x if (y.gt.sqeps .and. y.le.0.5d0) datanh = x*(1.0d0 + 1 dcsevl (8.d0*x*x-1.d0, atnhcs, nterms) ) if (y.gt.0.5d0) datanh = 0.5d0*dlog ((1.0d0+x)/(1.0d0-x) ) c return end function initds (dos, nos, eta) c june 1977 edition. w. fullerton, c3, los alamos scientific lab. c c initialize the double precision orthogonal series dos so that initds c is the number of terms needed to insure the error is no larger than c eta. ordinarily eta will be chosen to be one-tenth machine precision. c c input arguments -- c dos dble prec array of nos coefficients in an orthogonal series. c nos number of coefficients in dos. c eta requested accuracy of series. c double precision dos(nos) c if (nos.lt.1) call seteru ( 1 35hinitds number of coefficients lt 1, 35, 2, 2) c err = 0. do 10 ii=1,nos i = nos + 1 - ii err = err + abs(sngl(dos(i))) if (err.gt.eta) go to 20 10 continue c 20 if (i.eq.nos) call seteru (28hinitds eta may be too small, 28, 1 1, 2) initds = i c return end double precision function dsqrt (x) c june 1977 edition. w. fullerton, c3, los alamos scientific lab. double precision x, sqrt2(3), y, d9pak, d1mach external alog, d1mach, d9pak data sqrt2(1) / 0.7071067811 8654752440 0844362104 85 d0 / data sqrt2(2) / 1.0 d0 / data sqrt2(3) / 1.4142135623 7309504880 1688724209 70 d0 / c data niter / 0 / c if (niter.eq.0) niter = 1.443*alog(-0.104*alog(0.1*sngl(d1mach(3)) 1 )) + 1.0 c if (x.le.0.d0) go to 20 c call d9upak (x, y, n) ixpnt = n/2 irem = n - 2*ixpnt + 2 c c the approximation below has accuracy of 4.16 digits. z = y dsqrt = .261599e0 + z*(1.114292e0 + z*(-.516888e0 + z*.141067e0)) c do 10 iter=1,niter dsqrt = dsqrt + 0.5d0*(y - dsqrt*dsqrt) / dsqrt 10 continue c dsqrt = d9pak (sqrt2(irem)*dsqrt, ixpnt) return c 20 if (x.lt.0.d0) call seteru (21hdsqrt x is negative, 21, 1, 1) dsqrt = 0.0d0 return c end subroutine seteru (messg, nmessg, nerr, iopt) common /cseter/ iunflo integer messg(1) data iunflo / 0 / c if (iopt.ne.0) call seterr (messg, nmessg, nerr, iopt) if (iopt.ne.0) return c if (iunflo.le.0) return call seterr (messg, nmessg, nerr, 1) c return end double precision function dcsevl (x, a, n) c c evaluate the n-term chebyshev series a at x. adapted from c r. broucke, algorithm 446, c.a.c.m., 16, 254 (1973). c c input arguments -- c x dble prec value at which the series is to be evaluated. c a dble prec array of n terms of a chebyshev series. in eval- c uating a, only half the first coef is summed. c n number of terms in array a. c double precision a(n), x, twox, b0, b1, b2 c if (n.lt.1) call seteru (28hdcsevl number of terms le 0, 28, 2,2) if (n.gt.1000) call seteru (31hdcsevl number of terms gt 1000, 1 31, 3, 2) if (x.lt.(-1.1d0) .or. x.gt.1.1d0) call seteru ( 1 25hdcsevl x outside (-1,+1), 25, 1, 1) c twox = 2.0d0*x b1 = 0.d0 b0 = 0.d0 do 10 i=1,n b2 = b1 b1 = b0 ni = n - i + 1 b0 = twox*b1 - b2 + a(ni) 10 continue c dcsevl = 0.5d0 * (b0-b2) c return end double precision function dlog (x) c june 1977 edition. w. fullerton, c3, los alamos scientific lab. double precision x, alncs(11), center(4), alncen(5), aln2, y, t, 1 t2, xn, dcsevl, d1mach external d1mach, dcsevl, initds c c series for aln on the interval 0. to 3.46021e-03 c with weighted error 4.15e-32 c log weighted error 31.38 c significant figures required 31.21 c decimal places required 31.90 c data aln cs( 1) / +.1334719987 7973881561 6893860471 87 d+1 / data aln cs( 2) / +.6937562832 8411286281 3724383542 25 d-3 / data aln cs( 3) / +.4293403902 0450834506 5592108036 62 d-6 / data aln cs( 4) / +.2893384779 5432594580 4664403875 87 d-9 / data aln cs( 5) / +.2051251753 0340580901 7418134477 26 d-12 / data aln cs( 6) / +.1503971705 5497386574 6151533199 99 d-15 / data aln cs( 7) / +.1129454069 5636464284 5216133333 33 d-18 / data aln cs( 8) / +.8635578867 1171868881 9466666666 66 d-22 / data aln cs( 9) / +.6695299053 4350370613 3333333333 33 d-25 / data aln cs( 10) / +.5249155744 8151466666 6666666666 66 d-28 / data aln cs( 11) / +.4153054068 0362666666 6666666666 66 d-31 / c data center(1) / 1.0d0 / data center(2) / 1.25d0 / data center(3) / 1.50d0 / data center(4) / 1.75d0 / c data alncen( 1) / 0.0d0 / data alncen( 2) / +.2231435513 1420975576 6295090309 83 d+0 / data alncen( 3) / +.4054651081 0816438197 8013115464 34 d+0 / data alncen( 4) / +.5596157879 3542268627 0888500526 82 d+0 / data alncen( 5) / +.6931471805 5994530941 7232121458 17 d+0 / c c aln2 = alog(2.0) - 0.625 data aln2 / 0.0681471805 5994530941 7232121458 18d0 / data nterms / 0 / c if (nterms.eq.0) nterms = initds (alncs, 11, 28.9*sngl(d1mach(3))) c if (x.le.0.d0) call seteru ( 1 29hdlog x is zero or negative, 29, 1, 2) c call d9upak (x, y, n) c xn = n - 1 y = 2.0d0*y ntrval = 4.0d0*y - 2.5d0 c if (ntrval.eq.5) t = ((y-1.0d0)-1.0d0) / (y+2.0d0) if (ntrval.lt.5) t = (y-center(ntrval)) / (y+center(ntrval)) t2 = t*t dlog = 0.625d0*xn + (aln2*xn + alncen(ntrval) + 2.0d0*t + 1 t*t2*dcsevl(578.d0*t2-1.0d0, alncs, nterms) ) c return end function alog (x) c june 1977 edition. w. fullerton, c3, los alamos scientific lab. dimension alncs(6), center(4), alncen(5) external csevl, inits, r1mach c c series for aln on the interval 0. to 3.46021d-03 c with weighted error 1.50e-16 c log weighted error 15.82 c significant figures required 15.65 c decimal places required 16.21 c data aln cs( 1) / 1.3347199877 973882e0 / data aln cs( 2) / .0006937562 83284112e0 / data aln cs( 3) / .0000004293 40390204e0 / data aln cs( 4) / .0000000002 89338477e0 / data aln cs( 5) / .0000000000 00205125e0 / data aln cs( 6) / .0000000000 00000150e0 / c data center(1) / 1.0 / data center(2) / 1.25 / data center(3) / 1.50 / data center(4) / 1.75 / c data alncen( 1) / 0.0e0 / data alncen( 2) / +.2231435513 14209755 e+0 / data alncen( 3) / +.4054651081 08164381 e+0 / data alncen( 4) / +.5596157879 35422686 e+0 / data alncen( 5) / +.6931471805 59945309 e+0 / c c aln2 = alog(2.0) - 0.625 data aln2 / 0.0681471805 59945309e0 / data nterms / 0 / c if (nterms.eq.0) nterms = inits (alncs, 6, 28.9*r1mach(3)) c if (x.le.0.) call seteru ( 1 29halog x is zero or negative, 29, 1, 2) c call r9upak (x, y, n) c xn = n - 1 y = 2.0*y ntrval = 4.0*y - 2.5 if (ntrval.eq.5) t = ((y-1.0)-1.0) / (y+2.0) if (ntrval.lt.5) t = (y-center(ntrval))/(y+center(ntrval)) t2 = t*t c alog = 0.625*xn + (aln2*xn + alncen(ntrval) + 2.0*t + 1 t*t2*csevl(578.0*t2-1.0, alncs, nterms) ) c return end subroutine d9upak (x, y, n) c august 1980 portable edition. w. fullerton, los alamos scientific lab c c unpack floating point number x so that x = y * 2.0**n, where c 0.5 .le. abs(y) .lt. 1.0 . c double precision x, y, absx c absx = dabs(x) n = 0 y = 0.0d0 if (x.eq.0.0d0) return c 10 if (absx.ge.0.5d0) go to 20 n = n - 1 absx = absx*2.0d0 go to 10 c 20 if (absx.lt.1.0d0) go to 30 n = n + 1 absx = absx*0.5d0 go to 20 c 30 y = dsign (absx, x) return c end double precision function d9pak (y, n) c december 1979 edition. w. fullerton, c3, los alamos scientific lab. c c pack a base 2 exponent into floating point number x. this routine is c almost the inverse of d9upak. it is not exactly the inverse, because c dabs(x) need not be between 0.5 and 1.0. if both d9pak and 2.d0**n c were known to be in range we could compute c d9pak = x * 2.0d0**n c double precision y, aln2b, aln210, d1mach external d1mach, i1mach data nmin, nmax / 2*0 / data aln210 / 3.321928094 8873623478 7031942948 9 d0 / c if (nmin.ne.0) go to 10 aln2b = 1.0d0 if (i1mach(10).ne.2) aln2b = d1mach(5)*aln210 nmin = aln2b*dble(float(i1mach(15))) nmax = aln2b*dble(float(i1mach(16))) c 10 call d9upak (y, d9pak, ny) c nsum = n + ny if (nsum.lt.nmin) go to 40 if (nsum.gt.nmax) call seteru ( 1 31hd9pak packed number overflows, 31, 1, 2) c if (nsum.eq.0) return if (nsum.gt.0) go to 30 c 20 d9pak = 0.5d0*d9pak nsum = nsum + 1 if (nsum.ne.0) go to 20 return c 30 d9pak = 2.0d0*d9pak nsum = nsum - 1 if (nsum.ne.0) go to 30 return c 40 call seteru (32hd9pak packed number underflows, 32, 1, 0) d9pak = 0.0d0 return c end subroutine seterr (messg, nmessg, nerr, iopt) c c this version modified by w. fullerton to dump if iopt = 1 and c not recovering. c seterr sets lerror = nerr, optionally prints the message and dumps c according to the following rules... c c if iopt = 1 and recovering - just remember the error. c if iopt = 1 and not recovering - print, dump and stop. c if iopt = 2 - print, dump and stop. c c input c c messg - the error message. c nmessg - the length of the message, in characters. c nerr - the error number. must have nerr non-zero. c iopt - the option. must have iopt=1 or 2. c c error states - c c 1 - message length not positive. c 2 - cannot have nerr=0. c 3 - an unrecovered error followed by another error. c 4 - bad value for iopt. c c only the first 72 characters of the message are printed. c c the error handler calls a subroutine named fdump to produce a c symbolic dump. to complete the package, a dummy version of fdump c is supplied, but it should be replaced by a locally written version c which at least gives a trace-back. c integer messg(1) external i1mach, i8save c c the unit for error messages. c iwunit=i1mach(4) c if (nmessg.ge.1) go to 10 c c a message of non-positive length is fatal. c write(iwunit,9000) 9000 format(52h1error 1 in seterr - message length not positive.) go to 60 c c nw is the number of words the message occupies. c 10 nw=(min0(nmessg,72)-1)/i1mach(6)+1 c if (nerr.ne.0) go to 20 c c cannot turn the error state off using seterr. c write(iwunit,9001) 9001 format(42h1error 2 in seterr - cannot have nerr=0// 1 34h the current error message follows///) call e9rint(messg,nw,nerr,.true.) itemp=i8save(1,1,.true.) go to 50 c c set lerror and test for a previous unrecovered error. c 20 if (i8save(1,nerr,.true.).eq.0) go to 30 c write(iwunit,9002) 9002 format(23h1error 3 in seterr -, 1 48h an unrecovered error followed by another error.// 2 48h the previous and current error messages follow.///) call eprint call e9rint(messg,nw,nerr,.true.) go to 50 c c save this message in case it is not recovered from properly. c 30 call e9rint(messg,nw,nerr,.true.) c if (iopt.eq.1 .or. iopt.eq.2) go to 40 c c must have iopt = 1 or 2. c write(iwunit,9003) 9003 format(42h1error 4 in seterr - bad value for iopt// 1 34h the current error message follows///) go to 50 c c test for recovery. c 40 if (iopt.eq.2) go to 50 c if (i8save(2,0,.false.).eq.1) return c c call eprint c stop c 50 call eprint 60 call fdump stop c end function inits (os, nos, eta) c april 1977 version. w. fullerton, c3, los alamos scientific lab. c c initialize the orthogonal series so that inits is the number of terms c needed to insure the error is no larger than eta. ordinarily, eta c will be chosen to be one-tenth machine precision. c c input arguments -- c os array of nos coefficients in an orthogonal series. c nos number of coefficients in os. c eta requested accuracy of series. c dimension os(nos) c if (nos.lt.1) call seteru ( 1 35hinits number of coefficients lt 1, 35, 2, 2) c err = 0. do 10 ii=1,nos i = nos + 1 - ii err = err + abs(os(i)) if (err.gt.eta) go to 20 10 continue c 20 if (i.eq.nos) call seteru (28hinits eta may be too small, 28, 1 1, 2) inits = i c return end subroutine r9upak (x, y, n) c august 1980 portable edition. w. fullerton, los alamos scientific lab c c unpack floating point number x so that x = y * 2.0**n, where c 0.5 .le. abs(y) .lt. 1.0 . c absx = abs(x) n = 0 y = 0.0 if (x.eq.0.0) return c 10 if (absx.ge.0.5) go to 20 n = n - 1 absx = absx*2.0 go to 10 c 20 if (absx.lt.1.0) go to 30 n = n + 1 absx = absx*0.5 go to 20 c 30 y = sign (absx, x) return c end function csevl (x, cs, n) c april 1977 version. w. fullerton, c3, los alamos scientific lab. c c evaluate the n-term chebyshev series cs at x. adapted from c r. broucke, algorithm 446, c.a.c.m., 16, 254 (1973). also see fox c and parker, chebyshev polys in numerical analysis, oxford press, p.56. c c input arguments -- c x value at which the series is to be evaluated. c cs array of n terms of a chebyshev series. in eval- c uating cs, only half the first coef is summed. c n number of terms in array cs. c dimension cs(1) c if (n.lt.1) call seteru (28hcsevl number of terms le 0, 28, 2,2) if (n.gt.1000) call seteru (31hcsevl number of terms gt 1000, 1 31, 3, 2) if (x.lt.(-1.1) .or. x.gt.1.1) call seteru ( 1 25hcsevl x outside (-1,+1), 25, 1, 1) c b1 = 0. b0 = 0. twox = 2.*x do 10 i=1,n b2 = b1 b1 = b0 ni = n + 1 - i b0 = twox*b1 - b2 + cs(ni) 10 continue c csevl = 0.5 * (b0-b2) c return end subroutine e9rint(messg,nw,nerr,save) c c this routine stores the current error message or prints the old one, c if any, depending on whether or not save = .true. . c integer messg(nw) logical save external i1mach, i8save c c messgp stores at least the first 72 characters of the previous c message. its length is machine dependent and must be at least c c 1 + 71/(the number of characters stored per integer word). c integer messgp(36),fmt(14),ccplus c c start with no previous message. c data messgp(1)/1h1/, nwp/0/, nerrp/0/ c c set up the format for printing the error message. c the format is simply (a1,14x,72axx) where xx=i1mach(6) is the c number of characters stored per integer word. c data ccplus / 1h+ / c data fmt( 1) / 1h( / data fmt( 2) / 1ha / data fmt( 3) / 1h1 / data fmt( 4) / 1h, / data fmt( 5) / 1h1 / data fmt( 6) / 1h4 / data fmt( 7) / 1hx / data fmt( 8) / 1h, / data fmt( 9) / 1h7 / data fmt(10) / 1h2 / data fmt(11) / 1ha / data fmt(12) / 1hx / data fmt(13) / 1hx / data fmt(14) / 1h) / c if (.not.save) go to 20 c c save the message. c nwp=nw nerrp=nerr do 10 i=1,nw 10 messgp(i)=messg(i) c go to 30 c 20 if (i8save(1,0,.false.).eq.0) go to 30 c c print the message. c iwunit=i1mach(4) write(iwunit,9000) nerrp 9000 format(7h error ,i4,4h in ) c call s88fmt(2,i1mach(6),fmt(12)) write(iwunit,fmt) ccplus,(messgp(i),i=1,nwp) c 30 return c end integer function i8save(isw,ivalue,set) c c if (isw = 1) i8save returns the current error number and c sets it to ivalue if set = .true. . c c if (isw = 2) i8save returns the current recovery switch and c sets it to ivalue if set = .true. . c logical set c integer iparam(2) c iparam(1) is the error number and iparam(2) is the recovery switch. c c start execution error free and with recovery turned off. c data iparam(1) /0/, iparam(2) /2/ c i8save=iparam(isw) if (set) iparam(isw)=ivalue c return c end subroutine eprint c c this subroutine prints the last error message, if any. c integer messg(1) c call e9rint(messg,1,1,.false.) return c end subroutine s88fmt( n, w, ifmt ) c c s88fmt replaces ifmt(1), ... , ifmt(n) with c the characters corresponding to the n least significant c digits of w. c integer n,w,ifmt(n) c integer nt,wt,digits(10) c data digits( 1) / 1h0 / data digits( 2) / 1h1 / data digits( 3) / 1h2 / data digits( 4) / 1h3 / data digits( 5) / 1h4 / data digits( 6) / 1h5 / data digits( 7) / 1h6 / data digits( 8) / 1h7 / data digits( 9) / 1h8 / data digits(10) / 1h9 / c nt = n wt = w c 10 if (nt .le. 0) return idigit = mod( wt, 10 ) ifmt(nt) = digits(idigit+1) wt = wt/10 nt = nt - 1 go to 10 c end REAL FUNCTION R1MACH(I) C C SINGLE-PRECISION MACHINE CONSTANTS C C R1MACH(1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. C C R1MACH(2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C C R1MACH(3) = B**(-T), THE SMALLEST RELATIVE SPACING. C C R1MACH(4) = B**(1-T), THE LARGEST RELATIVE SPACING. C C R1MACH(5) = LOG10(B) C C TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT, C THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY C REMOVING THE C FROM COLUMN 1. C ON RARE MACHINES A STATIC STATEMENT MAY NEED TO BE ADDED. C (BUT PROBABLY MORE SYSTEMS PROHIBIT IT THAN REQUIRE IT.) C C FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), THE FIRST C SET OF CONSTANTS BELOW SHOULD BE APPROPRIATE. C C WHERE POSSIBLE, DECIMAL, OCTAL OR HEXADECIMAL CONSTANTS ARE USED C TO SPECIFY THE CONSTANTS EXACTLY. SOMETIMES THIS REQUIRES USING C EQUIVALENT INTEGER ARRAYS. IF YOUR COMPILER USES HALF-WORD C INTEGERS BY DEFAULT (SOMETIMES CALLED INTEGER*2), YOU MAY NEED TO C CHANGE INTEGER TO INTEGER*4 OR OTHERWISE INSTRUCT YOUR COMPILER C TO USE FULL-WORD INTEGERS IN THE NEXT 5 DECLARATIONS. C INTEGER SMALL(2) INTEGER LARGE(2) INTEGER RIGHT(2) INTEGER DIVER(2) INTEGER LOG10(2) INTEGER SC C REAL RMACH(5) C 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)) C C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T C 3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T C PC 7300), AND 8087 BASED MICROS (E.G. IBM PC AND AT&T 6300). C C DATA SMALL(1) / 8388608 / C DATA LARGE(1) / 2139095039 / C DATA RIGHT(1) / 864026624 / C DATA DIVER(1) / 872415232 / C DATA LOG10(1) / 1050288283 /, SC/987/ C C MACHINE CONSTANTS FOR AMDAHL MACHINES. C C DATA SMALL(1) / 1048576 / C DATA LARGE(1) / 2147483647 / C DATA RIGHT(1) / 990904320 / C DATA DIVER(1) / 1007681536 / C DATA LOG10(1) / 1091781651 /, SC/987/ C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA RMACH(1) / Z400800000 / C DATA RMACH(2) / Z5FFFFFFFF / C DATA RMACH(3) / Z4E9800000 / C DATA RMACH(4) / Z4EA800000 / C DATA RMACH(5) / Z500E730E8 /, SC/987/ C C MACHINE CONSTANTS FOR THE BURROUGHS 5700/6700/7700 SYSTEMS. C C DATA RMACH(1) / O1771000000000000 / C DATA RMACH(2) / O0777777777777777 / C DATA RMACH(3) / O1311000000000000 / C DATA RMACH(4) / O1301000000000000 / C DATA RMACH(5) / O1157163034761675 /, SC/987/ C C MACHINE CONSTANTS FOR FTN4 ON THE CDC 6000/7000 SERIES. C C DATA RMACH(1) / 00564000000000000000B / C DATA RMACH(2) / 37767777777777777776B / C DATA RMACH(3) / 16414000000000000000B / C DATA RMACH(4) / 16424000000000000000B / C DATA RMACH(5) / 17164642023241175720B /, SC/987/ C C MACHINE CONSTANTS FOR FTN5 ON THE CDC 6000/7000 SERIES. C C DATA RMACH(1) / O"00564000000000000000" / C DATA RMACH(2) / O"37767777777777777776" / C DATA RMACH(3) / O"16414000000000000000" / C DATA RMACH(4) / O"16424000000000000000" / C DATA RMACH(5) / O"17164642023241175720" /, SC/987/ C C MACHINE CONSTANTS FOR CONVEX C-1. C C DATA RMACH(1) / '00800000'X / C DATA RMACH(2) / '7FFFFFFF'X / C DATA RMACH(3) / '34800000'X / C DATA RMACH(4) / '35000000'X / C DATA RMACH(5) / '3F9A209B'X /, SC/987/ C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C C DATA RMACH(1) / 200034000000000000000B / C DATA RMACH(2) / 577767777777777777776B / C DATA RMACH(3) / 377224000000000000000B / C DATA RMACH(4) / 377234000000000000000B / C DATA RMACH(5) / 377774642023241175720B /, SC/987/ C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200. C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE - C STATIC RMACH(5) C C DATA SMALL/20K,0/,LARGE/77777K,177777K/ C DATA RIGHT/35420K,0/,DIVER/36020K,0/ C DATA LOG10/40423K,42023K/, SC/987/ C C MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7. C C DATA SMALL(1),SMALL(2) / '20000000, '00000201 / C DATA LARGE(1),LARGE(2) / '37777777, '00000177 / C DATA RIGHT(1),RIGHT(2) / '20000000, '00000352 / C DATA DIVER(1),DIVER(2) / '20000000, '00000353 / C DATA LOG10(1),LOG10(2) / '23210115, '00000377 /, SC/987/ C C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. C 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 THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86. C C DATA RMACH(1) / Z00100000 / C DATA RMACH(2) / Z7FFFFFFF / C DATA RMACH(3) / Z3B100000 / C DATA RMACH(4) / Z3C100000 / C DATA RMACH(5) / Z41134413 /, SC/987/ C C MACHINE CONSTANTS FOR THE INTERDATA 8/32 C WITH THE UNIX SYSTEM FORTRAN 77 COMPILER. C C FOR THE INTERDATA FORTRAN VII COMPILER REPLACE C THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S. C C DATA RMACH(1) / Z'00100000' / C DATA RMACH(2) / Z'7EFFFFFF' / C DATA RMACH(3) / Z'3B100000' / C DATA RMACH(4) / Z'3C100000' / C DATA RMACH(5) / Z'41134413' /, SC/987/ C C MACHINE CONSTANTS FOR THE PDP-10 (KA OR KI PROCESSOR). C C DATA RMACH(1) / "000400000000 / C DATA RMACH(2) / "377777777777 / C DATA RMACH(3) / "146400000000 / C DATA RMACH(4) / "147400000000 / C DATA RMACH(5) / "177464202324 /, SC/987/ C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C 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 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 PDP-11 FORTRANS SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1),SMALL(2) / 128, 0 / C DATA LARGE(1),LARGE(2) / 32767, -1 / C DATA RIGHT(1),RIGHT(2) / 13440, 0 / C DATA DIVER(1),DIVER(2) / 13568, 0 / C DATA LOG10(1),LOG10(2) / 16282, 8347 /, SC/987/ C C DATA SMALL(1),SMALL(2) / O000200, O000000 / C DATA LARGE(1),LARGE(2) / O077777, O177777 / C DATA RIGHT(1),RIGHT(2) / O032200, O000000 / C DATA DIVER(1),DIVER(2) / O032400, O000000 / C DATA LOG10(1),LOG10(2) / O037632, O020233 /, SC/987/ C C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000. C C DATA SMALL(1) / $00800000 / C DATA LARGE(1) / $7F7FFFFF / C DATA RIGHT(1) / $33800000 / C DATA DIVER(1) / $34000000 / C DATA LOG10(1) / $3E9A209B /, SC/987/ C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. C 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 C MACHINE CONSTANTS FOR THE VAX UNIX F77 COMPILER. C DATA SMALL(1) / 128 / DATA LARGE(1) / -32769 / DATA RIGHT(1) / 13440 / DATA DIVER(1) / 13568 / DATA LOG10(1) / 547045274 /, SC/987/ C C MACHINE CONSTANTS FOR THE VAX-11 WITH C FORTRAN IV-PLUS COMPILER. C C DATA RMACH(1) / Z00000080 / C DATA RMACH(2) / ZFFFF7FFF / C DATA RMACH(3) / Z00003480 / C DATA RMACH(4) / Z00003500 / C DATA RMACH(5) / Z209B3F9A /, SC/987/ C C MACHINE CONSTANTS FOR VAX/VMS VERSION 2.2. C C DATA RMACH(1) / '80'X / C DATA RMACH(2) / 'FFFF7FFF'X / C DATA RMACH(3) / '3480'X / C DATA RMACH(4) / '3500'X / C DATA RMACH(5) / '209B3F9A'X /, SC/987/ C C *** ISSUE STOP 778 IF ALL DATA STATEMENTS ARE COMMENTED... IF (SC .NE. 987) STOP 778 IF (I .LT. 1 .OR. I .GT. 5) GOTO 999 R1MACH = RMACH(I) RETURN 999 WRITE(*,1999) I 1999 FORMAT(' R1MACH - I OUT OF BOUNDS',I10) STOP END INTEGER FUNCTION I1MACH(I) C C I/O UNIT NUMBERS. C C I1MACH( 1) = THE STANDARD INPUT UNIT. C C I1MACH( 2) = THE STANDARD OUTPUT UNIT. C C I1MACH( 3) = THE STANDARD PUNCH UNIT. C C I1MACH( 4) = THE STANDARD ERROR MESSAGE UNIT. C C WORDS. C C I1MACH( 5) = THE NUMBER OF BITS PER INTEGER STORAGE UNIT. C C I1MACH( 6) = THE NUMBER OF CHARACTERS PER CHARACTER STORAGE UNIT. C FOR FORTRAN 77, THIS IS ALWAYS 1. FOR FORTRAN 66, C CHARACTER STORAGE UNIT = INTEGER STORAGE UNIT. C C INTEGERS. C C ASSUME INTEGERS ARE REPRESENTED IN THE S-DIGIT, BASE-A FORM C C SIGN ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) ) C C WHERE 0 .LE. X(I) .LT. A FOR I=0,...,S-1. C C I1MACH( 7) = A, THE BASE. C C I1MACH( 8) = S, THE NUMBER OF BASE-A DIGITS. C C I1MACH( 9) = A**S - 1, THE LARGEST MAGNITUDE. C C FLOATING-POINT NUMBERS. C C ASSUME FLOATING-POINT NUMBERS ARE REPRESENTED IN THE T-DIGIT, C BASE-B FORM C C SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) ) C C WHERE 0 .LE. X(I) .LT. B FOR I=1,...,T, C 0 .LT. X(1), AND EMIN .LE. E .LE. EMAX. C C I1MACH(10) = B, THE BASE. C C SINGLE-PRECISION C C I1MACH(11) = T, THE NUMBER OF BASE-B DIGITS. C C I1MACH(12) = EMIN, THE SMALLEST EXPONENT E. C C I1MACH(13) = EMAX, THE LARGEST EXPONENT E. C C DOUBLE-PRECISION C C I1MACH(14) = T, THE NUMBER OF BASE-B DIGITS. C C I1MACH(15) = EMIN, THE SMALLEST EXPONENT E. C C I1MACH(16) = EMAX, THE LARGEST EXPONENT E. C C TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT, C THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY C REMOVING THE C FROM COLUMN 1. ALSO, THE VALUES OF C I1MACH(1) - I1MACH(4) SHOULD BE CHECKED FOR CONSISTENCY C WITH THE LOCAL OPERATING SYSTEM. FOR FORTRAN 77, YOU MAY WISH C TO ADJUST THE DATA STATEMENT SO IMACH(6) IS SET TO 1, AND C THEN TO COMMENT OUT THE EXECUTABLE TEST ON I .EQ. 6 BELOW. C ON RARE MACHINES A STATIC STATEMENT MAY NEED TO BE ADDED. C (BUT PROBABLY MORE SYSTEMS PROHIBIT IT THAN REQUIRE IT.) C C FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), THE FIRST C SET OF CONSTANTS BELOW SHOULD BE APPROPRIATE, EXCEPT PERHAPS C FOR IMACH(1) - IMACH(4). C INTEGER IMACH(16),OUTPUT,SANITY C EQUIVALENCE (IMACH(4),OUTPUT) C C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T C 3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T C PC 7300), AND 8087 BASED MICROS (E.G. IBM PC AND AT&T 6300). C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 128 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 /, SANITY/987/ C C MACHINE CONSTANTS FOR AMDAHL MACHINES. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 /, SANITY/987/ C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA IMACH( 1) / 7 / C DATA IMACH( 2) / 2 / C DATA IMACH( 3) / 2 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 33 / C DATA IMACH( 9) / Z1FFFFFFFF / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -256 / C DATA IMACH(13) / 255 / C DATA IMACH(14) / 60 / C DATA IMACH(15) / -256 / C DATA IMACH(16) / 255 /, SANITY/987/ C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 39 / C DATA IMACH( 9) / O0007777777777777 / C DATA IMACH(10) / 8 / C DATA IMACH(11) / 13 / C DATA IMACH(12) / -50 / C DATA IMACH(13) / 76 / C DATA IMACH(14) / 26 / C DATA IMACH(15) / -50 / C DATA IMACH(16) / 76 /, SANITY/987/ C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 39 / C DATA IMACH( 9) / O0007777777777777 / C DATA IMACH(10) / 8 / C DATA IMACH(11) / 13 / C DATA IMACH(12) / -50 / C DATA IMACH(13) / 76 / C DATA IMACH(14) / 26 / C DATA IMACH(15) / -32754 / C DATA IMACH(16) / 32780 /, SANITY/987/ C C MACHINE CONSTANTS FOR FTN4 ON THE CDC 6000/7000 SERIES. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 60 / C DATA IMACH( 6) / 10 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 48 / C DATA IMACH( 9) / 00007777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -929 / C DATA IMACH(13) / 1070 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -929 / C DATA IMACH(16) / 1069 /, SANITY/987/ C C MACHINE CONSTANTS FOR FTN5 ON THE CDC 6000/7000 SERIES. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 60 / C DATA IMACH( 6) / 10 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 48 / C DATA IMACH( 9) / O"00007777777777777777" / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -929 / C DATA IMACH(13) / 1070 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -929 / C DATA IMACH(16) / 1069 /, SANITY/987/ C C MACHINE CONSTANTS FOR CONVEX C-1. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) /-1024 / C DATA IMACH(16) / 1023 /, SANITY/987/ C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 102 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 63 / C DATA IMACH( 9) / 777777777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -8189 / C DATA IMACH(13) / 8190 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -8099 / C DATA IMACH(16) / 8190 /, SANITY/987/ C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200. C C DATA IMACH( 1) / 11 / C DATA IMACH( 2) / 12 / C DATA IMACH( 3) / 8 / C DATA IMACH( 4) / 10 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) /32767 / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 /, SANITY/987/ C C MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / 3 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 23 / C DATA IMACH( 9) / 8388607 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 23 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 38 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 /, SANITY/987/ C C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 43 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / O377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 63 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 /, SANITY/987/ C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / Z7FFFFFFF / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 /, SANITY/987/ C C MACHINE CONSTANTS FOR THE INTERDATA 8/32 C WITH THE UNIX SYSTEM FORTRAN 77 COMPILER. C C FOR THE INTERDATA FORTRAN VII COMPILER REPLACE C THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / Z'7FFFFFFF' / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 62 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 62 /, SANITY/987/ C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 5 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / "377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 54 / C DATA IMACH(15) / -101 / C DATA IMACH(16) / 127 /, SANITY/987/ C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 5 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / "377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 62 / C DATA IMACH(15) / -128 / C DATA IMACH(16) / 127 /, SANITY/987/ C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 32-BIT INTEGER ARITHMETIC. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 /, SANITY/987/ C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 16-BIT INTEGER ARITHMETIC. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) / 32767 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 /, SANITY/987/ C C MACHINE CONSTANTS FOR THE PRIME 50 SERIES SYSTEMS C WTIH 32-BIT INTEGERS AND 64V MODE INSTRUCTIONS, C SUPPLIED BY IGOR BRAY. C C DATA IMACH( 1) / 1 / C DATA IMACH( 2) / 1 / C DATA IMACH( 3) / 2 / C DATA IMACH( 4) / 1 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / :17777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 23 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / +127 / C DATA IMACH(14) / 47 / C DATA IMACH(15) / -32895 / C DATA IMACH(16) / +32637 /, SANITY/987/ C C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000. C C DATA IMACH( 1) / 0 / C DATA IMACH( 2) / 0 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 0 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 1 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 128 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 /, SANITY/987/ C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. C C NOTE THAT THE PUNCH UNIT, I1MACH(3), HAS BEEN SET TO 7 C WHICH IS APPROPRIATE FOR THE UNIVAC-FOR SYSTEM. C IF YOU HAVE THE UNIVAC-FTN SYSTEM, SET IT TO 1. C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / O377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 60 / C DATA IMACH(15) /-1024 / C DATA IMACH(16) / 1023 /, SANITY/987/ C C MACHINE CONSTANTS FOR VAX. C DATA IMACH( 1) / 5 / DATA IMACH( 2) / 6 / DATA IMACH( 3) / 7 / DATA IMACH( 4) / 6 / DATA IMACH( 5) / 32 / DATA IMACH( 6) / 4 / DATA IMACH( 7) / 2 / DATA IMACH( 8) / 31 / DATA IMACH( 9) / 2147483647 / DATA IMACH(10) / 2 / DATA IMACH(11) / 24 / DATA IMACH(12) / -127 / DATA IMACH(13) / 127 / DATA IMACH(14) / 56 / DATA IMACH(15) / -127 / DATA IMACH(16) / 127 /, SANITY/987/ C C *** ISSUE STOP 777 IF ALL DATA STATEMENTS ARE COMMENTED... IF (SANITY .NE. 987) STOP 777 IF (I .LT. 1 .OR. I .GT. 16) GO TO 10 C I1MACH = IMACH(I) C/6S C/7S IF(I.EQ.6) I1MACH=1 C/ RETURN 10 WRITE(OUTPUT,1999) I 1999 FORMAT(' I1MACH - I OUT OF BOUNDS',I10) STOP END subroutine fdump return end subroutine stoda (neq, y, yh, nyh, yh1, ewt, savf, acor, 1 wm, iwm, f, jac, pjac, slvs) clll. optimize external f, jac, pjac, slvs integer neq, nyh, iwm integer iownd, ialth, ipup, lmax, meo, nqnyh, nslp, 1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, 2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu integer iownd2, icount, irflag, jtyp, mused, mxordn, mxords integer i, i1, iredo, iret, j, jb, m, ncf, newq integer lm1, lm1p1, lm2, lm2p1, nqm1, nqm2 double precision y, yh, yh1, ewt, savf, acor, wm double precision conit, crate, el, elco, hold, rmax, tesco, 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround double precision rownd2, pdest, pdlast, ratio, cm1, cm2, 1 pdnorm double precision dcon, ddn, del, delp, dsm, dup, exdn, exsm, exup, 1 r, rh, rhdn, rhsm, rhup, told, vmnorm double precision alpha, dm1, dm2, exm1, exm2, pdh, pnorm, rate, 1 rh1, rh1it, rh2, rm, sm1 dimension neq(1), y(1), yh(nyh,1), yh1(1), ewt(1), savf(1), 1 acor(1), wm(1), iwm(1) dimension sm1(12) common /ls0001/ conit, crate, el(13), elco(13,12), 1 hold, rmax, tesco(3,12), 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, iownd(14), 3 ialth, ipup, lmax, meo, nqnyh, nslp, 4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu common /lsa001/ rownd2, pdest, pdlast, ratio, cm1(12), cm2(5), 1 pdnorm, 2 iownd2(3), icount, irflag, jtyp, mused, mxordn, mxords data sm1/0.5d0, 0.575d0, 0.55d0, 0.45d0, 0.35d0, 0.25d0, 1 0.20d0, 0.15d0, 0.10d0, 0.075d0, 0.050d0, 0.025d0/ c----------------------------------------------------------------------- c stoda performs one step of the integration of an initial value c problem for a system of ordinary differential equations. c note.. stoda is independent of the value of the iteration method c indicator miter, when this is .ne. 0, and hence is independent c of the type of chord method used, or the jacobian structure. c communication with stoda is done with the following variables.. c c y = an array of length .ge. n used as the y argument in c all calls to f and jac. c neq = integer array containing problem size in neq(1), and c passed as the neq argument in all calls to f and jac. c yh = an nyh by lmax array containing the dependent variables c and their approximate scaled derivatives, where c lmax = maxord + 1. yh(i,j+1) contains the approximate c j-th derivative of y(i), scaled by h**j/factorial(j) c (j = 0,1,...,nq). on entry for the first step, the first c two columns of yh must be set from the initial values. c nyh = a constant integer .ge. n, the first dimension of yh. c yh1 = a one-dimensional array occupying the same space as yh. c ewt = an array of length n containing multiplicative weights c for local error measurements. local errors in y(i) are c compared to 1.0/ewt(i) in various error tests. c savf = an array of working storage, of length n. c acor = a work array of length n, used for the accumulated c corrections. on a successful return, acor(i) contains c the estimated one-step local error in y(i). c wm,iwm = real and integer work arrays associated with matrix c operations in chord iteration (miter .ne. 0). c pjac = name of routine to evaluate and preprocess jacobian matrix c and p = i - h*el0*jac, if a chord method is being used. c it also returns an estimate of norm(jac) in pdnorm. c slvs = name of routine to solve linear system in chord iteration. c ccmax = maximum relative change in h*el0 before pjac is called. c h = the step size to be attempted on the next step. c h is altered by the error control algorithm during the c problem. h can be either positive or negative, but its c sign must remain constant throughout the problem. c hmin = the minimum absolute value of the step size h to be used. c hmxi = inverse of the maximum absolute value of h to be used. c hmxi = 0.0 is allowed and corresponds to an infinite hmax. c hmin and hmxi may be changed at any time, but will not c take effect until the next change of h is considered. c tn = the independent variable. tn is updated on each step taken. c jstart = an integer used for input only, with the following c values and meanings.. c 0 perform the first step. c .gt.0 take a new step continuing from the last. c -1 take the next step with a new value of h, c n, meth, miter, and/or matrix parameters. c -2 take the next step with a new value of h, c but with other inputs unchanged. c on return, jstart is set to 1 to facilitate continuation. c kflag = a completion code with the following meanings.. c 0 the step was succesful. c -1 the requested error could not be achieved. c -2 corrector convergence could not be achieved. c -3 fatal error in pjac or slvs. c a return with kflag = -1 or -2 means either c abs(h) = hmin or 10 consecutive failures occurred. c on a return with kflag negative, the values of tn and c the yh array are as of the beginning of the last c step, and h is the last step size attempted. c maxord = the maximum order of integration method to be allowed. c maxcor = the maximum number of corrector iterations allowed. c msbp = maximum number of steps between pjac calls (miter .gt. 0). c mxncf = maximum number of convergence failures allowed. c meth = current method. c meth = 1 means adams method (nonstiff) c meth = 2 means bdf method (stiff) c meth may be reset by stoda. c miter = corrector iteration method. c miter = 0 means functional iteration. c miter = jt .gt. 0 means a chord iteration corresponding c to jacobian type jt. (the lsoda argument jt is c communicated here as jtyp, but is not used in stoda c except to load miter following a method switch.) c miter may be reset by stoda. c n = the number of first-order differential equations. c----------------------------------------------------------------------- kflag = 0 told = tn ncf = 0 ierpj = 0 iersl = 0 jcur = 0 icf = 0 delp = 0.0d0 if (jstart .gt. 0) go to 200 if (jstart .eq. -1) go to 100 if (jstart .eq. -2) go to 160 c----------------------------------------------------------------------- c on the first call, the order is set to 1, and other variables are c initialized. rmax is the maximum ratio by which h can be increased c in a single step. it is initially 1.e4 to compensate for the small c initial h, but then is normally equal to 10. if a failure c occurs (in corrector convergence or error test), rmax is set at 2 c for the next increase. c cfode is called to get the needed coefficients for both methods. c----------------------------------------------------------------------- lmax = maxord + 1 nq = 1 l = 2 ialth = 2 rmax = 10000.0d0 rc = 0.0d0 el0 = 1.0d0 crate = 0.7d0 hold = h nslp = 0 ipup = miter iret = 3 c initialize switching parameters. meth = 1 is assumed initially. ----- icount = 20 irflag = 0 pdest = 0.0d0 pdlast = 0.0d0 ratio = 5.0d0 call cfode (2, elco, tesco) do 10 i = 1,5 10 cm2(i) = tesco(2,i)*elco(i+1,i) call cfode (1, elco, tesco) do 20 i = 1,12 20 cm1(i) = tesco(2,i)*elco(i+1,i) go to 150 c----------------------------------------------------------------------- c the following block handles preliminaries needed when jstart = -1. c ipup is set to miter to force a matrix update. c if an order increase is about to be considered (ialth = 1), c ialth is reset to 2 to postpone consideration one more step. c if the caller has changed meth, cfode is called to reset c the coefficients of the method. c if h is to be changed, yh must be rescaled. c if h or meth is being changed, ialth is reset to l = nq + 1 c to prevent further changes in h for that many steps. c----------------------------------------------------------------------- 100 ipup = miter lmax = maxord + 1 if (ialth .eq. 1) ialth = 2 if (meth .eq. mused) go to 160 call cfode (meth, elco, tesco) ialth = l iret = 1 c----------------------------------------------------------------------- c the el vector and related constants are reset c whenever the order nq is changed, or at the start of the problem. c----------------------------------------------------------------------- 150 do 155 i = 1,l 155 el(i) = elco(i,nq) nqnyh = nq*nyh rc = rc*el(1)/el0 el0 = el(1) conit = 0.5d0/dfloat(nq+2) go to (160, 170, 200), iret c----------------------------------------------------------------------- c if h is being changed, the h ratio rh is checked against c rmax, hmin, and hmxi, and the yh array rescaled. ialth is set to c l = nq + 1 to prevent a change of h for that many steps, unless c forced by a convergence or error test failure. c----------------------------------------------------------------------- 160 if (h .eq. hold) go to 200 rh = h/hold h = hold iredo = 3 go to 175 170 rh = dmax1(rh,hmin/dabs(h)) 175 rh = dmin1(rh,rmax) rh = rh/dmax1(1.0d0,dabs(h)*hmxi*rh) c----------------------------------------------------------------------- c if meth = 1, also restrict the new step size by the stability region. c if this reduces h, set irflag to 1 so that if there are roundoff c problems later, we can assume that is the cause of the trouble. c----------------------------------------------------------------------- if (meth .eq. 2) go to 178 irflag = 0 pdh = dmax1(dabs(h)*pdlast,0.000001d0) if (rh*pdh*1.00001d0 .lt. sm1(nq)) go to 178 rh = sm1(nq)/pdh irflag = 1 178 continue r = 1.0d0 do 180 j = 2,l r = r*rh do 180 i = 1,n 180 yh(i,j) = yh(i,j)*r h = h*rh rc = rc*rh ialth = l if (iredo .eq. 0) go to 690 c----------------------------------------------------------------------- c this section computes the predicted values by effectively c multiplying the yh array by the pascal triangle matrix. c rc is the ratio of new to old values of the coefficient h*el(1). c when rc differs from 1 by more than ccmax, ipup is set to miter c to force pjac to be called, if a jacobian is involved. c in any case, pjac is called at least every msbp steps. c----------------------------------------------------------------------- 200 if (dabs(rc-1.0d0) .gt. ccmax) ipup = miter if (nst .ge. nslp+msbp) ipup = miter tn = tn + h i1 = nqnyh + 1 do 215 jb = 1,nq i1 = i1 - nyh cdir$ ivdep do 210 i = i1,nqnyh 210 yh1(i) = yh1(i) + yh1(i+nyh) 215 continue pnorm = vmnorm (n, yh1, ewt) c----------------------------------------------------------------------- c up to maxcor corrector iterations are taken. a convergence test is c made on the r.m.s. norm of each correction, weighted by the error c weight vector ewt. the sum of the corrections is accumulated in the c vector acor(i). the yh array is not altered in the corrector loop. c----------------------------------------------------------------------- 220 m = 0 rate = 0.0d0 del = 0.0d0 do 230 i = 1,n 230 y(i) = yh(i,1) call f (neq, tn, y, savf) nfe = nfe + 1 if (ipup .le. 0) go to 250 c----------------------------------------------------------------------- c if indicated, the matrix p = i - h*el(1)*j is reevaluated and c preprocessed before starting the corrector iteration. ipup is set c to 0 as an indicator that this has been done. c----------------------------------------------------------------------- call pjac (neq, y, yh, nyh, ewt, acor, savf, wm, iwm, f, jac) ipup = 0 rc = 1.0d0 nslp = nst crate = 0.7d0 if (ierpj .ne. 0) go to 430 250 do 260 i = 1,n 260 acor(i) = 0.0d0 270 if (miter .ne. 0) go to 350 c----------------------------------------------------------------------- c in the case of functional iteration, update y directly from c the result of the last function evaluation. c----------------------------------------------------------------------- do 290 i = 1,n savf(i) = h*savf(i) - yh(i,2) 290 y(i) = savf(i) - acor(i) del = vmnorm (n, y, ewt) do 300 i = 1,n y(i) = yh(i,1) + el(1)*savf(i) 300 acor(i) = savf(i) go to 400 c----------------------------------------------------------------------- c in the case of the chord method, compute the corrector error, c and solve the linear system with that as right-hand side and c p as coefficient matrix. c----------------------------------------------------------------------- 350 do 360 i = 1,n 360 y(i) = h*savf(i) - (yh(i,2) + acor(i)) call slvs (wm, iwm, y, savf) if (iersl .lt. 0) go to 430 if (iersl .gt. 0) go to 410 del = vmnorm (n, y, ewt) do 380 i = 1,n acor(i) = acor(i) + y(i) 380 y(i) = yh(i,1) + el(1)*acor(i) c----------------------------------------------------------------------- c test for convergence. if m.gt.0, an estimate of the convergence c rate constant is stored in crate, and this is used in the test. c c we first check for a change of iterates that is the size of c roundoff error. if this occurs, the iteration has converged, and a c new rate estimate is not formed. c in all other cases, force at least two iterations to estimate a c local lipschitz constant estimate for adams methods. c on convergence, form pdest = local maximum lipschitz constant c estimate. pdlast is the most recent nonzero estimate. c----------------------------------------------------------------------- 400 continue if (del .le. 100.0d0*pnorm*uround) go to 450 if (m .eq. 0 .and. meth .eq. 1) go to 405 if (m .eq. 0) go to 402 rm = 1024.0d0 if (del .le. 1024.0d0*delp) rm = del/delp rate = dmax1(rate,rm) crate = dmax1(0.2d0*crate,rm) 402 dcon = del*dmin1(1.0d0,1.5d0*crate)/(tesco(2,nq)*conit) if (dcon .gt. 1.0d0) go to 405 pdest = dmax1(pdest,rate/dabs(h*el(1))) if (pdest .ne. 0.0d0) pdlast = pdest go to 450 405 continue m = m + 1 if (m .eq. maxcor) go to 410 if (m .ge. 2 .and. del .gt. 2.0d0*delp) go to 410 delp = del call f (neq, tn, y, savf) nfe = nfe + 1 go to 270 c----------------------------------------------------------------------- c the corrector iteration failed to converge. c if miter .ne. 0 and the jacobian is out of date, pjac is called for c the next try. otherwise the yh array is retracted to its values c before prediction, and h is reduced, if possible. if h cannot be c reduced or mxncf failures have occurred, exit with kflag = -2. c----------------------------------------------------------------------- 410 if (miter .eq. 0 .or. jcur .eq. 1) go to 430 icf = 1 ipup = miter go to 220 430 icf = 2 ncf = ncf + 1 rmax = 2.0d0 tn = told i1 = nqnyh + 1 do 445 jb = 1,nq i1 = i1 - nyh cdir$ ivdep do 440 i = i1,nqnyh 440 yh1(i) = yh1(i) - yh1(i+nyh) 445 continue if (ierpj .lt. 0 .or. iersl .lt. 0) go to 680 if (dabs(h) .le. hmin*1.00001d0) go to 670 if (ncf .eq. mxncf) go to 670 rh = 0.25d0 ipup = miter iredo = 1 go to 170 c----------------------------------------------------------------------- c the corrector has converged. jcur is set to 0 c to signal that the jacobian involved may need updating later. c the local error test is made and control passes to statement 500 c if it fails. c----------------------------------------------------------------------- 450 jcur = 0 if (m .eq. 0) dsm = del/tesco(2,nq) if (m .gt. 0) dsm = vmnorm (n, acor, ewt)/tesco(2,nq) if (dsm .gt. 1.0d0) go to 500 c----------------------------------------------------------------------- c after a successful step, update the yh array. c decrease icount by 1, and if it is -1, consider switching methods. c if a method switch is made, reset various parameters, c rescale the yh array, and exit. if there is no switch, c consider changing h if ialth = 1. otherwise decrease ialth by 1. c if ialth is then 1 and nq .lt. maxord, then acor is saved for c use in a possible order increase on the next step. c if a change in h is considered, an increase or decrease in order c by one is considered also. a change in h is made only if it is by a c factor of at least 1.1. if not, ialth is set to 3 to prevent c testing for that many steps. c----------------------------------------------------------------------- kflag = 0 iredo = 0 nst = nst + 1 hu = h nqu = nq mused = meth do 460 j = 1,l do 460 i = 1,n 460 yh(i,j) = yh(i,j) + el(j)*acor(i) icount = icount - 1 if (icount .ge. 0) go to 488 if (meth .eq. 2) go to 480 c----------------------------------------------------------------------- c we are currently using an adams method. consider switching to bdf. c if the current order is greater than 5, assume the problem is c not stiff, and skip this section. c if the lipschitz constant and error estimate are not polluted c by roundoff, go to 470 and perform the usual test. c otherwise, switch to the bdf methods if the last step was c restricted to insure stability (irflag = 1), and stay with adams c method if not. when switching to bdf with polluted error estimates, c in the absence of other information, double the step size. c c when the estimates are ok, we make the usual test by computing c the step size we could have (ideally) used on this step, c with the current (adams) method, and also that for the bdf. c if nq .gt. mxords, we consider changing to order mxords on switching. c compare the two step sizes to decide whether to switch. c the step size advantage must be at least ratio = 5 to switch. c----------------------------------------------------------------------- if (nq .gt. 5) go to 488 if (dsm .gt. 100.0d0*pnorm*uround .and. pdest .ne. 0.0d0) 1 go to 470 if (irflag .eq. 0) go to 488 rh2 = 2.0d0 nqm2 = min0(nq,mxords) go to 478 470 continue exsm = 1.0d0/dfloat(l) rh1 = 1.0d0/(1.2d0*dsm**exsm + 0.0000012d0) rh1it = 2.0d0*rh1 pdh = pdlast*dabs(h) if (pdh*rh1 .gt. 0.00001d0) rh1it = sm1(nq)/pdh rh1 = dmin1(rh1,rh1it) if (nq .le. mxords) go to 474 nqm2 = mxords lm2 = mxords + 1 exm2 = 1.0d0/dfloat(lm2) lm2p1 = lm2 + 1 dm2 = vmnorm (n, yh(1,lm2p1), ewt)/cm2(mxords) rh2 = 1.0d0/(1.2d0*dm2**exm2 + 0.0000012d0) go to 476 474 dm2 = dsm*(cm1(nq)/cm2(nq)) rh2 = 1.0d0/(1.2d0*dm2**exsm + 0.0000012d0) nqm2 = nq 476 continue if (rh2 .lt. ratio*rh1) go to 488 c the switch test passed. reset relevant quantities for bdf. ---------- 478 rh = rh2 icount = 20 meth = 2 miter = jtyp pdlast = 0.0d0 nq = nqm2 l = nq + 1 go to 170 c----------------------------------------------------------------------- c we are currently using a bdf method. consider switching to adams. c compute the step size we could have (ideally) used on this step, c with the current (bdf) method, and also that for the adams. c if nq .gt. mxordn, we consider changing to order mxordn on switching. c compare the two step sizes to decide whether to switch. c the step size advantage must be at least 5/ratio = 1 to switch. c if the step size for adams would be so small as to cause c roundoff pollution, we stay with bdf. c----------------------------------------------------------------------- 480 continue exsm = 1.0d0/dfloat(l) if (mxordn .ge. nq) go to 484 nqm1 = mxordn lm1 = mxordn + 1 exm1 = 1.0d0/dfloat(lm1) lm1p1 = lm1 + 1 dm1 = vmnorm (n, yh(1,lm1p1), ewt)/cm1(mxordn) rh1 = 1.0d0/(1.2d0*dm1**exm1 + 0.0000012d0) go to 486 484 dm1 = dsm*(cm2(nq)/cm1(nq)) rh1 = 1.0d0/(1.2d0*dm1**exsm + 0.0000012d0) nqm1 = nq exm1 = exsm 486 rh1it = 2.0d0*rh1 pdh = pdnorm*dabs(h) if (pdh*rh1 .gt. 0.00001d0) rh1it = sm1(nqm1)/pdh rh1 = dmin1(rh1,rh1it) rh2 = 1.0d0/(1.2d0*dsm**exsm + 0.0000012d0) if (rh1*ratio .lt. 5.0d0*rh2) go to 488 alpha = dmax1(0.001d0,rh1) dm1 = (alpha**exm1)*dm1 if (dm1 .le. 1000.0d0*uround*pnorm) go to 488 c the switch test passed. reset relevant quantities for adams. -------- rh = rh1 icount = 20 meth = 1 miter = 0 pdlast = 0.0d0 nq = nqm1 l = nq + 1 go to 170 c c no method switch is being made. do the usual step/order selection. -- 488 continue ialth = ialth - 1 if (ialth .eq. 0) go to 520 if (ialth .gt. 1) go to 700 if (l .eq. lmax) go to 700 do 490 i = 1,n 490 yh(i,lmax) = acor(i) go to 700 c----------------------------------------------------------------------- c the error test failed. kflag keeps track of multiple failures. c restore tn and the yh array to their previous values, and prepare c to try the step again. compute the optimum step size for this or c one lower order. after 2 or more failures, h is forced to decrease c by a factor of 0.2 or less. c----------------------------------------------------------------------- 500 kflag = kflag - 1 tn = told i1 = nqnyh + 1 do 515 jb = 1,nq i1 = i1 - nyh cdir$ ivdep do 510 i = i1,nqnyh 510 yh1(i) = yh1(i) - yh1(i+nyh) 515 continue rmax = 2.0d0 if (dabs(h) .le. hmin*1.00001d0) go to 660 if (kflag .le. -3) go to 640 iredo = 2 rhup = 0.0d0 go to 540 c----------------------------------------------------------------------- c regardless of the success or failure of the step, factors c rhdn, rhsm, and rhup are computed, by which h could be multiplied c at order nq - 1, order nq, or order nq + 1, respectively. c in the case of failure, rhup = 0.0 to avoid an order increase. c the largest of these is determined and the new order chosen c accordingly. if the order is to be increased, we compute one c additional scaled derivative. c----------------------------------------------------------------------- 520 rhup = 0.0d0 if (l .eq. lmax) go to 540 do 530 i = 1,n 530 savf(i) = acor(i) - yh(i,lmax) dup = vmnorm (n, savf, ewt)/tesco(3,nq) exup = 1.0d0/dfloat(l+1) rhup = 1.0d0/(1.4d0*dup**exup + 0.0000014d0) 540 exsm = 1.0d0/dfloat(l) rhsm = 1.0d0/(1.2d0*dsm**exsm + 0.0000012d0) rhdn = 0.0d0 if (nq .eq. 1) go to 550 ddn = vmnorm (n, yh(1,l), ewt)/tesco(1,nq) exdn = 1.0d0/dfloat(nq) rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0) c if meth = 1, limit rh according to the stability region also. -------- 550 if (meth .eq. 2) go to 560 pdh = dmax1(dabs(h)*pdlast,0.000001d0) if (l .lt. lmax) rhup = dmin1(rhup,sm1(l)/pdh) rhsm = dmin1(rhsm,sm1(nq)/pdh) if (nq .gt. 1) rhdn = dmin1(rhdn,sm1(nq-1)/pdh) pdest = 0.0d0 560 if (rhsm .ge. rhup) go to 570 if (rhup .gt. rhdn) go to 590 go to 580 570 if (rhsm .lt. rhdn) go to 580 newq = nq rh = rhsm go to 620 580 newq = nq - 1 rh = rhdn if (kflag .lt. 0 .and. rh .gt. 1.0d0) rh = 1.0d0 go to 620 590 newq = l rh = rhup if (rh .lt. 1.1d0) go to 610 r = el(l)/dfloat(l) do 600 i = 1,n 600 yh(i,newq+1) = acor(i)*r go to 630 610 ialth = 3 go to 700 c if meth = 1 and h is restricted by stability, bypass 10 percent test. 620 if (meth .eq. 2) go to 622 if (rh*pdh*1.00001d0 .ge. sm1(newq)) go to 625 622 if (kflag .eq. 0 .and. rh .lt. 1.1d0) go to 610 625 if (kflag .le. -2) rh = dmin1(rh,0.2d0) c----------------------------------------------------------------------- c if there is a change of order, reset nq, l, and the coefficients. c in any case h is reset according to rh and the yh array is rescaled. c then exit from 690 if the step was ok, or redo the step otherwise. c----------------------------------------------------------------------- if (newq .eq. nq) go to 170 630 nq = newq l = nq + 1 iret = 2 go to 150 c----------------------------------------------------------------------- c control reaches this section if 3 or more failures have occured. c if 10 failures have occurred, exit with kflag = -1. c it is assumed that the derivatives that have accumulated in the c yh array have errors of the wrong order. hence the first c derivative is recomputed, and the order is set to 1. then c h is reduced by a factor of 10, and the step is retried, c until it succeeds or h reaches hmin. c----------------------------------------------------------------------- 640 if (kflag .eq. -10) go to 660 rh = 0.1d0 rh = dmax1(hmin/dabs(h),rh) h = h*rh do 645 i = 1,n 645 y(i) = yh(i,1) call f (neq, tn, y, savf) nfe = nfe + 1 do 650 i = 1,n 650 yh(i,2) = h*savf(i) ipup = miter ialth = 5 if (nq .eq. 1) go to 200 nq = 1 l = 2 iret = 3 go to 150 c----------------------------------------------------------------------- c all returns are made through this section. h is saved in hold c to allow the caller to change h on the next step. c----------------------------------------------------------------------- 660 kflag = -1 go to 720 670 kflag = -2 go to 720 680 kflag = -3 go to 720 690 rmax = 10.0d0 700 r = 1.0d0/tesco(2,nqu) do 710 i = 1,n 710 acor(i) = acor(i)*r 720 hold = h jstart = 1 return c----------------------- end of subroutine stoda ----------------------- end double precision function vmnorm (n, v, w) clll. optimize c----------------------------------------------------------------------- c this function routine computes the weighted max-norm c of the vector of length n contained in the array v, with weights c contained in the array w of length n.. c vmnorm = max(i=1,...,n) abs(v(i))*w(i) c----------------------------------------------------------------------- integer n, i double precision v, w, vm dimension v(n), w(n) vm = 0.0d0 do 10 i = 1,n 10 vm = dmax1(vm,dabs(v(i))*w(i)) vmnorm = vm return c----------------------- end of function vmnorm ------------------------ end subroutine prja (neq, y, yh, nyh, ewt, ftem, savf, wm, iwm, 1 f, jac) clll. optimize external f, jac integer neq, nyh, iwm integer iownd, iowns, 1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, 2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu integer iownd2, iowns2, jtyp, mused, mxordn, mxords integer i, i1, i2, ier, ii, j, j1, jj, lenp, 1 mba, mband, meb1, meband, ml, ml3, mu, np1 double precision y, yh, ewt, ftem, savf, wm double precision rowns, 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround double precision rownd2, rowns2, pdnorm double precision con, fac, hl0, r, r0, srur, yi, yj, yjj, 1 vmnorm, fnorm, bnorm dimension neq(1), y(1), yh(nyh,1), ewt(1), ftem(1), savf(1), 1 wm(1), iwm(1) common /ls0001/ rowns(209), 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, 3 iownd(14), iowns(6), 4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu common /lsa001/ rownd2, rowns2(20), pdnorm, 1 iownd2(3), iowns2(2), jtyp, mused, mxordn, mxords c----------------------------------------------------------------------- c prja is called by stoda to compute and process the matrix c p = i - h*el(1)*j , where j is an approximation to the jacobian. c here j is computed by the user-supplied routine jac if c miter = 1 or 4 or by finite differencing if miter = 2 or 5. c j, scaled by -h*el(1), is stored in wm. then the norm of j (the c matrix norm consistent with the weighted max-norm on vectors given c by vmnorm) is computed, and j is overwritten by p. p is then c subjected to lu decomposition in preparation for later solution c of linear systems with p as coefficient matrix. this is done c by dgefa if miter = 1 or 2, and by dgbfa if miter = 4 or 5. c c in addition to variables described previously, communication c with prja uses the following.. c y = array containing predicted values on entry. c ftem = work array of length n (acor in stoda). c savf = array containing f evaluated at predicted y. c wm = real work space for matrices. on output it contains the c lu decomposition of p. c storage of matrix elements starts at wm(3). c wm also contains the following matrix-related data.. c wm(1) = sqrt(uround), used in numerical jacobian increments. c iwm = integer work space containing pivot information, starting at c iwm(21). iwm also contains the band parameters c ml = iwm(1) and mu = iwm(2) if miter is 4 or 5. c el0 = el(1) (input). c pdnorm= norm of jacobian matrix. (output). c ierpj = output error flag, = 0 if no trouble, .gt. 0 if c p matrix found to be singular. c jcur = output flag = 1 to indicate that the jacobian matrix c (or approximation) is now current. c this routine also uses the common variables el0, h, tn, uround, c miter, n, nfe, and nje. c----------------------------------------------------------------------- nje = nje + 1 ierpj = 0 jcur = 1 hl0 = h*el0 go to (100, 200, 300, 400, 500), miter c if miter = 1, call jac and multiply by scalar. ----------------------- 100 lenp = n*n do 110 i = 1,lenp 110 wm(i+2) = 0.0d0 call jac (neq, tn, y, 0, 0, wm(3), n) con = -hl0 do 120 i = 1,lenp 120 wm(i+2) = wm(i+2)*con go to 240 c if miter = 2, make n calls to f to approximate j. -------------------- 200 fac = vmnorm (n, savf, ewt) r0 = 1000.0d0*dabs(h)*uround*dfloat(n)*fac if (r0 .eq. 0.0d0) r0 = 1.0d0 srur = wm(1) j1 = 2 do 230 j = 1,n yj = y(j) r = dmax1(srur*dabs(yj),r0/ewt(j)) y(j) = y(j) + r fac = -hl0/r call f (neq, tn, y, ftem) do 220 i = 1,n 220 wm(i+j1) = (ftem(i) - savf(i))*fac y(j) = yj j1 = j1 + n 230 continue nfe = nfe + n 240 continue c compute norm of jacobian. -------------------------------------------- pdnorm = fnorm (n, wm(3), ewt)/dabs(hl0) c add identity matrix. ------------------------------------------------- np1 = n + 1 j = 3 do 250 i = 1,n wm(j) = wm(j) + 1.0d0 250 j = j + np1 c do lu decomposition on p. -------------------------------------------- call dgefa (wm(3), n, n, iwm(21), ier) if (ier .ne. 0) ierpj = 1 return c dummy block only, since miter is never 3 in this routine. ------------ 300 return c if miter = 4, call jac and multiply by scalar. ----------------------- 400 ml = iwm(1) mu = iwm(2) ml3 = ml + 3 mband = ml + mu + 1 meband = mband + ml lenp = meband*n do 410 i = 1,lenp 410 wm(i+2) = 0.0d0 call jac (neq, tn, y, ml, mu, wm(ml3), meband) con = -hl0 do 420 i = 1,lenp 420 wm(i+2) = wm(i+2)*con go to 570 c if miter = 5, make mband calls to f to approximate j. ---------------- 500 ml = iwm(1) mu = iwm(2) mband = ml + mu + 1 mba = min0(mband,n) meband = mband + ml meb1 = meband - 1 srur = wm(1) fac = vmnorm (n, savf, ewt) r0 = 1000.0d0*dabs(h)*uround*dfloat(n)*fac if (r0 .eq. 0.0d0) r0 = 1.0d0 do 560 j = 1,mba do 530 i = j,n,mband yi = y(i) r = dmax1(srur*dabs(yi),r0/ewt(i)) 530 y(i) = y(i) + r call f (neq, tn, y, ftem) do 550 jj = j,n,mband y(jj) = yh(jj,1) yjj = y(jj) r = dmax1(srur*dabs(yjj),r0/ewt(jj)) fac = -hl0/r i1 = max0(jj-mu,1) i2 = min0(jj+ml,n) ii = jj*meb1 - ml + 2 do 540 i = i1,i2 540 wm(ii+i) = (ftem(i) - savf(i))*fac 550 continue 560 continue nfe = nfe + mba 570 continue c compute norm of jacobian. -------------------------------------------- pdnorm = bnorm (n, wm(3), meband, ml, mu, ewt)/dabs(hl0) c add identity matrix. ------------------------------------------------- ii = mband + 2 do 580 i = 1,n wm(ii) = wm(ii) + 1.0d0 580 ii = ii + meband c do lu decomposition of p. -------------------------------------------- call dgbfa (wm(3), meband, n, ml, mu, iwm(21), ier) if (ier .ne. 0) ierpj = 1 return c----------------------- end of subroutine prja ------------------------ end double precision function bnorm (n, a, nra, ml, mu, w) clll. optimize c----------------------------------------------------------------------- c this function computes the norm of a banded n by n matrix, c stored in the array a, that is consistent with the weighted max-norm c on vectors, with weights stored in the array w. c ml and mu are the lower and upper half-bandwidths of the matrix. c nra is the first dimension of the a array, nra .ge. ml+mu+1. c in terms of the matrix elements a(i,j), the norm is given by.. c bnorm = max(i=1,...,n) ( w(i) * sum(j=1,...,n) abs(a(i,j))/w(j) ) c----------------------------------------------------------------------- integer n, nra, ml, mu integer i, i1, jlo, jhi, j double precision a, w double precision an, sum dimension a(nra,n), w(n) an = 0.0d0 do 20 i = 1,n sum = 0.0d0 i1 = i + mu + 1 jlo = max0(i-ml,1) jhi = min0(i+mu,n) do 10 j = jlo,jhi 10 sum = sum + dabs(a(i1-j,j))/w(j) an = dmax1(an,sum*w(i)) 20 continue bnorm = an return c----------------------- end of function bnorm ------------------------- end double precision function fnorm (n, a, w) clll. optimize c----------------------------------------------------------------------- c this function computes the norm of a full n by n matrix, c stored in the array a, that is consistent with the weighted max-norm c on vectors, with weights stored in the array w.. c fnorm = max(i=1,...,n) ( w(i) * sum(j=1,...,n) abs(a(i,j))/w(j) ) c----------------------------------------------------------------------- integer n, i, j double precision a, w, an, sum dimension a(n,n), w(n) an = 0.0d0 do 20 i = 1,n sum = 0.0d0 do 10 j = 1,n 10 sum = sum + dabs(a(i,j))/w(j) an = dmax1(an,sum*w(i)) 20 continue fnorm = an return c----------------------- end of function fnorm ------------------------- end