subroutine newton_rc ( fx, ido, n, x ) c*********************************************************************72 c cc NEWTON_RC solves a small system of nonlinear equations. c c Discussion: c c NEWTON_RC uses Newton's method, with the jacobian approximated via c finite differences. c c NEWTON_RC uses "reverse communication". That is, NEWTON_RC does not call c any user subroutines. Instead, the user repeatedly calls NEWTON_RC, c and the user and NEWTON_RC communicate via values of the parameter IDO. c c To begin using NEWTON_RC, set IDO=0, set X to an approximate root, c set FX(1) to a small positive tolerance, and call NEWTON_RC. c c NEWTON_RC will return with IDO=2. Evaluate the residual at the c point X that is returned by NEWTON_RC, store that value in FX, c and call NEWTON_RC back with IDO=2. c c NEWTON_RC will return with IDO=1. Again, evaluate the residual c at the point X returned by NEWTON_RC, store that value in FX< c and call NEWTON_RC back with IDO=1. c c This process will be repeated until NEWTON_RC returns with c IDO=0, meaning that the value X is a good approximation to c the root, or IDO=-1, in which case the algorithm failed. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 07 January 2009 c c Author: c c John Burkardt c c Parameters: c c Input, double precision FX(N), should be set to the c residual of the nonlinear equations if NEWTON_RC returned c with IDO=1 or 2. c c Input/output, integer IDO. c On input, the values of IDO mean: c * 0, this is a new problem. c N contains the number of equations. c X contains a starting guess for the solution. c FX(1) contains a tolerance for the residuals. c * 1 or 2, the user has evaluated the residual c as requested by NEWTON_RC. c N and X are unchanged. c FX contains the value of the residual of X. c On output from NEWTON_RC, the values of IDO mean: c * 0, the output value X is a good approximation to the solution. c * 1, please evaluate the residual of the output value X, c store that value in FX, and call NEWTON_RC with IDO=1. c The jacobian is being approximated. c * 2, please evaluate the residual of the output value X, c store that value in FX, and call NEWTON_RC with IDO=2. c * -1, the algorithm failed. c c Input, integer N, the number of equations. c c Input/output, double precision X(N). c implicit none integer maxn parameter ( maxn = 10 ) integer n double precision delx(maxn) double precision delxj double precision eps double precision fprime(maxn,maxn) double precision fx(n) double precision fxold(maxn) double precision fxnrm double precision fxtol integer i integer idamax integer ido integer info integer ipivot(maxn) integer j integer job integer ncall double precision r8_epsilon double precision x(n) double precision xold save delxj save eps save fprime save fxold save fxtol save j save ncall save xold ncall = ncall + 1 if ( 500 .lt. ncall ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'NEWTON_RC - Fatal error!' write ( *, '(a,i6,a)' ) ' NEWTON_RC called ', ncall, ' times.' stop end if c c IDO=0, user is starting a new problem. c if ( ido .eq. 0 ) then if ( maxn .lt. n ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'NEWTON_RC - Fatal error!' write ( *, '(a)' ) ' N is too large!' write ( *, '(a,i6)' ) ' Input N = ', n write ( *, '(a,i6)' ) ' Maximum legal N = ', maxn ido=-1 return end if eps = sqrt ( r8_epsilon ( ) ) if ( fx(1) .le. 0.0 ) then fxtol = eps else fxtol = fx(1) end if ncall = 0 ido = 2 c c IDO=1, user is returning F(X+delX(J)) for a jacobian estimation c else if ( ido .eq. 1 ) then do i = 1, n fprime(i,j) = ( fx(i) - fxold(i) ) / delxj end do x(j) = xold if ( j .lt. n ) then j = j + 1 delxj = eps*(abs(x(j))+1.0) xold = x(j) x(j) = x(j)+delxj else call dgefa(fprime,maxn,n,ipivot,info) if(info.ne.0)then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'NEWTON_RC - Fatal error!' write ( *, '(a)' ) ' Estimated jacobian is singular.' ido=-1 return end if job=0 do i=1,n delx(i)=-fx(i) enddo call dgesl(fprime,maxn,n,ipivot,delx,job) do i=1,n x(i)=x(i)+delx(i) enddo ido=2 end if c c IDO=2, user is returning F(X) for a convergence evaluation. c else if (ido.eq.2)then i = idamax ( n, fx, 1 ) fxnrm = abs ( fx(i) ) write ( *, '(a,i3,a,g14.6)' ) '|FX(',i,')|=', fxnrm if(fxnrm.le.fxtol)then ido=0 else ido=1 do i=1,n fxold(i)=fx(i) end do j=1 delxj=eps*(abs(x(j))+1.0) xold=x(j) x(j)=x(j)+delxj end if c c Unexpected value of IDO. c else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'NEWTON_RC - Fatal error.' write ( *, '(a,i6)' ) ' Unexpected value of IDO = ', ido stop end if return end function r8_epsilon ( ) c*********************************************************************72 c cc R8_EPSILON returns the R8 roundoff unit. c c Discussion: c c The roundoff unit is a number R which is a power of 2 with the c property that, to the precision of the computer's arithmetic, c 1 .lt. 1 + R c but c 1 = ( 1 + R / 2 ) c c FORTRAN90 provides the superior library routine c c EPSILON ( X ) c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 01 September 2012 c c Author: c c John Burkardt c c Parameters: c c Output, double precision R8_EPSILON, the R8 roundoff unit. c implicit none double precision r8_epsilon r8_epsilon = 2.220446049250313D-016 return end subroutine daxpy ( n, da, dx, incx, dy, incy ) c*********************************************************************72 c cc DAXPY computes constant times a vector plus a vector. c c Discussion: c c This routine uses double precision real arithmetic. c c This routine uses unrolled loops for increments equal to one. c c Modified: c c 08 July 2007 c c Author: c c Jack Dongarra c c Reference: c c Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, c LINPACK User's Guide, c SIAM, 1979, c ISBN13: 978-0-898711-72-1, c LC: QA214.L56. c c Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, c Basic Linear Algebra Subprograms for FORTRAN usage, c ACM Transactions on Mathematical Software, c Volume 5, Number 3, pages 308-323, 1979. c c Parameters: c c Input, integer N, the number of elements in DX and DY. c c Input, double precision DA, the multiplier of DX. c c Input, double precision DX(*), the first vector. c c Input, integer INCX, the increment between successive entries of DX. c c Input/output, double precision DY(*), the second vector. c On output, DY(*) has been replaced by DY(*) + DA * DX(*). c c Input, integer INCY, the increment between successive entries of DY. c implicit none double precision da double precision dx(*) double precision dy(*) integer i integer incx integer incy integer ix integer iy integer m integer n if ( n .le. 0 ) then return end if if ( da .eq. 0.0d0 ) then return end if 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 i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy end do 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 i = 1,m dy(i) = dy(i) + da*dx(i) end do if( n .lt. 4 ) return 40 continue do i = m + 1, 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) end do return end function ddot ( n, dx, incx, dy, incy ) c*********************************************************************72 c cc DDOT forms the dot product of two vectors. c c Discussion: c c This routine uses double precision real arithmetic. c c This routine uses unrolled loops for increments equal to one. c c Modified: c c 07 July 2007 c c Author: c c Jack Dongarra c c Reference: c c Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, c LINPACK User's Guide, c SIAM, 1979, c ISBN13: 978-0-898711-72-1, c LC: QA214.L56. c c Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, c Basic Linear Algebra Subprograms for FORTRAN usage, c ACM Transactions on Mathematical Software, c Volume 5, Number 3, pages 308-323, 1979. c c Parameters: c c Input, integer N, the number of entries in the vectors. c c Input, double precision DX(*), the first vector. c c Input, integer INCX, the increment between successive entries in DX. c c Input, double precision DY(*), the second vector. c c Input, integer INCY, the increment between successive entries in DY. c c Output, double precision DDOT, the sum of the product of the c corresponding entries of DX and DY. c implicit none double precision ddot double precision dx(*),dy(*),dtemp integer i,incx,incy,ix,iy,m,n 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 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 i = 1,n dtemp = dtemp + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy end do 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 i = 1,m dtemp = dtemp + dx(i)*dy(i) end do if( n .lt. 5 ) go to 60 40 continue do i = m+1, 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) end do 60 ddot = dtemp return end subroutine dgefa(a,lda,n,ipvt,info) c*********************************************************************72 c cc 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 integer lda,n,ipvt(1),info double precision a(lda,1) double precision t integer idamax,j,k,kp1,l,nm1 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 dgesl(a,lda,n,ipvt,b,job) c*********************************************************************72 c cc DGESL solves a linear system factored by DGEFA. 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 integer lda,n,ipvt(1),job double precision a(lda,1),b(1) double precision ddot,t integer k,kb,l,nm1 c nm1 = n - 1 if (job .ne. 0) go to 50 c c job = 0 , solve a * x = b c first solve l*y = b c if (nm1 .lt. 1) go to 30 do 20 k = 1, nm1 l = ipvt(k) t = b(l) if (l .eq. k) go to 10 b(l) = b(k) b(k) = t 10 continue call daxpy(n-k,t,a(k+1,k),1,b(k+1),1) 20 continue 30 continue c c now solve u*x = y c do 40 kb = 1, n k = n + 1 - kb b(k) = b(k)/a(k,k) t = -b(k) call daxpy(k-1,t,a(1,k),1,b(1),1) 40 continue go to 100 50 continue c c job = nonzero, solve trans(a) * x = b c first solve trans(u)*y = b c do 60 k = 1, n t = ddot(k-1,a(1,k),1,b(1),1) b(k) = (b(k) - t)/a(k,k) 60 continue c c now solve trans(l)*x = y c if (nm1 .lt. 1) go to 90 do 80 kb = 1, nm1 k = n - kb b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1) l = ipvt(k) if (l .eq. k) go to 70 t = b(l) b(l) = b(k) b(k) = t 70 continue 80 continue 90 continue 100 continue return end subroutine dscal ( n, da, dx, incx ) c*********************************************************************72 c cc DSCAL scales a vector by a constant. c c Discussion: c c This routine uses double precision real arithmetic. c c Modified: c c 07 July 2007 c c Author: c c Jack Dongarra c c Reference: c c Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, c LINPACK User's Guide, c SIAM, 1979, c ISBN13: 978-0-898711-72-1, c LC: QA214.L56. c c Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, c Basic Linear Algebra Subprograms for FORTRAN usage, c ACM Transactions on Mathematical Software, c Volume 5, Number 3, pages 308-323, 1979. c c Parameters: c c Input, integer N, the number of entries in the vector. c c Input, double precision SA, the multiplier. c c Input/output, double precision X(*), the vector to be scaled. c c Input, integer INCX, the increment between successive entries of X. c implicit none double precision da,dx(*) integer i,incx,m,n,nincx if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do i = 1,nincx,incx dx(i) = da*dx(i) end do 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 i = 1,m dx(i) = da*dx(i) end do if( n .lt. 5 ) return 40 continue do i = m+1, 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) end do return end subroutine dswap ( n, dx, incx, dy, incy ) c*********************************************************************72 c cc DSWAP interchanges two vectors. c c Discussion: c c This routine uses double precision real arithmetic. c c Modified: c c 07 July 2007 c c Author: c c Jack Dongarra c c Reference: c c Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, c LINPACK User's Guide, c SIAM, 1979, c ISBN13: 978-0-898711-72-1, c LC: QA214.L56. c c Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, c Basic Linear Algebra Subprograms for FORTRAN usage, c ACM Transactions on Mathematical Software, c Volume 5, Number 3, pages 308-323, 1979. c c Parameters: c c Input, integer N, the number of entries in the vectors. c c Input/output, double precision X(*), one of the vectors to swap. c c Input, integer INCX, the increment between successive entries of X. c c Input/output, double precision Y(*), one of the vectors to swap. c c Input, integer INCY, the increment between successive elements of Y. c implicit none double precision dx(*),dy(*),dtemp integer i,incx,incy,ix,iy,m,n 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 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 i = 1,n dtemp = dx(ix) dx(ix) = dy(iy) dy(iy) = dtemp ix = ix + incx iy = iy + incy end do return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,3) if( m .eq. 0 ) go to 40 do i = 1,m dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp end do if( n .lt. 3 ) return 40 continue do i = m+1, n, 3 dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp dtemp = dx(i + 1) dx(i + 1) = dy(i + 1) dy(i + 1) = dtemp dtemp = dx(i + 2) dx(i + 2) = dy(i + 2) dy(i + 2) = dtemp end do return end function idamax ( n, dx, incx ) c*********************************************************************72 c cc IDAMAX finds the index of element having maximum absolute value. c c Discussion: c c This routine uses double precision real arithmetic. c c Modified: c c 07 July 2007 c c Author: c c Jack Dongarra c c Reference: c c Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, c LINPACK User's Guide, c SIAM, 1979, c ISBN13: 978-0-898711-72-1, c LC: QA214.L56. c c Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, c Basic Linear Algebra Subprograms for FORTRAN usage, c ACM Transactions on Mathematical Software, c Volume 5, Number 3, pages 308-323, 1979. c c Parameters: c c Input, integer N, the number of entries in the vector. c c Input, double precision X(*), the vector to be examined. c c Input, integer INCX, the increment between successive entries of SX. c c Output, integer IDAMAX, the index of the element of SX of maximum c absolute value. c implicit none double precision dx(*),dmax integer idamax integer i,incx,ix,n idamax = 0 if( n.lt.1 .or. incx.le.0 ) return idamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 dmax = dabs(dx(1)) ix = ix + incx do i = 2,n if(dabs(dx(ix)).le.dmax) go to 5 idamax = i dmax = dabs(dx(ix)) 5 ix = ix + incx end do return c c code for increment equal to 1 c 20 dmax = dabs(dx(1)) do i = 2,n if( dmax .lt. dabs(dx(i)) ) then idamax = i dmax = dabs(dx(i)) end if end do return end subroutine timestamp ( ) c*********************************************************************72 c cc TIMESTAMP prints out the current YMDHMS date as a timestamp. c c Discussion: c c This FORTRAN77 version is made available for cases where the c FORTRAN90 version cannot be used. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 12 January 2007 c c Author: c c John Burkardt c c Parameters: c c None c implicit none character * ( 8 ) ampm integer d character * ( 8 ) date integer h integer m integer mm character * ( 9 ) month(12) integer n integer s character * ( 10 ) time integer y save month data month / & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' / call date_and_time ( date, time ) read ( date, '(i4,i2,i2)' ) y, m, d read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm if ( h .lt. 12 ) then ampm = 'AM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h .lt. 12 ) then ampm = 'PM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, & '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, month(m), y, h, ':', n, ':', s, '.', mm, ampm return end