********************************************************** Line 1 C C October 24 2011 C ========================================== C C Main Program DESCON C ========================== C C MINPACK2 APPLICATIONS C Minimal Surface Area Problem C ========================================== C C DESCON is a subroutine dedicated to compute the minimizer of C a differentiable function with a large number of variables. C It is supposed that we have the algebraic expression of the C function and its gradient. C C DESCON implements an accelerated conjugate gradient algorithm C that at each iteration both the descent and the conjugacy C conditions are guaranteed. C C The algorithm introduces the modified Wolfe line search in which C the parameter in the second Wolfe condition is modified at C each iteration. C C The algorithm is described in: C N. Andrei, Another conjugate gradient algorithm with C guaranteed descent and conjugacy conditions for large-scale C unconstrained optimization. C ICI Technical Report, November 9, 2011. C C---------------------------------------------------------------- C Neculai Andrei C Research Institute for Informatics C Center for Advanced Modeling and Optimization C and C Academy of Romanian Scientists ***************************************************************** * integer nexp,n,nx,ny,iter,irs, fgcnt,lscnt integer maxiter, maxfg integer itert, irstot, fgtot,lstot, stoptest integer*4 gh,gm,gs,gc, ght,gmt,gst,gct, timpexp double precision epsg, dc,cc real*8 f, gnorm, timp, proc real*8 epsm integer*2 iyear, imonth, iday character*40 numef logical angle, powell common /acca/epsm C LOCAL ARRAYS double precision x(5000000), grad(5000000) real*8 bottom(2000), top(2000), left(2000), right(2000) * Output files: open(unit=4,file='a7.out',status='unknown') open(unit=5,file='a7.rez',status='unknown') open(unit=7,file='a7.sol',status='unknown') * Epsilon machine computation (Please see the common /acca/) epsm = 0.111e-26 c-------------------------------------------------------------------- c Initial values of parameters: c dc - descent condition (w) c cc - conjucagy condition parameter (v) dc = 0.875d0 cc = 0.05d0 c------------------------------------------------------------------- c Restart parameters: angle or Powell angle = .false. powell = .true. nexp = 1 numef = 'Minimal Surface Area' c c Discretization steps. Parameter of the problem. c do nx = 1000,1000 ny=nx n = nx * ny c------------------------------------------------------------------- c stoptest = option parameter for selection the * stopping criterion: * 1: if(ginf .lt. epsg) * 2: if(gnorm .le. epsg) * where: * ginf = infinite norm of gradient g(xk), * gnorm = norm 2 of gradient g(xk). * stoptest = 1 c------------------------------------------------------------------------- c c Some convergence constant epsg = (10.d0)**(-6) write(4,18) 18 format(4x,'***************************************************',/, * 4x,'* DESCON Neculai Andrei ***',/, * 4x,'* Another Accelerated Conjugate Gradient with ***',/, * 4x,'* Guaranteed Descent and Conjugacy conditions. ***',/, * 4x,'* --------------------------------------------- ***',/, * 4x,'* Project: FCGA ***',/, * 4x,'* The Fastest Conjugate Gradient Algorithm ***',/, * 4x,'***************************************************',/) write(4,19) dc, cc 19 format(4x,'Parameters: Descent=',e20.13,4x,'Conjugacy=',e20.13,/) call getdat(iyear, imonth, iday) write(4,22) imonth, iday, iyear 22 format(4x,'Date: --- Month:',i2,' Day:',i2,' Year: ',i4,/) if(powell) then write(4,20) 20 format(4x,'DESCON. Powell restart.') end if if(angle) then write(4,21) 21 format(4x,'DESCON. Angle restart.') end if write(4,26) * *------------------------* *--------------------------------------------- | Here Start Experiments | * *------------------------* write(4,941) nexp, numef, stoptest write(*,941) nexp, numef, stoptest 941 format(/,5x,i2,4x,'DLDC Algorithm. ','Function:',a40,/, * 11x,'stoptest=',i2,/) * *---------------------------------------------------------------- * write(4,25) 25 format(7x,'n',3x,'iter',3x,'irs',2x,'fgcnt',1x,'lscnt', * 3x,'time(c)',8x,'fxnew',14x,'gnorm',9x,'ib',3x,'it') write(4,26) 26 format(1x,91('-')) * *---------------------------------------------------------------- * itert=0 irstot=0 fgtot=0 lstot=0 timp=0.d0 * C Maximum number of iterations and function evaluations maxiter = 10000 maxfg = 15000 call dmsabc(nx,ny,bottom,top,left,right) * Initial point call dmsafg(nx,ny,x,f,grad,'XS',bottom,top,left,right) * Call the solver call gettim(gh,gm,gs,gc) call descon(n,x,epsg,maxiter,maxfg,f,gnorm,stoptest,dc,cc, * iter,irs,fgcnt,lscnt, angle, powell, * ibeta, itheta, nexp, nx,ny) call gettim(ght,gmt,gst,gct) call exetim(gh,gm,gs,gc, ght,gmt,gst,gct) *------------------------------------ Compute elapsed time in centeseconds timpexp = ght*360000 + gmt*6000 + gst*100 + gct C Write statistics C *---------------------------------------------------------- Write *.out write(4,99) n,iter,irs,fgcnt,lscnt,timpexp,f,gnorm,ibeta,itheta 99 format(1x,i7,1x,i6,1x,i5,1x,i6,1x,i5,1x, * i9,e20.13,e20.13,i3,2x,i3) *---------------------------------------------------------- Write *.rez if(n .eq. 1000) then write(5,101) nexp,n,iter,fgcnt,timpexp,f 101 format(i2,i6,2x,i6,2x,i6,2x,i6,2x,e20.13) else write(5,102) n,iter,fgcnt,timpexp,f 102 format(2x,i6,2x,i6,2x,i6,2x,i6,2x,e20.13) end if *------------------------------------------------------------------ itert = itert + iter irstot= irstot + irs fgtot = fgtot + fgcnt lstot = lstot + lscnt timp = timp + float(timpexp) write(*,110) n, iter, irs, fgcnt, lscnt,timpexp,f,gnorm 110 format(1x,i7,1x,i6,1x,i5,1x,i6,1x,i5,1x,i9,' c',e13.6,e13.6) proc = float(irstot)*100.d0/float(itert) write(4,26) write(4,150) itert, irstot, fgtot, lstot, timp/100.d0, proc write(*,155) itert, irstot, fgtot, lstot, timp/100.d0, proc 150 format(3x,'TOTAL',1x,i6,1x,i5,1x,i6,1x,i5,3x,f7.2,' (seconds)', * 4x,'proc= ',f6.2,'%',/) 155 format(3x,'TOTAL',1x,i6,1x,i5,1x,i6,1x,i5,3x,f7.2,' (seconds)', * 4x,'proc= ',f6.2,'%',/) c---------------------------------------------------------- End do nx end do c write(5,160) 160 format(1x,55('-')) write(5,161) 161 format(2x,' n',2x,' iter',' fgcnt time(c)',10x,'fx',/) write(5,19) dc, cc write(7,175)(x(i),i=1,n) 175 format(1000(1x,f10.8)) stop end c********************************************************* END main program c *** Conjugate Gradient Algorithms Project *** C =============================================== c c * * *------------------------------------------------------------------- * Subroutine DESCON * ===================== * * * Neculai Andrei * Research Institute for Informatics * Center for Advanced Modeling and Optimization * 8-10, Averescu Avenue, Bucharest 1, Romania * E-mail: nandrei@ici.ro * voice: 666.58.70 * and * Academy of Romanian Scientists * Science and Information Technology Section * 54, Splaiul Independentei, Bucharest 5, Romania * * * * /-----------------------------------------------------------------\ * | DESCON is a subroutine dedicated to compute the minimizer of | * | a differentiable function with a large number of variables. | * | | * | This subroutine is accompanied by subroutine "LineSearch" which | * | implements the Wolfe line search. Both these subroutines belong | * | to DESCON package. | * | | * | The user must provide a subroutine to evaluate the function | * | value and its gradient in a point. The name of this subroutine | * | is EVALFG. | * | The algebraic expression of the functions considered in this | * | program, as well as their Fortran code are presented into the | * | paper: | * | N. Andrei, "An unconstrained optimization test functions | * | collection", Advanced Modeling and Optimization, vol.10, No.1, | * | (2008) pp.147-161. | * | http://www.ici.ro/camo/journal/v10n1.htm | * | | * | There are some facilities for the user to specify: | * | 1) The termination criterion. | * | 2) Convergence tolerance for gradient. | * | 3) Convergence tolerance for function value. | * | 4) Maximum number of iterations in LineSearch subroutine. | * | 5) The parameter w from the descent condition. | * | 6) The parameter v from the conjugacy condition. | * | | * |-----------------------------------------------------------------| * | | * | The calling sequence of DESCON is: | * | | * | subroutine descon(n,x,epsg,maxiter,maxfg,f,gnorm,stoptest,dc,| * | cc,iter,irs,fgcnt,lscnt,angle, powell, | * | ibeta, itheta, nexp) | * | | * |Input parameters: | * |================= | * |n (integer) number of variables. | * |x (double) starting guess, length n. On output | * | contains the solution. | * |epsg (double) convergence tolerance for gradient. | * |maxiter (integer) maximum number of iterations. | * |maxfg (integer) maxumum number of function and its gradient | * | evaluations. | * |stoptest = option parameter for selection of | * | stopping criterion: | * | if stoptest = 1 then consider the following test: | * | if(ginf .le. epsg) | * | if stoptest = 2 then consider the following test: | * | if(gnorm .le. epsg) | * | where: | * | ginf = infinite norm of gradient g(xk), | * | gnorm = norm-2 of gradient g(xk). | * |dc (double) parameter for sufficient descent condition | * |cc (double) parameter for conjugacy condition | * |angle (logical) parameter specifying the angle criterion of | * | restart. | * |powell (logical) parameter specifying the Powell criterion of| * | restart. | * |nexp (integer) parameter specifying the number of the | * | problem considered in a train of experiments| * | | * | | * |Output parameters: | * |================== | * |f (double) function value in final (optimal) point. | * |gnorm (double) norm-2 of gradient in final point. | * |iter (integer) number of iterations to get the final point.| * |irs (integer) number of restart iterations. | * |fgcnt (integer) number of function evaluations. | * |lscnt (integer) number of line searches. | * |ibeta (integer) number of iterations in which beta is zero | * |itheta (integer) number of iterations in which theta is one | * |-----------------------------------------------------------------| * | | * | | * |Calling subroutines: | * |==================== | * |Subroutine DESCON is calling two subroutines: | * |EVALFG an user subroutine (function and gradient), | * |LINESEARCH a package subroutine. | * | | * |The user must supply a subroutine with the function and its | * |gradient: | * | call evalfg(n,x,fx,grad,nexp) | * |where: | * | n (integer) number of variables. | * | x (double) the current iteration. | * | fx (double) function value in point x. | * | grad (double) array with gradient of function in point x. | * | nexp (integer) parameter specifying the number of the | * | problem considered in a train of experiments. | * | | * | | * |The Wolfe line search is implemented in the subroutine | * |LINESEARCH, which belongs to the package. | * |The calling sequence of this subroutine is as follows: | * | call LineSearch (n,x,f,d,gtd,dnorm,alpha,xnew,fnew,gnew,sigma, | * | fgcnt,lscnt,lsflag, nexp) | * |where: | * | n (integer) number of variables. | * | x (double) the current iteration. | * | f (double) function value in current point. | * | d (double) array with search direction. | * | gtd (double) scalar: grad'*d. | * | dnorm (double) 2 norm of d. | * | alpha (double) step length (given by the LineSearch). | * | xnew (double) array with the new estimation of variables. | * | fnew (double) function value in xnew. | * | gnew (double) array with gradient in xnew. | * | sigma (double) parameter sigma in the second Wolfe line | * | search condition. (input parameter) | * | fgcnt (integer) number of function evaluations. | * | lscnt (integer) number of line searches. | * | lsflag (integer) parameter for abnormal Line Search | * | Termination. If the # of iterations in | * | LineSearch is greater than 20 then lsflag=1 | * | nexp (integer) parameter specifying the number of the | * | problem considered in a train of experiments | * | | * |Subroutine LINESEARCH is the same as that used in CONMIN package | * |by Shanno and Phua, and in SCALCG package by Andrei. | * | | * |-----------------------------------------------------------------| * | Neculai Andrei, 2011 | * \-----------------------------------------------------------------/ * * c c******************************************************************** subroutine descon(n,x,epsg,maxiter,maxfg,f,gnorm,stoptest,dc,cc, * iter,irs,fgcnt,lscnt,angle, powell, * ibeta, itheta, nexp, nx,ny) parameter(ia=1000000) C SCALAR ARGUMENTS integer n,iter,irs,fgcnt,lscnt,maxiter,maxfg integer stoptest, nexp double precision epsg, f,gnorm logical angle, powell C ARRAY ARGUMENTS double precision x(n) C LOCAL SCALARS integer i,lsflag double precision fnew,alpha,gtg, beta, + dnorm,dnormprev, ginf, + gtd, gtgp, stg, acc,bdc, + yts, ytg, theta, delta, dc,cc, epsm, bn,sigma C LOCAL ARRAYS double precision xnew(ia),g(ia),gnew(ia),d(ia), + y(ia), s(ia) real*8 bottom(2000), top(2000), left(2000), right(2000) common /acca/epsm * Initialization n5 = mod(n,5) n6 = n5 + 1 iter = 0 irs = 0 fgcnt = 0 lscnt = 0 ibeta = 0 itheta = 0 c call evalfg(n,x,f,g, nexp) call dmsafg(nx,ny,x,f,g,'FG',bottom,top,left,right) fgcnt = fgcnt + 1 gtg = 0.0d0 do i = 1,n5 d(i) = - g(i) gtg = gtg + g(i) ** 2 end do do i = n6,n,5 d(i) = -g(i) d(i+1) = -g(i+1) d(i+2) = -g(i+2) d(i+3) = -g(i+3) d(i+4) = -g(i+4) gtg = gtg + g(i)**2 + g(i+1)**2 + g(i+2)**2 + * g(i+3)**2 + g(i+4)**2 end do gnorm = sqrt( gtg ) gtd = -gtg dnorm = gnorm if ( gnorm .gt. 0.0d0 ) then alpha = 1.0d0 / dnorm end if sigma = 0.8d0 C ------------------------------ Main loop -------------------- C==================================================================== 110 continue c------------------------------------------------ STOP TEST if(iter .eq. 0) go to 91 if(stoptest .eq. 1) then ginf=dabs(g(1)) do i=2,n5 if(dabs(g(i)) .gt. ginf) ginf = dabs(g(i)) end do do i=n6,n,5 if(dabs(g(i)) .gt. ginf) ginf = dabs(g(i)) if(dabs(g(i+1)) .gt. ginf) ginf = dabs(g(i+1)) if(dabs(g(i+2)) .gt. ginf) ginf = dabs(g(i+2)) if(dabs(g(i+3)) .gt. ginf) ginf = dabs(g(i+3)) if(dabs(g(i+4)) .gt. ginf) ginf = dabs(g(i+4)) end do if(ginf .le. epsg) go to 999 end if c if(stoptest .eq. 2) then if(gnorm .le. epsg) go to 999 end if c 91 continue c------------------------------------------------ END STOP TEST * Increment iteration iter = iter + 1 if(iter .gt. maxiter) go to 999 * Line search * Determine the step length ALPHA and the new point XNEW, as well as * the function value in xnew, FNEW, and the gradient in xnew, GNEW. call LineSearch(n,x,f,d,gtd,dnorm,alpha,xnew,fnew,gnew, + sigma,fgcnt,lscnt,lsflag, n5,n6,nexp,nx,ny) if(fgcnt .gt. maxfg) go to 999 alpha = alpha * *---------------------------------- Acceleration section C ==================== c c (We use an/bn only if b is different from zero) c (an=gtd) bn=0.d0 do i = 1,n5 bn = bn + (g(i)-gnew(i))*d(i) end do do i = n6,n,5 bn= bn + (g(i)-gnew(i))*d(i)+ * (g(i+1)-gnew(i+1))*d(i+1)+ * (g(i+2)-gnew(i+2))*d(i+2)+ * (g(i+3)-gnew(i+3))*d(i+3)+ * (g(i+4)-gnew(i+4))*d(i+4) end do * if(dabs(bn) .gt. epsm) then do i=1,n5 xnew(i) = x(i) + (gtd/bn)*alpha*d(i) end do do i=n6,n,5 xnew(i) = x(i) + (gtd/bn)*alpha*d(i) xnew(i+1) = x(i+1) + (gtd/bn)*alpha*d(i+1) xnew(i+2) = x(i+2) + (gtd/bn)*alpha*d(i+2) xnew(i+3) = x(i+3) + (gtd/bn)*alpha*d(i+3) xnew(i+4) = x(i+4) + (gtd/bn)*alpha*d(i+4) end do c call evalfg(n,xnew,fnew,gnew, nexp) call dmsafg(nx,ny,xnew,fnew,gnew,'FG',bottom,top,left,right) fgcnt = fgcnt + 1 if(fgcnt .gt. maxfg) go to 999 end if * * *---------------------------------- Prepare some scalar products * * gtg = 0.d0 gtgp= 0.d0 ytg = 0.d0 stg = 0.d0 yts = 0.d0 do i = 1,n5 s(i) = xnew(i) - x(i) y(i) = gnew(i) - g(i) ytg = ytg + gnew(i)*y(i) stg = stg + gnew(i)*s(i) yts = yts + y(i)*s(i) gtgp = gtgp + gnew(i)*g(i) x(i) = xnew(i) g(i) = gnew(i) gtg = gtg + g(i) * g(i) end do do i = n6,n,5 s(i) = xnew(i) - x(i) s(i+1) = xnew(i+1) - x(i+1) s(i+2) = xnew(i+2) - x(i+2) s(i+3) = xnew(i+3) - x(i+3) s(i+4) = xnew(i+4) - x(i+4) y(i) = gnew(i) - g(i) y(i+1) = gnew(i+1) - g(i+1) y(i+2) = gnew(i+2) - g(i+2) y(i+3) = gnew(i+3) - g(i+3) y(i+4) = gnew(i+4) - g(i+4) ytg = ytg + gnew(i)*y(i)+gnew(i+1)*y(i+1)+gnew(i+2)*y(i+2)+ * gnew(i+3)*y(i+3)+gnew(i+4)*y(i+4) stg = stg + gnew(i)*s(i)+gnew(i+1)*s(i+1)+gnew(i+2)*s(i+2)+ * gnew(i+3)*s(i+3)+gnew(i+4)*s(i+4) yts = yts + y(i)*s(i)+y(i+1)*s(i+1)+y(i+2)*s(i+2)+ * y(i+3)*s(i+3)+y(i+4)*s(i+4) gtgp = gtgp + gnew(i)*g(i)+gnew(i+1)*g(i+1)+gnew(i+2)*g(i+2)+ * gnew(i+3)*g(i+3)+gnew(i+4)*g(i+4) x(i) = xnew(i) x(i+1) = xnew(i+1) x(i+2) = xnew(i+2) x(i+3) = xnew(i+3) x(i+4) = xnew(i+4) g(i) = gnew(i) g(i+1) = gnew(i+1) g(i+2) = gnew(i+2) g(i+3) = gnew(i+3) g(i+4) = gnew(i+4) gtg = gtg + g(i)*g(i)+g(i+1)*g(i+1)+g(i+2)*g(i+2)+ * g(i+3)*g(i+3)+g(i+4)*g(i+4) end do * gnorm= sqrt( gtg ) * f = fnew dnormprev = dnorm * * Sigma computation. * * sigma = gtg/(dabs(ytg)+gtg) * * * * -------------------------------- Delta, Theta and Beta computation * ================================= * * Delta computation. () * delta = stg*ytg-gtg*yts * acc (conjugacy condition) (ak) * bdc (descent condition) (bk) acc = cc*stg + ytg bdc = dc*yts*gtg + ytg*stg c * Beta and theta computation * If |delta| is greater than epsilon machine, then compute * beta and theta by means of the algorithm. * Otherwise, set theta=1 and beta = 0. * * Now theta computation: if(dabs(delta) .gt. epsm) then theta = (acc/ytg)*(1.d0+yts*gtg/delta)-bdc/delta else theta = 1.d0 itheta=itheta+1 end if c c * For beta computation we consider 3 possibilities: * 1) beta as in algorithm DESCON, * 2) beta as suggested by Dai and Liao+ (i.e. similar to PRP+), * 3) beta truncated (see Hager and Zhang). * * Therefore, the algorithm can be implemented in 4 variants as follows: * * a) (2.1),(2.2) and (2.12) for theta and (2.13) for beta computation. * * b) (2.1),(2.2) and (2.12) for theta and (7.29) for beta computation. * * c) (2.1),(2.2) and (2.12) for theta and (2.13) for beta comptation * with truncation. * * d) (2.1) (2.2) and (2.12) for theta and (7.29) for beta computation * with truncation. * * Please see ICI Technical Report, November 9, 2011. * * Numerical experiments proved that all these variants have similar * performances. * They do not improve significantly the performances of the algorithm. * However, variant a) is slightly better. c c if(yts .ne. 0.d0 .and. dabs(delta) .gt. epsm) then beta = ytg/yts-(bdc*ytg)/(yts*delta)+acc*gtg/delta cc beta = dmax1(ytg/yts,0.d0) - cc * (bdc*ytg)/(yts*delta) + cc * (acc*gtg)/delta cc beta = dmax1(beta,-1.d0/(dnorm*dmin1(0.1d0,gnorm))) else beta = 0.d0 ibeta=ibeta+1 end if * * --------------------------------- Direction computation * ======================= * dnorm = 0.0d0 gtd = 0.0d0 do i = 1,n5 d(i) = -theta*g(i) + beta*s(i) dnorm = dnorm + d(i) ** 2 gtd = gtd + g(i) * d(i) end do do i = n6,n,5 d(i) =-theta * g(i) + beta * s(i) d(i+1) =-theta * g(i+1) + beta * s(i+1) d(i+2) =-theta * g(i+2) + beta * s(i+2) d(i+3) =-theta * g(i+3) + beta * s(i+3) d(i+4) =-theta * g(i+4) + beta * s(i+4) dnorm = dnorm + d(i)**2+d(i+1)**2+d(i+2)**2+ * d(i+3)**2+d(i+4)**2 gtd = gtd + g(i)*d(i)+g(i+1)*d(i+1)+g(i+2)*d(i+2)+ * g(i+3)*d(i+3)+g(i+4)*d(i+4) end do dnorm = sqrt( dnorm ) * *------------------------------------ END direction computation * * RESTART CRITERIA * ================ * * Angle Restart Test * if(angle) then if ( gtd .gt. -1.0d-03 * gnorm * dnorm ) then irs = irs + 1 do i = 1,n d(i) = -g(i) end do dnorm = gnorm gtd = -gtg end if end if * * Beale-Powell restart test * if(powell) then if(dabs(gtgp) .gt. 0.2d0*dabs(gtg)) then irs = irs + 1 do i = 1,n5 d(i) = -g(i) end do do i = n6,n,5 d(i) = -g(i) d(i+1) = -g(i+1) d(i+2) = -g(i+2) d(i+3) = -g(i+3) d(i+4) = -g(i+4) end do dnorm = gnorm gtd = -gtg end if end if *------------------------------------------- Prepare first trial c of steplength if(dnorm .ne. 0.d0) then alpha = alpha * dnormprev / dnorm else alpha = 1.d0 end if c go to 110 *------------------------------------------ End of main loop 999 continue * return end *-------------------------------------------- END DESCON subroutine c****************************************************************** subroutine LineSearch (n,x,f,d,gtd,dnorm,alpha,xnew,fnew,gnew, + sigma,fgcnt,lscnt,lsflag, n5,n6,nexp, nx,ny) C This is the one-dimensional line search used in CONMIN C SCALAR ARGUMENTS integer n,fgcnt,lscnt,lsflag, nexp double precision f,gtd,dnorm,alpha,fnew C ARRAY ARGUMENTS double precision x(n),d(n),xnew(n),gnew(n) real*8 bottom(2000), top(2000), left(2000), right(2000) C LOCAL SCALARS integer i,lsiter double precision alphap,alphatemp,fp,dp,gtdnew,a,b,sigma lsflag = 0 * Maximum number of LineSearch is max$ls (now=20) max$ls=10 nexp=1 alphap = 0.0d0 fp = f dp = gtd do i = 1,n5 xnew(i) = x(i) + alpha * d(i) end do do i = n6,n,5 xnew(i) = x(i) + alpha * d(i) xnew(i+1) = x(i+1) + alpha * d(i+1) xnew(i+2) = x(i+2) + alpha * d(i+2) xnew(i+3) = x(i+3) + alpha * d(i+3) xnew(i+4) = x(i+4) + alpha * d(i+4) end do c1 c call evalfg(n,xnew,fnew,gnew, nexp) call dmsafg(nx,ny,xnew,fnew,gnew,'FG',bottom,top,left,right) fgcnt = fgcnt + 1 gtdnew = 0.0d0 do i = 1,n5 gtdnew = gtdnew + gnew(i) * d(i) end do do i = n6,n,5 gtdnew = gtdnew + gnew(i)*d(i)+gnew(i+1)*d(i+1)+gnew(i+2)*d(i+2)+ * gnew(i+3)*d(i+3)+gnew(i+4)*d(i+4) end do lsiter = 0 10 if ( alpha * dnorm .gt. 1.0d-30 .and. lsiter .lt. max$ls .and. + .not. ( gtdnew .eq. 0.0d0 .and. fnew .lt. f ) .and. + ( ( fnew .gt. f + 1.0d-04 * alpha * gtd .or. + dabs( gtdnew / gtd ) .gt. sigma ) .or. ( lsiter .eq. 0 .and. + dabs( gtdnew / gtd ) .gt. 0.5d0 ) ) ) then 20 if ( alpha * dnorm .gt. 1.0d-30 .and. fnew .gt. f .and. + gtdnew .lt. 0.0d0 ) then alpha = alpha / 3.0d0 do i = 1,n5 xnew(i) = x(i) + alpha * d(i) end do do i = n6,n,5 xnew(i) = x(i) + alpha * d(i) xnew(i+1) = x(i+1) + alpha * d(i+1) xnew(i+2) = x(i+2) + alpha * d(i+2) xnew(i+3) = x(i+3) + alpha * d(i+3) xnew(i+4) = x(i+4) + alpha * d(i+4) end do c2 c call evalfg(n,xnew,fnew,gnew, nexp) call dmsafg(nx,ny,xnew,fnew,gnew,'FG',bottom,top,left,right) fgcnt = fgcnt + 1 gtdnew = 0.0d0 do i = 1,n5 gtdnew = gtdnew + gnew(i) * d(i) end do do i = n6,n,5 gtdnew = gtdnew + gnew(i)*d(i)+gnew(i+1)*d(i+1)+ * gnew(i+2)*d(i+2)+gnew(i+3)*d(i+3)+ * gnew(i+4)*d(i+4) end do alphap = 0.0d0 fp = f dp = gtd goto 20 end if a = dp + gtdnew - 3.0d0 * ( fp - fnew ) / ( alphap - alpha ) b = a ** 2 - dp * gtdnew if ( b .gt. 0.0d0 ) then b = sqrt( b ) else b = 0.0d0 end if alphatemp = alpha - ( alpha - alphap ) * ( gtdnew + b - a ) / + ( gtdnew - dp + 2.0d0 * b ) if ( gtdnew / dp .le. 0.0d0 ) then if ( 0.99d0 * max( alpha, alphap ) .lt. alphatemp .or. + alphatemp .lt. 1.01d0 * min( alpha, alphap ) ) then alphatemp = ( alpha + alphap ) / 2.0d0 end if else if ( gtdnew .lt. 0.0d0 .and. + alphatemp .lt. 1.01d0 * max( alpha, alphap ) ) then alphatemp = 2.0d0 * max( alpha, alphap ) end if if ( ( gtdnew .gt. 0.0d0 .and. + alphatemp .gt. 0.99d0 * min( alpha, alphap ) ) .or. + alphatemp .lt. 0.0d0 ) then alphatemp = min( alpha, alphap ) / 2.0d0 end if end if alphap = alpha fp = fnew dp = gtdnew alpha = alphatemp do i = 1,n5 xnew(i) = x(i) + alpha * d(i) end do do i = n6,n,5 xnew(i) = x(i) + alpha * d(i) xnew(i+1) = x(i+1) + alpha * d(i+1) xnew(i+2) = x(i+2) + alpha * d(i+2) xnew(i+3) = x(i+3) + alpha * d(i+3) xnew(i+4) = x(i+4) + alpha * d(i+4) end do c3 c call evalfg(n,xnew,fnew,gnew, nexp) call dmsafg(nx,ny,xnew,fnew,gnew,'FG',bottom,top,left,right) fgcnt = fgcnt + 1 gtdnew = 0.0d0 do i = 1,n5 gtdnew = gtdnew + gnew(i) * d(i) end do do i = n6,n,5 gtdnew = gtdnew + gnew(i)*d(i)+gnew(i+1)*d(i+1)+ * gnew(i+2)*d(i+2)+gnew(i+3)*d(i+3)+ * gnew(i+4)*d(i+4) end do lsiter = lsiter + 1 goto 10 end if if ( lsiter .ge. max$ls ) then lsflag = 1 end if if ( lsiter .ne. 0 ) then lscnt = lscnt + 1 end if return end *---------------------------------- End LineSearch subroutine *----------------------------------------------------------- * Date created : May 30, 1995 * Date last modified : May 30, 1995 * * Subroutine for execution time computation. * *----------------------------------------------------------- * subroutine exetim(tih,tim,tis,tic, tfh,tfm,tfs,tfc) * integer*4 tih,tim,tis,tic integer*4 tfh,tfm,tfs,tfc * integer*4 ti,tf integer*4 ch,cm,cs data ch,cm,cs/360000,6000,100/ * ti=tih*ch+tim*cm+tis*cs+tic tf=tfh*ch+tfm*cm+tfs*cs+tfc tf=tf-ti tfh=tf/ch tf=tf-tfh*ch tfm=tf/cm tf=tf-tfm*cm tfs=tf/cs tfc=tf-tfs*cs * return end *---------------------------------------------- End of EXETIM * * * *** * ***** * *** * * c Function and its Gradient c subroutine dmsafg(nx,ny,x,f,fgrad,task,bottom,top,left,right) character*(*) task integer nx, ny double precision f double precision x(nx*ny), fgrad(nx*ny), bottom(nx+2), top(nx+2), + left(ny+2), right(ny+2) c ********** c c Subroutine dmsafg c c This subroutine computes the function and gradient of the c minimal surface area problem. c c The subroutine statement is c c subroutine dmsafg(nx,ny,x,f,fgrad,task,bottom,top,left,right) c c where c c nx is an integer variable. c On entry nx is the number of grid points in the first c coordinate direction. c On exit nx is unchanged. c c ny is an integer variable. c On entry ny is the number of grid points in the second c coordinate direction. c On exit ny is unchanged. c c x is a double precision array of dimension nx*ny. c On entry x specifies the vector x if task = 'F', 'G', or 'FG'. c Otherwise x need not be specified. c On exit x is unchanged if task = 'F', 'G', or 'FG'. Otherwise c x is set according to task. c c f is a double precision variable. c On entry f need not be specified. c On exit f is set to the function evaluated at x if task = 'F' c or 'FG'. c c fgrad is a double precision array of dimension nx*ny. c On entry fgrad need not be specified. c On exit fgrad contains the gradient evaluated at x if c task = 'G' or 'FG'. c c task is a character variable. c On entry task specifies the action of the subroutine: c c task action c ---- ------ c 'F' Evaluate the function at x. c 'G' Evaluate the gradient at x. c 'FG' Evaluate the function and the gradient at x. c 'XS' Set x to the standard starting point xs. c c On exit task is unchanged. c c bottom is a double precision array of dimension nx + 2. c On entry bottom must contain boundary data beginning c with the lower left corner of the domain. c On exit bottom is unchanged. c c top is a double precision array of dimension nx + 2. c On entry top must contain boundary data beginning with c the upper left corner of the domain. c On exit top is unchanged. c c left is a double precision array of dimension ny + 2. c On entry left must contain boundary data beginning with c the lower left corner of the domain. c On exit left is unchanged. c c right is a double precision array of dimension ny + 2. c On entry right must contain boundary data beginning with c the lower right corner of the domain. c On exit right is unchanged. c c MINPACK-2 Project. November 1993. c Argonne National Laboratory and University of Minnesota. c Brett M. Averick. c c ********** double precision one, p5, two, zero parameter (zero=0.0d0,p5=0.5d0,one=1.0d0,two=2.0d0) logical feval, geval integer i, j, k double precision alphaj, area, betai, dvdx, dvdy, fl, fu, hx, hy, + v, vb, vl, vr, vt, xline, yline c Initialize. hx = one/dble(nx+1) hy = one/dble(ny+1) area = p5*hx*hy c Compute the standard starting point if task = 'XS'. if (task .eq. 'XS') then do 20 j = 1, ny alphaj = dble(j)*hy do 10 i = 1, nx k = nx*(j-1) + i betai = dble(i)*hx yline = alphaj*top(i+1) + (one-alphaj)*bottom(i+1) xline = betai*right(j+1) + (one-betai)*left(j+1) x(k) = (yline+xline)/two 10 continue 20 continue return end if if (task .eq. 'F' .or. task .eq. 'FG') then feval = .true. else feval = .false. end if if (task .eq. 'G' .or. task .eq. 'FG') then geval = .true. else geval = .false. end if c Evaluate the function if task = 'F', the gradient if task = 'G', c or both if task = 'FG'. if (feval) f = zero if (geval) then do 30 k = 1, nx*ny fgrad(k) = zero 30 continue end if c Computation of the function and gradient over the lower c triangular elements. do 50 j = 0, ny do 40 i = 0, nx k = nx*(j-1) + i if (i .ge. 1 .and. j .ge. 1) then v = x(k) else if (j .eq. 0) v = bottom(i+1) if (i .eq. 0) v = left(j+1) end if if (i .lt. nx .and. j .gt. 0) then vr = x(k+1) else if (i .eq. nx) vr = right(j+1) if (j .eq. 0) vr = bottom(i+2) end if if (i .gt. 0 .and. j .lt. ny) then vt = x(k+nx) else if (i .eq. 0) vt = left(j+2) if (j .eq. ny) vt = top(i+1) end if dvdx = (vr-v)/hx dvdy = (vt-v)/hy fl = sqrt(one+dvdx**2+dvdy**2) if (feval) f = f + fl if (geval) then if (i .ge. 1 .and. j .ge. 1) + fgrad(k) = fgrad(k) - (dvdx/hx+dvdy/hy)/fl if (i .lt. nx .and. j .gt. 0) + fgrad(k+1) = fgrad(k+1) + (dvdx/hx)/fl if (i .gt. 0 .and. j .lt. ny) + fgrad(k+nx) = fgrad(k+nx) + (dvdy/hy)/fl end if 40 continue 50 continue c Computation of the function and the gradient over the upper c triangular elements. do 70 j = 1, ny + 1 do 60 i = 1, nx + 1 k = nx*(j-1) + i if (i .le. nx .and. j .gt. 1) then vb = x(k-nx) else if (j .eq. 1) vb = bottom(i+1) if (i .eq. nx+1) vb = right(j) end if if (i .gt. 1 .and. j .le. ny) then vl = x(k-1) else if (j .eq. ny+1) vl = top(i) if (i .eq. 1) vl = left(j+1) end if if (i .le. nx .and. j .le. ny) then v = x(k) else if (i .eq. nx+1) v = right(j+1) if (j .eq. ny+1) v = top(i+1) end if dvdx = (v-vl)/hx dvdy = (v-vb)/hy fu = sqrt(one+dvdx**2+dvdy**2) if (feval) f = f + fu if (geval) then if (i .le. nx .and. j .gt. 1) + fgrad(k-nx) = fgrad(k-nx) - (dvdy/hy)/fu if (i .gt. 1 .and. j .le. ny) + fgrad(k-1) = fgrad(k-1) - (dvdx/hx)/fu if (i .le. nx .and. j .le. ny) + fgrad(k) = fgrad(k) + (dvdx/hx+dvdy/hy)/fu end if 60 continue 70 continue c Scale the function and the gradient. if (feval) f = area*f if (geval) then do 80 k = 1, nx*ny fgrad(k) = area*fgrad(k) 80 continue end if end subroutine dmsabc(nx,ny,bottom,top,left,right) integer nx, ny double precision bottom(nx+2), top(nx+2), left(ny+2), right(ny+2) c ********** c c Subroutine dmsabc c c This subroutine computes Enneper's boundary conditions for the c minimal surface area problem on the unit square centered at the c origin. c c The subroutine statement is c c subroutine dmsabc(nx,ny,hx,hy,bottom,top,left,right) c c where c c nx is an integer variable. c On entry nx is the number of grid points in the first c coordinate direction. c On exit nx is unchanged. c c ny is an integer variable. c On entry ny is the number of grid points in the second c coordinate direction. c On exit ny is unchanged. c c bottom is a double precision array of dimension nx + 2. c On entry bottom need not be specified. c On exit bottom contains boundary values for the bottom c boundary of the domain. c c top is a double precision array of dimension nx + 2. c On entry top need not be specified. c On exit top contains boundary values for the top boundary of c the domain. c left is a double precision array of dimension ny + 2. c On entry left need not be specified. c On exit left contains boundary values for the left boundary c of the domain. c c right is a double precision array of dimension ny + 2. c On entry right need not be specified. c On exit right contains boundary values for the right boundary c of the domain. c c MINPACK-2 Project. November 1993. c Argonne National Laboratory and University of Minnesota. c Brett M. Averick. c c ********** integer maxit double precision b, l, one, r, t, three, tol, two parameter (one=1.0d0,two=2.0d0,three=3.0d0) parameter (maxit=5,tol=1.0d-10) parameter (b=-.50d0,t=.50d0,l=-.50d0,r=.50d0) integer i, j, k, limit double precision det, fnorm, hx, hy, xt, yt double precision nf(2), njac(2,2), u(2) c Compute Enneper's boundary conditions: bottom, top, left, then c right. Enneper's boundary values are obtained by defining c bv(x,y) = u**2 - v**2 where u and v are the unique solutions of c x = u + u*(v**2) - (u**3)/3, y = -v - (u**2)*v + (v**3)/3. hx = (r-l)/dble(nx+1) hy = (t-b)/dble(ny+1) do 40 j = 1, 4 if (j .eq. 1) then yt = b xt = l limit = nx + 2 else if (j .eq. 2) then yt = t xt = l limit = nx + 2 else if (j .eq. 3) then yt = b xt = l limit = ny + 2 else if (j .eq. 4) then yt = b xt = r limit = ny + 2 end if c Use Newton's method to solve xt = u + u*(v**2) - (u**3)/3, c yt = -v - (u**2)*v + (v**3)/3. do 30 i = 1, limit u(1) = xt u(2) = -yt do 10 k = 1, maxit nf(1) = u(1) + u(1)*u(2)**2 - u(1)**3/three - xt nf(2) = -u(2) - u(1)**2*u(2) + u(2)**3/three - yt fnorm = sqrt(nf(1)*nf(1)+nf(2)*nf(2)) if (fnorm .le. tol) go to 20 njac(1,1) = one + u(2)**2 - u(1)**2 njac(1,2) = two*u(1)*u(2) njac(2,1) = -two*u(1)*u(2) njac(2,2) = -one - u(1)**2 + u(2)**2 det = njac(1,1)*njac(2,2) - njac(1,2)*njac(2,1) u(1) = u(1) - (njac(2,2)*nf(1)-njac(1,2)*nf(2))/det u(2) = u(2) - (njac(1,1)*nf(2)-njac(2,1)*nf(1))/det 10 continue 20 continue if (j .eq. 1) then bottom(i) = u(1)*u(1) - u(2)*u(2) xt = xt + hx else if (j .eq. 2) then top(i) = u(1)*u(1) - u(2)*u(2) xt = xt + hx else if (j .eq. 3) then left(i) = u(1)*u(1) - u(2)*u(2) yt = yt + hy else if (j .eq. 4) then right(i) = u(1)*u(1) - u(2)*u(2) yt = yt + hy end if 30 continue 40 continue end C December 22, 2006 -- Neculai Andrei -- C ************* Last line for Minimal Surface Area Problem. C Neculai Andrei C Last line DESCON package C================================================================ * Last Line