********************************************************** Line 1 C C December 2, 2011 C ========================================== C C C Main Program THREECG C ============================ C C Optimal Design with Composite Materials C =========================================== C C *** Application 4.4 from MINPACK-2 C *** B.M. Averick, R.G. Carter, J.J. More, G.L. Xue C *** The MINPACK-2 test problem collection. C *** Argonne National Laboratory, Preprint MCS-P153-0692 C *** June 1992. pp.36-39. C C C THREECG 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 THREECG implements an accelerated conjugate gradient algorithm C with three terms, that at each iteration both the descent and C the conjugacy conditions are guaranteed. C C THREECG implements an algorithm which is a modificatio of the C Hestenes and Stiefel conjugate gradient algorithm, or of that C of Hager and Zhang in such a way that the search direction is C descent and it satisfies the conjugacy condition. C These properties are independent of the line search. C C C The algorithm is described in: C N. Andrei, A simple three-term conjugate gradient algorithm for C unconstrained optimization. C Journal of Computational and Applied Mathematics, vol.241 (2013) C pp. 19-29. C C Remark: C This paper is dedicated to the memory of Neculai Gheorghiu C (1930-2009), my professor of Mathematical Analysis, Faculty of C Mathematics, "Alexandru Ioan Cuza" University, Iasi, Romania. C C------------------------------------------------------------------------- C C The algorithm is defined as: C C x(k+1) = x(k) + alpha*d(k) C C where alpha is computed by the Wolfe line search conditions. C C C The search direction of this algorithm is computed as: C C d(k+1) = -g(k+1) - delta(k)*s(k) - etha(k)*y(k) C C where: C C ||y(k)||^2 s(k)T*g(k+1) y(k)T*g(k+1) C delta(k) = (1 + ----------) * ------------ - ------------ , C y(k)T*s(k) y(k)T*s(k) y(k)T*s(k) C C C s(k)T*g(k+1) C etha(k) = ------------. C y(k)T*s(k) 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,iter,irs, fgcnt,lscnt,maxiter,maxfg integer itert, irstot, fgtot,lstot, stoptest integer*4 gh,gm,gs,gc, ght,gmt,gst,gct, timpexp integer nx, ny double precision epsg, epsm real*8 f, gnorm, timp, proc integer*2 iyear, imonth, iday character*40 numef logical angle, powell real*8 lambda common /acca/epsm * LOCAL ARRAYS double precision x(1000000), g(1000000) * Output files: open(unit=4,file='tha3.out',status='unknown') open(unit=5,file='tha3.rez',status='unknown') * Epsilon machine computation (Please see the common /acca/) epsm = 0.111e-30 *------------------------------------------------------------------- * Restart parameters: angle or Powell angle = .false. powell = .true. *------------------------------------------------------------------- * 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 *------------------------------------------------------------------------- * * Maximum number of iterations * Maximum number of function and gradient evaluations maxiter = 10000 maxfg = 15000 * Some convergence constant epsg = (10.d0)**(-6) *------------------------------------------------------------------------- write(4,15) 15 format(4x,'***************************************************',/, * 4x,'* THREECG ***',/, * 4x,'* A simple three-term conjugate gradient ***',/, * 4x,'* algorithm with Guaranteed Descent and ***',/, * 4x,'* Conjugacy conditions ***',/, * 4x,'* --------------------------------------------- ***',/, * 4x,'* Project: FCGA ***',/, * 4x,'* The Fastest Conjugate Gradient Algorithm ***',/, * 4x,'* ***',/, * 4x,'* Dr. Neculai Andrei ***',/, * 4x,'* Research Institute for Informatics ***',/, * 4x,'* Bucharest - Romania ***',/, * 4x,'***************************************************',/) call getdat(iyear, imonth, iday) write(4,21) imonth, iday, iyear 21 format(4x,'Date: --- Month:',i2,' Day:',i2,' Year: ',i4,/) if(powell) then write(4,22) 22 format(4x,'THREECG. Powell restart.') end if if(angle) then write(4,23) 23 format(4x,'THREECG. Angle restart.') end if write(4,26) 26 format(1x,91('-')) * *------------------------* *--------------------------------------------- | Here Start Experiments | * *------------------------* numef = 'Optimal Design with Composite Materials ' write(4,101) nexp, numef, stoptest write(*,101) nexp, numef, stoptest 101 format(/,5x,i2,4x,'THREECG Algorithm. ','Function:',a40,/, * 11x,'stoptest=',i2,/) * *---------------------------------------------------------------- * write(4,102) 102 format(7x,'n',3x,'iter',3x,'irs',2x,'fgcnt',1x,'lscnt', * 3x,'time(c)',8x,'fxnew',14x,'gnorm') write(4,26) * *---------------------------------------------------------------- * itert=0 irstot=0 fgtot=0 lstot=0 timp=0.d0 * *-------------------------------------- Dimension of the problem * lambda=0.008d0 nexp = 1 nx=1000 ny=1000 n=nx*ny * Initial guess call dodcfg(nx,ny,x,f,g,'XS',lambda) * Call the solver call gettim(gh,gm,gs,gc) call threecg(n,x,epsg,maxiter,maxfg,f,gnorm,stoptest, * iter,irs,fgcnt,lscnt, angle, powell,nexp) 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 * Write statistics * *---------------------------------------------------------- Write *.out write(4,201) n,iter,irs,fgcnt,lscnt,timpexp,f,gnorm 201 format(1x,i7,1x,i6,1x,i5,1x,i6,1x,i5,1x, * i9,e20.13,e20.13) *---------------------------------------------------------- Write *.rez if(n .eq. 1000) then write(5,301) nexp,n,iter,fgcnt,timpexp,f 301 format(i2,i6,2x,i6,2x,i6,2x,i6,2x,e20.13) else write(5,302) n,iter,fgcnt,timpexp,f 302 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(*,401) n, iter, irs, fgcnt, lscnt,timpexp,f,gnorm 401 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,501) itert, irstot, fgtot, lstot, timp/100.d0, proc write(*,502) itert, irstot, fgtot, lstot, timp/100.d0, proc 501 format(3x,'TOTAL',1x,i6,1x,i5,1x,i6,1x,i5,3x,f7.2,' (seconds)', * 4x,'proc= ',f6.2,'%',/) 502 format(3x,'TOTAL',1x,i6,1x,i5,1x,i6,1x,i5,3x,f7.2,' (seconds)', * 4x,'proc= ',f6.2,'%',/) * write(5,601) 601 format(1x,55('-')) write(5,602) 602 format(2x,' n',2x,' iter',' fgcnt time(c)',10x,'fx',/) stop end ********************************************************** END main program * *** Conjugate Gradient Algorithms Project *** * =============================================== * * * * *------------------------------------------------------------------- * Subroutine THREECG * ====================== * * * 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 * * * * /-----------------------------------------------------------------\ * | THREECG 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 THREECG 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. | * | | * |-----------------------------------------------------------------| * | | * | The calling sequence of THREECG is: | * | | * | subroutine threecg(n,x,epsg,maxiter,maxfg,f,gnorm,stoptest, | * | iter,irs,fgcnt,lscnt,angle, powell, | * | 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). | * |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. | * |-----------------------------------------------------------------| * | | * | | * |Calling subroutines: | * |==================== | * |Subroutine THREECG 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 | * | | * |-----------------------------------------------------------------| * | Neculai Andrei, 2011 | * \-----------------------------------------------------------------/ * * * ********************************************************************* subroutine threecg(n,x,epsg,maxiter,maxfg,f,gnorm,stoptest, * iter,irs,fgcnt,lscnt,angle, powell,nexp) parameter(ia=1000000) * SCALAR ARGUMENTS integer n,iter,irs,fgcnt,lscnt,maxiter,maxfg integer stoptest, nexp double precision epsg, f,gnorm logical angle, powell common /acca/epsm * ARRAY ARGUMENTS double precision x(n) * LOCAL SCALARS integer i,lsflag double precision fnew,alpha,gtg, + dnorm,dnormprev, ginf, + gtd, gtgp, stg, yty, epsm, + yts, ytg, etha, delta, bn,sigma, lambda integer nx, ny * LOCAL ARRAYS double precision xnew(ia),g(ia),gnew(ia),d(ia), + y(ia), s(ia) * Initialization c Optimal Design with Composite Materials c ======================================= lambda=0.008d0 nx=1000 ny=1000 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 dodcfg(nx,ny,x,f,g,'FG',lambda) 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 * Initial value of parameter sigma. sigma = 0.8d0 * -------------------------------- Main loop -------------------- *==================================================================== 110 continue *------------------------------------ STOP test section * ================= 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 * if(stoptest .eq. 2) then if(gnorm .le. epsg) go to 999 end if * 91 continue *---------------------------------- Increment iteration section * =========================== iter = iter + 1 if(iter .gt. maxiter) go to 999 *---------------------------------- Line search section * =================== * * 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) if(fgcnt .gt. maxfg) go to 999 * *---------------------------------- Acceleration section * ==================== * * (We use an/bn only if b is different from zero) * (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 dodcfg(nx,ny,xnew,fnew,gnew,'FG',lambda) 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 yty = 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) yty = yty + y(i)*y(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) yty = yty + y(i)*y(i)+y(i+1)*y(i+1)+y(i+2)*y(i+2)+ * y(i+3)*y(i+3)+y(i+4)*y(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 * * -------------------------------- Delta and etha computation * ========================== * * * Now delta computation: if(dabs(yts) .gt. epsm) then delta = (1.d0+yty/yts)*stg/yts-ytg/yts else delta = 1.d0 end if * * Now etha computation * if(yts .gt. epsm) then etha = stg/yts else etha = 0.d0 end if * * --------------------------------- Direction computation * ===================== * dnorm = 0.0d0 gtd = 0.0d0 do i = 1,n5 d(i) = -g(i) - delta*s(i) - etha*y(i) dnorm = dnorm + d(i) ** 2 gtd = gtd + g(i) * d(i) end do do i = n6,n,5 d(i) = -g(i) - delta * s(i) - etha * y(i) d(i+1) = -g(i+1) - delta * s(i+1) - etha * y(i+1) d(i+2) = -g(i+2) - delta * s(i+2) - etha * y(i+2) d(i+3) = -g(i+3) - delta * s(i+3) - etha * y(i+3) d(i+4) = -g(i+4) - delta * s(i+4) - etha * y(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 * of steplength * =================== if(dnorm .ne. 0.d0) then alpha = alpha * dnormprev / dnorm else alpha = 1.d0 end if * go to 110 *------------------------------------------ End of main loop * ================ 999 continue * return end *-------------------------------------------- END THREECG subroutine c****************************************************************** subroutine LineSearch (n,x,f,d,gtd,dnorm,alpha,xnew,fnew,gnew, + sigma,fgcnt,lscnt,lsflag, n5,n6,nexp) 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) C LOCAL SCALARS integer i,lsiter double precision alphap,alphatemp,fp,dp,gtdnew,a,b,sigma real*8 lambda integer nx, ny nexp=1 nx=1000 ny=1000 lambda=0.008d0 lsflag = 0 * Maximum number of LineSearch is max$ls (now=6) max$ls=8 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 dodcfg(nx,ny,xnew,fnew,gnew,'FG',lambda) 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.50d0 ) ) ) 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 dodcfg(nx,ny,xnew,fnew,gnew,'FG',lambda) 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 dodcfg(nx,ny,xnew,fnew,gnew,'FG',lambda) 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 * * * ***** * ********* * ***** * * *---------------------------------------------------------------- subroutine dodcfg(nx,ny,x,f,fgrad,task,lambda) character*(*) task integer nx, ny double precision f, lambda double precision x(nx*ny), fgrad(nx*ny) c ********** c c Subroutine dodcfg c c This subroutine computes the function and gradient of the c optimal design with composite materials problem. c c The subroutine statement is c c subroutine dodcfg(nx,ny,x,f,fgrad,task,lambda) 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 lambda is a double precision variable. c On entry lambda is the Lagrange multiplier. c On exit lambda is unchanged. c c Subprograms called c c MINPACK-supplied ... dodcps c c MINPACK-2 Project. November 1993. c Argonne National Laboratory and University of Minnesota. c Brett M. Averick. c c ********** double precision mu1, mu2, one, p5, two, zero parameter (zero=0.0d0,p5=0.5d0,one=1.0d0,two=2.0d0) parameter (mu1=one,mu2=two) logical feval, geval integer i, j, k double precision area, dpsi, dpsip, dvdx, dvdy, gradv, hx, hxhy, + hy, temp, t1, t2, v, vb, vl, vr, vt external dodcps c Initialization. hx = one/dble(nx+1) hy = one/dble(ny+1) hxhy = hx*hy area = p5*hxhy c Compute the break points. t1 = sqrt(two*lambda*mu1/mu2) t2 = sqrt(two*lambda*mu2/mu1) c Compute the standard starting point if task = 'XS'. if (task .eq. 'XS') then do 20 j = 1, ny temp = dble(min(j,ny-j+1))*hy do 10 i = 1, nx k = nx*(j-1) + i x(k) = -(min(dble(min(i,nx-i+1))*hx,temp))**2 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 the gradient over the lower c triangular elements. do 50 j = 0, ny do 40 i = 0, nx k = nx*(j-1) + i v = zero vr = zero vt = zero if (j .ge. 1 .and. i .ge. 1) v = x(k) if (i .lt. nx .and. j .gt. 0) vr = x(k+1) if (i .gt. 0 .and. j .lt. ny) vt = x(k+nx) dvdx = (vr-v)/hx dvdy = (vt-v)/hy gradv = dvdx**2 + dvdy**2 if (feval) then call dodcps(gradv,mu1,mu2,t1,t2,dpsi,0,lambda) f = f + dpsi end if if (geval) then call dodcps(gradv,mu1,mu2,t1,t2,dpsip,1,lambda) if (i .ge. 1 .and. j .ge. 1) + fgrad(k) = fgrad(k) - two*(dvdx/hx+dvdy/hy)*dpsip if (i .lt. nx .and. j .gt. 0) + fgrad(k+1) = fgrad(k+1) + two*(dvdx/hx)*dpsip if (i .gt. 0 .and. j .lt. ny) + fgrad(k+nx) = fgrad(k+nx) + two*(dvdy/hy)*dpsip 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 vb = zero vl = zero v = zero if (i .le. nx .and. j .gt. 1) vb = x(k-nx) if (i .gt. 1 .and. j .le. ny) vl = x(k-1) if (i .le. nx .and. j .le. ny) v = x(k) dvdx = (v-vl)/hx dvdy = (v-vb)/hy gradv = dvdx**2 + dvdy**2 if (feval) then call dodcps(gradv,mu1,mu2,t1,t2,dpsi,0,lambda) f = f + dpsi end if if (geval) then call dodcps(gradv,mu1,mu2,t1,t2,dpsip,1,lambda) if (i .le. nx .and. j .gt. 1) + fgrad(k-nx) = fgrad(k-nx) - two*(dvdy/hy)*dpsip if (i .gt. 1 .and. j .le. ny) + fgrad(k-1) = fgrad(k-1) - two*(dvdx/hx)*dpsip if (i .le. nx .and. j .le. ny) + fgrad(k) = fgrad(k) + two*(dvdx/hx+dvdy/hy)*dpsip end if 60 continue 70 continue c Scale the function. if (feval) f = area*f c Integrate v over the domain. if (feval) then temp = zero do 80 k = 1, nx*ny temp = temp + x(k) 80 continue f = f + hxhy*temp end if if (geval) then do 90 k = 1, nx*ny fgrad(k) = area*fgrad(k) + hxhy 90 continue end if end subroutine dodcps(t,mu1,mu2,t1,t2,result,option,lambda) integer option double precision t, mu1, mu2, t1, t2, result, lambda c ********** c c This subroutine computes the function psi(t) and the scaled c functions psi'(t)/t and psi''(t)/t for the optimal design c with composite materials problem. c c The subroutine statement is c c subroutine dodcps(t,mu1,mu2,t1,t2,result,option,lambda) c c where c c t is a double precision variable. c On entry t is the variable t c On exit t is unchanged c c mu1 is a double precision variable. c On entry mu1 is the reciprocal shear modulus of material 1. c On exit mu1 is unchanged. c c mu2 is a double precision variable. c On entry mu2 is the reciprocal shear modulus of material 2. c On exit mu2 is unchanged. c c t1 is a double precision variable. c On entry t1 is the first breakpoint. c On exit t1 is unchanged. c c t2 is a double precision variable. c On entry t2 is the second breakpoint. c On exit t2 is unchanged. c c result is a double precision variable. c On entry result need not be specified. c On exit result is set according to task. c c option is an integer variable. c On entry option specifies the action of the subroutine: c c if option = 0 then evaluate the function psi(t). c if option = 1 then evaluate the scaled function psi'(t)/t. c if option = 2 then evaluate the scaled function psi''(t)/t. c c On option task is unchanged. c c lambda is a double precision variable c On entry lambda is the Lagrange multiplier. c On exit lambda 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 p25, p5, zero parameter (zero=0.0d0,p25=0.25d0,p5=0.5d0) double precision sqrtt sqrtt = sqrt(t) if (option .eq. 0) then if (sqrtt .le. t1) then result = p5*mu2*t else if (sqrtt .gt. t1 .and. sqrtt .lt. t2) then result = mu2*t1*sqrtt - lambda*mu1 else if (sqrtt .ge. t2) then result = p5*mu1*t + lambda*(mu2-mu1) end if else if (option .eq. 1) then if (sqrtt .le. t1) then result = p5*mu2 else if (sqrtt .gt. t1 .and. sqrtt .lt. t2) then result = p5*mu2*t1/sqrtt else if (sqrtt .ge. t2) then result = p5*mu1 end if else if (option .eq. 2) then if (sqrtt .le. t1) then result = zero else if (sqrtt .gt. t1 .and. sqrtt .lt. t2) then result = -p25*mu2*t1/(sqrtt*t) else if (sqrtt .ge. t2) then result = zero end if end if end *====================================================== c END OPTIMAL DESIGN WITH COMPOSITE MATERIALS. C Neculai Andrei C Last line THREECG package C================================================================ * Last Line