MODULE ddgbundle_mod ! DDG-Bundle method USE r_precision, ONLY : prec ! Precision for reals. IMPLICIT NONE ! MODULE ddgbundle_mod includes the following subroutines (S) and functions (F). PUBLIC :: & init_ddgb, & ! S Initialization for DDG-bundle subroutine. ddgbundle ! S DDG-bundle subroutine for derivative free, nonsmooth, ! large-scale optimization. Contains: ! S restar Initialization and reinitialization. ! F sclpar Selection of scaling parameter. PRIVATE :: & lls_armijo, & ! S Armijo line search. lls_cubic, & ! S Wolfe line search. Contains: ! F qint Quadratic interpolation. ! F cint Cubic interpolation. aggre1, & ! S Simplified discrete gradient aggregation. aggre2, & ! S Discrete gradient aggregation. dgrad, & ! S Computation of discrete gradient wprint, & ! S Printout the error and warning messages. rprint ! S Printout the results. CONTAINS !*********************************************************************** !* * !* * SUBROUTINE init_ddgb * * !* * !* Initialization for diagonal discrete gradient bundle subroutine * !* for large-scale unconstrained nonsmooth optimization. * !* * !*********************************************************************** SUBROUTINE init_ddgb(iterm) USE param, ONLY : zero,half,one,small,large USE initializat, ONLY : & n, & ! Number of variables. mcu, & ! Maximum number of stored corrections, mcu >= 1. iprint, & ! Printout specification (see initializat for more details). iscale, & ! Selection of the diagonal matrix updating ! (see initializat for more details). tolf, & ! Tolerance for change of function values. tolf2, & ! Second tolerance for change of function values. tolb, & ! Tolerance for the function value. tolg, & ! Tolerance for the first termination criterion. told, & ! Tolerance for the diagonal elements. xmax, & ! Maximum stepsize, 1 < XMAX. eta, & ! Distance measure parameter, eta >= 0. epsl, & ! Line search parameter, 0 < epsl < 0.25. mtesf, & ! Maximum number of iterations with changes of ! function values smaller than tolf. mit, & ! Maximun number of iterations. mfe ! Maximun number of function evaluations. IMPLICIT NONE ! Scalar Arguments INTEGER, INTENT(OUT) :: iterm ! Cause of termination: ! 0 - Everything is ok. ! -5 - Invalid input parameters. ! Initialization and error checking IF (iprint < -1) iprint = 1 ! Printout specification. iterm = 0 IF (n <= 0) THEN iterm = -5 CALL wprint(iterm,iprint,1) RETURN END IF IF (epsl >= 0.25_prec) THEN iterm = -5 CALL wprint(iterm,iprint,4) RETURN END IF IF (mcu <= 0) THEN iterm = -5 CALL wprint(iterm,iprint,2) RETURN END IF ! Default values IF (mit <= 0) mit = 10000 ! Maximum number of iterations. IF (mfe <= 0) mfe = n*mit ! Maximum number of function evaluations. IF (tolf <= zero) tolf = 1.0E-08_prec ! Tolerance for change of function values. IF (tolf2 == zero) tolf2 = 1.0E+04_prec ! Second tolerance for change of function values. IF (tolb == zero) tolb = -large + small ! Tolerance for the function value. IF (tolg <= zero) tolg = 1.0E-06_prec ! Tolerance for the termination criterion. IF (told <= zero) told = 0.01_prec ! Tolerance for the diagonal elements. IF (xmax <= zero) xmax = 1.5_prec ! Maximum stepsize. IF (eta < zero) eta = half ! Distance measure parameter IF (epsl <= zero) epsl = 1.0E-04_prec ! Line search parameter, IF (mtesf <= 0) mtesf = 10 ! Maximum number of iterations with changes ! of function values smaller than tolf. IF (iscale > 6 .OR. iscale < 0) iscale = 4 ! Selection of the diagonal matrix updating method. END SUBROUTINE init_ddgb !*********************************************************************** !* * !* * SUBROUTINE ddgbundle * * !* * !* Diagonal discrete gradient bundle subroutine. * !* * !*********************************************************************** SUBROUTINE ddgbundle(f,nit,nfe,noutit,iterm,strtim) USE param, ONLY : small,large,zero,half,one ! Parameters. USE exe_time, ONLY : getime ! Execution time. USE initializat, ONLY : & n, & ! Number of variables. x, & ! Vector of variables mcu, & ! Maximum number of stored corrections, mcu >= 1. iprint, & ! Printout specification (see initializat for more details). iscale, & ! Selection of the diagonal matrix updating method ! (see initializat for more details). tolf, & ! Tolerance for change of function values. tolf2, & ! Second tolerance for change of function values. tolb, & ! Tolerance for the function value. tolg, & ! Tolerance for the first termination criterion. told, & ! Tolerance for the diagonal elements of diag. xmax, & ! Maximum stepsize, 1 < XMAX. eta, & ! Distance measure parameter, eta >= 0. epsl, & ! Line search parameter, 0 < epsl < 0.25. mtesf, & ! Maximum number of iterations with changes of ! function values smaller than tolf. mit, & ! Maximun number of iterations. mfe, & ! Maximun number of function evaluations. time ! Maximum time USE obj_fun, ONLY : & myf ! Computation of the value of the objective function. USE subpro, ONLY : & vdot, & ! Dot product. vxdiag, & ! Vector is multiplied by a diagonal matrix y:=diag*x. xdiffy, & ! Difference of two vectors. copy, & ! Copying of a vector. copy2 ! Copying of two vectors. IMPLICIT NONE ! Scalar Arguments REAL(KIND=prec), INTENT(OUT) :: & f ! Value of the objective function. INTEGER, INTENT(OUT) :: & noutit, & ! Number of outer iterations. nit,nfe ! Number of iterations and function evaluations. INTEGER, INTENT(OUT) :: & iterm ! Cause of termination: ! 1 - The problem has been solved with desired accuracy. ! 2 - Changes in function values < tolf in mtesf ! subsequent iterations. ! 3 - Changes in function value < tolf*small*MAX(|f_k|,|f_(k-1)|,1), ! where small is the smallest positive number such that ! 1.0 + small > 1.0. ! 4 - Number of function calls > mfe. ! 5 - Number of iterations > mit. ! 6 - Time limit exceeded. ! 7 - f < tolb. ! 8 - Failure in attaining the demanded accuracy. ! -1 - Two consecutive restarts. ! -2 - Number of restarts > maximum number of restarts. ! -3 - Failure in function calculations ! (assigned by the user). ! -4 ! -5 - Invalid input parameters. ! -6 - Unspecified error. ! Local Arrays REAL(KIND=prec), DIMENSION(n) :: & xo, & ! Previous vector of variables. dg, & ! Discrete gradient of the objective function. dgp, & ! Previous discrete gradient of the ohjective function. dga, & ! Aggregate discrete gradient. g, & ! Direction for discrete gradient. s, & ! Difference of current and previous variables. u, & ! Difference of current and previous discrete gradients. d, & ! Direction vector. diag ! Diagonal matrix. REAL(KIND=prec), DIMENSION(n*mcu) :: & ssmatrix, & ! Matrix whose columns are stored differences of ! variables^2 sumatrix ! Matrix whose columns are stored differences of ! variables and discrete gradients. ! Local Scalars REAL(KIND=prec) :: & tmpx, & delta, & ! Outer iteration parameter. sl, & ! Size parameter for discrete gradient. slout, & ! Bound for the size parameter sl. alfn, & ! Locality measure. alfv, & ! Aggregate locality measure. epsr, & ! Line search parameter. dnorm, & ! Euclidean norm of the direction vector. dgnorm, & ! Euclidean norm of the aggregate discrete gradient. xnorm, & ! Stopping criterion. pxnorm, & ! Previous stopping criterion. p, & ! Directional derivative. tmax, & ! Maximum stepsize. t, & ! Stepsize. theta, & ! Correction parameter for stepsize. fo, & ! Previous value of the objective. stu, & ! stu = trans(s)*u. utu, & ! utu = trans(u)*u. sts, & ! sts = trans(s)*s. gamma_old, & ! Previous scaling parameter. gamma, & ! Scaling parameter. tmp ! Auxiliary value. !hihu uusi INTEGER :: i,j, & ngrad, & ! Number of discrete gradient calls. mcc, & ! Current number of stored corrections. inew, & ! Index for the circular arrays. iters, & ! Null step indicator. ! 0 - Null step. ! 1 - Serious step. nnk, & ! Consecutive null steps counter. ibun, & ! Index for the circular arrays in bundle updating. mal, & ! Current size of the bundle. neps, & ! Number of consecutive equal stopping criterions. ntesf, & ! Number of tests on function decrease. ncres, & ! Number of restarts. nres, & ! Number of consecutive restarts. nout ! Auxilary printout specification. ! Intrinsic Functions INTRINSIC ABS,MAX,SQRT ! Computational Time REAL strtim,elapstim ! Parameters REAL(KIND=prec), PARAMETER :: & fmin = -large, & ! Smallest acceptable value of the function. lengthd = 1.0E+20_prec, & ! Direction vector length (no need if you use Armijo line-search). tmin = 1.0E-12_prec, & ! Minimum stepsize. slmin = 1.0E-5_prec, & ! Minimum size parameter for the discrete gradient. pwt = 1.0E-8_prec, & ! Parameter for discrete gradient, pwt = z(sl)*alpha^j*e_j. sigma = 1.0E-1_prec ! Acceleration multiplier for outer iteration parameter. INTEGER, PARAMETER :: & maxeps = 20, & ! Maximum number of consecutive equal stopping criterions. maxnrs = 2000 ! Maximum number of restarts. IF (iprint > 3) WRITE (6,FMT='(1X,''Entry to DDG-Bundle:'')') ! Initialization nout = 0 nit = 0 nfe = 0 ntesf = 0 nres = 1 ncres = -1 neps = 0 iterm = 0 iters = 1 nnk = 0 alfn = zero alfv = zero stu = one sts = one utu = one tmax = xmax xnorm = large epsr = 0.25_prec+small IF (epsl+epsl >= epsr) THEN epsr = epsl+epsl + small IF (epsr >= half) THEN CALL wprint(iterm,iprint,-2) END IF END IF ! Outer iteration tmp = one/SQRT(REAL(n, kind=prec)) ! ngrad = 0 noutit = 0 slout = 1.25E+0_prec delta = (REAL(n*n, kind=prec)) outer_iteration: DO WHILE (delta > tolg*sigma) ! iterm is the latest cause of stopping ! if it is positive everything is ok. sl = max(2.0E-1_prec * slout,slmin) ! slout = sl noutit = noutit + 1 iterm = 0 iters = 1 ntesf = 0 neps = 0 nres = 1 ! First direction for the discrete gradient DO i=1,n g(i)=tmp END DO ! Computation of the value and the discrete gradient of the objective ! function and the search direction for the first iteration CALL myf(n,x,f,iterm) nfe = nfe + 1 CALL dgrad(n,x,sl,g,dg,f,ngrad,nfe,pwt,iterm) IF (iterm /= 0) EXIT outer_iteration CALL restar(n,mcc,inew,ibun,iters,dgp,dg,nnk,alfv,alfn,gamma,d,diag,mal,ncres) ! Start of the inner iteration inner_iteration: DO ! Computational time IF (time > 0.0E+00) THEN CALL getime(elapstim) IF (elapstim-strtim > time) THEN iterm = 6 EXIT outer_iteration ! Time limit exeeded, real exit END IF END IF ! Computation of norms IF (iters > 0) THEN dgnorm = vdot(n,dg,dg) dnorm = SQRT(vdot(n,d,d)) p = vdot(n,dg,d) ELSE dgnorm = vdot(n,dga,dga) dnorm = SQRT(vdot(n,d,d)) p = vdot(n,dga,d) END IF ! Test on descent direction IF (p+small*SQRT(dgnorm)*dnorm <= zero) THEN nres = 0 ELSE nres = nres + 1 IF (nres == 1) THEN CALL wprint(iterm,iprint,-3) CALL restar(n,mcc,inew,ibun,iters,dgp,dg,nnk,alfv,alfn,gamma,d,diag,mal,ncres) IF (ncres > maxnrs) THEN nout = maxnrs iterm = -2 EXIT outer_iteration ! real exit END IF CYCLE inner_iteration END IF nout = -1 iterm = -1 EXIT outer_iteration ! real exit END IF ! Stopping criterion nit = nit + 1 pxnorm = xnorm xnorm = -p + 2.0_prec*alfv ! Tests for termination IF(xnorm <= delta) THEN iterm = 1 EXIT inner_iteration END IF IF (nfe >= mfe) THEN ! too many function calls nout = mfe iterm = 4 EXIT outer_iteration ! real exit END IF IF (nit >= mit) THEN ! too many iterations nout = mit iterm = 5 EXIT outer_iteration ! real exit END IF IF (f <= tolb) THEN ! unbounded minimum iterm = 7 EXIT outer_iteration ! real exit END IF IF (iters == 0) THEN IF (ABS(xnorm - pxnorm) <= small) THEN neps = neps + 1 IF (neps > maxeps) THEN iterm = 8 EXIT inner_iteration ! try again with different delta END IF ELSE neps = 0 END IF ELSE neps = 0 END IF IF (pxnorm < xnorm .AND. nnk > 2) THEN CALL wprint(iterm,iprint,-4) END IF CALL rprint(n,noutit,nit,nfe,ngrad,x,f,xnorm,half*dgnorm+alfv,iterm,iprint) ! Preparation of line search fo = f IF (iters > 0) THEN CALL copy2(n,x,xo,dg,dgp) END IF IF (dnorm > zero) tmax = MAX(xmax/dnorm,sl) ! Initial step size IF (iters > 0) THEN t=min(2.0_prec,tmax) ELSE t=min(1.0_prec,tmax) END IF ! Armijo line search CALL lls_armijo(n,x,dg,d,xo,t,fo,f,p,alfn,sl,pwt,tmin,dnorm,xnorm, & epsl,epsr,eta,iters,nfe,ngrad,nnk,iterm) ! Line search with cubic interpolation ! theta = one ! IF (dnorm > lengthd) THEN ! theta=lengthd/dnorm ! END IF ! CALL lls_cubic(n,x,dg,d,xo,t,fo,f,p,alfn,sl,pwt,tmin,dnorm,xnorm, & ! theta,epsl,epsr,eta,iters,nfe,ngrad,nnk,iterm) IF (iterm /= 0) EXIT outer_iteration ! real exit IF (tolf2 >= 0) THEN IF (ABS(fo-f) <= tolf2*small*MAX(ABS(f),ABS(fo),one) & .AND. iters == 1) THEN iterm = 3 EXIT inner_iteration END IF END IF IF (ABS(fo-f) <= tolf) THEN ntesf = ntesf + 1 IF (ntesf >= mtesf .AND. iters == 1) THEN iterm = 2 EXIT inner_iteration END IF ELSE ntesf = 0 END IF ! Computation of variables difference CALL xdiffy(n,x,xo,s) ! Computation of aggregate values and gradients difference IF (iters == 0) THEN nnk = nnk + 1 IF (nnk == 1) THEN CALL xdiffy(n,dg,dgp,u) CALL aggre1(n,dg,dgp,dga,u,d,diag,alfn,alfv,gamma) ELSE CALL xdiffy(n,dg,dgp,u) CALL aggre2(n,dg,dgp,dga,d,diag,alfn,alfv,gamma) END IF CALL copy(n,xo,x) f = fo ELSE nnk = 0 CALL xdiffy(n,dg,dgp,u) END IF ! Serious step initialization IF (iters > 0) THEN alfn = zero alfv = zero END IF ! Direction finding gamma_old = gamma IF (iters > 0) THEN ! Serious step update and direction determination SELECT CASE(iscale) CASE(4) ! Herskovits diagonal update CALL update(n,mcu,mcc,inew,s,u,ssmatrix,sumatrix) ! original with sts and stu DO i = 1, n stu = zero sts = zero DO j = 1, mcc stu = stu + sumatrix((j-1)*n+i) sts = sts + ssmatrix((j-1)*n+i) END DO IF (stu > small) THEN diag(i) = sts/stu ! IF (diag(i) < 1.0E-2_prec) diag(i) = 1.0E-2_prec IF (diag(i) < 1.0E-2_prec) diag(i) = 1.0E-2_prec IF (diag(i) > told) diag(i) = told ! IF (diag(i) < told) THEN ! diag(i) = told ! END IF ELSE diag(i) = one END IF END DO CALL vxdiag(n,diag,dg,d) ! d = -one*d gamma = one CASE(5) ! Complement of Herskovits's diagonal update CALL update(n,mcu,mcc,inew,u,s,ssmatrix,sumatrix) ! complement with utu and stu DO i = 1, n stu = zero sts = zero DO j = 1, mcc stu = stu + sumatrix((j-1)*n+i) sts = sts + ssmatrix((j-1)*n+i) END DO IF (sts > small) THEN ! utu > 0, complement with stu/utu diag(i) = stu/sts IF (diag(i) < 1.0E-6_prec) diag(i) = 1.0E-6_prec IF (diag(i) > told) diag(i) = told ! IF (diag(i) < told) THEN ! diag(i) = told ! END IF ELSE diag(i) = one END IF END DO CALL vxdiag(n,diag,dg,d) ! d = -one*d gamma = one CASE(0,2) ! scaling with stu/utu stu = vdot(n,s,u) IF (stu > zero) THEN ! pos.def. matrix utu = vdot(n,u,u) gamma = sclpar(iscale,sts,stu,utu) diag = gamma ELSE ! diag = gamma_old ! skipping ! gamma = gamma_old ! diag = one ! setting diag = I gamma = one ! END IF d = -gamma*dg ! diag = gamma*I CASE(1,3) ! scaling with sts/stu stu = vdot(n,s,u) IF (stu > zero) THEN ! pos.def. matrix sts = vdot(n,s,s) gamma = sclpar(iscale,sts,stu,utu) diag = gamma ELSE ! diag = gamma_old !skipping ! gamma = gamma_old ! diag = one ! setting diag = I gamma = one ! END IF d = -gamma*dg ! diag = gamma*I CASE DEFAULT ! no scaling: diag = I ! no need to update diag or gamma d = -one*dg ! END SELECT ELSE ! Update and direction determination for null steps. SELECT CASE(iscale) CASE(4,5) ! Herskovits's diagonal update CALL vxdiag(n,diag,dga,d) ! diag = diag at null steps. d = -one*d ! diag = one ! diag = 1 at null steps. ! d = -one*dga ! IF (nnk == 1) THEN ! the first null step ! CALL update(n,mcu,mcc,inew,s,u,ssmatrix,sumatrix) ! original with sts and stu ! DO i = 1, n ! stu = zero ! sts = zero ! DO j = 1, mcc ! stu = stu + sumatrix((j-1)*n+i) ! sts = sts + ssmatrix((j-1)*n+i) ! END DO ! IF (stu > small) THEN ! stu > 0, original with sts/stu ! diag(i) = sts/stu ! IF (diag(i) < told) diag(i) = told ! ELSE ! diag(i) = one ! you could also skip this ! END IF ! END DO ! END IF ! CALL vxdiag(n,diag,dga,d) ! d = -one*d ! stu = vdot(n,s,u) ! diag = gamma*I at null steps if gamma < gamma_old ! IF (stu > zero) THEN ! pos.def. matrix ! sts = vdot(n,s,s) ! gamma = sclpar(iscale,sts,stu,utu) ! CALL update(n,mcu,mcc,inew,s,u,ssmatrix,sumatrix) ! IF (gamma < gamma_old) THEN ! d = -gamma*dga ! diag = gamma ! ELSE ! CALL vxdiag(n,diag,dga,d) ! diag = diag ! d = -one*d ! gamma = gamma_old ! END IF ! ELSE ! stu <= 0 ! CALL vxdiag(n,diag,dga,d) ! diag = diag ! d = -one*d ! END IF CASE(0,2) ! diag = 1 at null steps ! diag = one ! d = -one*dga ! diag = diag at null steps (skipping the update) d = -gamma*dga ! no need to update diag or gamma ! diag = gamma*I at null steps if gamma < gamma_old ! stu = vdot(n,s,u) ! IF (stu > zero) THEN ! pos.def. matrix ! utu = vdot(n,u,u) ! gamma = sclpar(iscale,sts,stu,utu) ! ELSE ! gamma = gamma_old !! gamma = one ! END IF ! IF (gamma > gamma_old) gamma = gamma_old !! IF (gamma > gamma_old) gamma = one ! d = -gamma*dga ! diag = gamma CASE(1,3) ! diag = 1 at null steps ! diag = one ! d = -one*dga ! diag = diag at null steps (skipping the update) d = -gamma*dga ! no need to update diag or gamma ! diag = gamma*I at null steps if gamma < gamma_old ! stu = vdot(n,s,u) ! IF (stu > zero) THEN ! pos.def. matrix ! sts = vdot(n,s,s) ! gamma = sclpar(iscale,sts,stu,utu) ! ELSE ! gamma = gamma_old !! gamma = one ! END IF ! IF (gamma > gamma_old) gamma = gamma_old !! IF (gamma > gamma_old) gamma = one ! d = -gamma*dga ! diag = gamma CASE DEFAULT ! no scaling: diag = I d = -one*dga ! no need to update diag or gamma END SELECT END IF END DO inner_iteration delta = MIN(delta*sigma,xnorm) END DO outer_iteration ! Printout the final results IF (iprint > 3) WRITE (6,FMT='(1X,''Exit from DDG-Bundle:'')') CALL wprint(iterm,iprint,nout) CALL rprint(n,noutit,nit,nfe,ngrad,x,f,xnorm,half*dgnorm+alfv,iterm,iprint) CONTAINS SUBROUTINE restar(n,mcc,inew,ibun,iters,dgp,dg,nnk,alfv,alfn,gamma,d,diag,mal,ncres) ! Initialization ! USE param, ONLY : zero,one ! given in host ! USE subpro, ONLY : copy ! given in host IMPLICIT NONE ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & dgp ! Basic discrete gradient of the objective function. REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: & dg ! Current (auxiliary) discrete gradient of the objective function. REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: & d, & ! Search direction. diag ! Diagonal varaible metric matrix. ! Scalar Arguments INTEGER, INTENT(IN) :: & n ! Number of variables INTEGER, INTENT(OUT) :: & mcc, & inew, & ibun, & ! Index for the circular arrays in bundle updating. nnk, & ! Consecutive null steps counter. mal ! Current size of the bundle. INTEGER, INTENT(INOUT) :: & iters, & ! Null step indicator. ! 0 - Null step. ! 1 - Serious step. ncres ! Number of restarts. REAL(KIND=prec), INTENT(OUT) :: & alfn, & ! Locality measure. alfv, & ! Aggregate locality measure. gamma ! Scaling parameter. ! Restart mcc = 0 inew = 1 ibun = 1 mal = 0 ncres = ncres + 1 IF (iters == 0) THEN CALL copy(n,dgp,dg) iters = 1 nnk = 0 alfv=zero alfn=zero END IF ! gamma = large gamma = one d = -one*dg diag = one END SUBROUTINE restar FUNCTION sclpar(iscale,sts,stu,utu) RESULT(spar) ! Calculation of the scaling parameter for DDG-Bundle. ! USE param, ONLY : small,one,half ! given in host ! USE initializat, ONLY : told ! tolerance for the diagonal elements. IMPLICIT NONE ! Scalar Arguments REAL(KIND=prec), INTENT(IN) :: & sts, & ! sts = trans(s)*s. stu, & ! stu = trans(s)*u. utu ! utu = trans(u)*u. REAL(KIND=prec) :: & spar ! Scaling parameter. INTEGER, INTENT(IN) :: & iscale ! Selection of the scaling: ! 0 - Scaling at every iteration with STU/UTU. ! 1 - Scaling at every iteration with STS/STU. ! 2 - Interval scaling with STU/UTU. ! 3 - Interval scaling with STS/STU. ! 6 - No scaling. ! Intrinsic Functions INTRINSIC SQRT ! Computation of scaling parameter. SELECT CASE(iscale) ! Scaling parameter = STU/UTU CASE(0,2) IF (utu < SQRT(small)) THEN spar = one RETURN ELSE spar = stu/utu END IF ! Scaling parameter = STS/STU CASE(1,3,4,5) IF (stu < SQRT(small)) THEN spar = one RETURN ELSE spar = sts/stu END IF ! No scaling CASE DEFAULT spar = one RETURN END SELECT ! Scaling SELECT CASE(iscale) ! Lower limit CASE(0,1) IF (spar < one) spar = one ! ! Interval scaling CASE(2) ! IF (spar < 0.1_prec .OR. spar > 10.0_prec) spar = one ! IF (spar < half .OR. spar > 5.0_prec) spar = one ! CASE(3) ! IF (spar < 0.1 .OR. spar > 10.0_prec) spar = one ! IF (spar < half .OR. spar > 5.0_prec) spar = one ! ! Heskovits diagonal (null step) CASE(4) IF (spar < told ) spar = told ! Scaling at every iteration CASE DEFAULT CONTINUE END SELECT END FUNCTION sclpar END SUBROUTINE ddgbundle !************************************************************************ !* * !* * SUBROUTINE lls_armijo * * !* * !* Special line search. * !* * !************************************************************************ SUBROUTINE lls_armijo(n,x,dg,d,xo,t,fo,f,p,alfn,sl,pwt,tmin,dnorm,& wk,epsl,epsr,eta,iters,nfe,ngrad,nnk,iterm) USE param, ONLY : zero,half,one USE obj_fun, ONLY : & myf ! Computation of the value of the objective function. USE subpro, ONLY : & scsum, & ! Sum of a vector and the scaled vector. vdot ! Dot product of vectors. IMPLICIT NONE ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & d, & ! Direction vector. xo ! Previous vector of variables. REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: & x ! Vector of variables. REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: & dg ! Discrete gradient of the objective function. ! Scalar Arguments REAL(KIND=prec), INTENT(INOUT) :: & epsl,& ! Linesearch parameter. epsr,& ! Linesearch parameter. p, & ! Directional derivative. t ! Stepsize REAL(KIND=prec), INTENT(IN) :: & sl, & pwt, & eta, & ! Distance measure parameter. fo, & ! Previous value of the objective function. wk, & ! Stopping parameter. dnorm, & ! Euclidean norm of the direction vector. tmin ! Lower limit for stepsize parameter. REAL(KIND=prec), INTENT(OUT) :: & f, & ! Value of the objective function. alfn ! Locality measure. INTEGER, INTENT(IN) :: & n, & ! Number of variables nnk ! Number of consequtive null steps. INTEGER, INTENT(INOUT) :: & nfe, & ! Number of function evaluations. ngrad ! Number of discrete gradient calls. INTEGER, INTENT(OUT) :: & iters, & ! Null step indicator. ! 0 - Null step. ! 1 - Serious step. iterm ! Cause of termination: ! 0 - Everything is ok. ! -3 - Failure in function or subgradient calculations. ! Local arrays REAL(KIND=prec), DIMENSION(n) :: & g, & ! Initial direction. xu ! trial point. ! Local Scalars REAL(KIND=prec) :: & fu, & ! Auxiliary values of the objective function. alpha, & ! Increase parameter. a_orig, & ! Original value of increase parameter. epslwk ! Auxiliary scalar. INTEGER :: i,nin, & ! Number of interpolations. maxnin ! Intrinsic Functions INTRINSIC ABS,MAX ! Initialization nin = 0 IF (nnk >= 1) THEN maxnin = 40 ! maxnin = 50 ELSE maxnin = 40 ! maxnin = 50 END IF alpha = 2.0_prec a_orig = alpha epslwk = epsl*wk ! Function evaluation at a new point DO i=1,n g(i)=d(i)/dnorm END DO CALL scsum(n,t,d,xo,x) CALL myf(n,x,f,iterm) nfe = nfe + 1 IF (iterm /= 0) RETURN CALL scsum(n,alpha*t,d,xo,xu) CALL myf(n,xu,fu,iterm) nfe = nfe + 1 IF (iterm /= 0) RETURN IF (f <= fo - t*epslwk) THEN ! descent condition IF(fu > fo - alpha*t*epslwk) THEN ! Serious step iters=1 CALL dgrad(n,x,sl,g,dg,f,ngrad,nfe,pwt,iterm) RETURN ELSE ! Increase t inct: DO nin=nin+1 f=fu x=xu alpha = alpha * a_orig CALL scsum(n,alpha*t,d,xo,xu) CALL myf(n,xu,fu,iterm) nfe = nfe + 1 IF( fu > fo - alpha*t*epslwk .OR. nin > maxnin) THEN ! Serious step iters=1 CALL dgrad(n,x,sl,g,dg,f,ngrad,nfe,pwt,iterm) RETURN END IF END DO inct END IF ELSE alpha = one ! Decrease t dect: DO nin=nin+1 IF (nin > maxnin) EXIT dect alpha = alpha / a_orig ! IF (alpha*t < tmin) THEN alpha = alpha * a_orig EXIT dect END IF CALL scsum(n,alpha*t,d,xo,x) CALL myf(n,x,f,iterm) IF (f <= fo - alpha*t*epslwk) THEN ! descent condition ! Serious step iters=1 CALL dgrad(n,x,sl,g,dg,f,ngrad,nfe,pwt,iterm) RETURN END IF END DO dect END IF ! Null step iters = 0 CALL dgrad(n,x,sl,g,dg,f,ngrad,nfe,pwt,iterm) IF (iterm /= 0) RETURN p = vdot(n,dg,d) ! alfn = MAX(ABS(fo-f+p*t),eta*(t*dnorm)**2) alfn = MAX(ABS(fo-f+alpha*t*p),eta*(alpha*t*dnorm)**2) END SUBROUTINE lls_armijo !************************************************************************ !* * !* * SUBROUTINE lls_cubic * * !* * !* Special line search. * !* * !************************************************************************ SUBROUTINE lls_cubic(n,x,dg,d,xo,t,fo,f,p,alfn,sl,pwt,tmin,dnorm,& wk,theta,epsl,epsr,eta,iters,nfe,ngrad,nnk,iterm) USE param, ONLY : zero,half,one USE obj_fun, ONLY : & myf ! Computation of the value of the objective function. USE subpro, ONLY : & scsum, & ! Sum of a vector and the scaled vector. vdot ! Dot product of vectors. IMPLICIT NONE ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & d, & ! Direction vector. xo ! Previous vector of variables. REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: & x ! Vector of variables. REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: & dg ! Discrete gradient of the objective function. ! Scalar Arguments REAL(KIND=prec), INTENT(INOUT) :: & epsl,& ! Linesearch parameter. epsr,& ! Linesearch parameter. p, & ! Directional derivative. t ! Stepsize REAL(KIND=prec), INTENT(IN) :: & sl, & pwt, & theta, & ! Scaling parameter. eta, & ! Distance measure parameter. fo, & ! Previous value of the objective function. wk, & ! Stopping parameter. dnorm, & ! Euclidean norm of the direction vector. tmin ! Lower limit for stepsize parameter. REAL(KIND=prec), INTENT(OUT) :: & f, & ! Value of the objective function. alfn ! Locality measure. INTEGER, INTENT(IN) :: & n, & ! Number of variables nnk ! Number of consequtive null steps. INTEGER, INTENT(INOUT) :: & nfe, & ! Number of function evaluations. ngrad ! Number of discrete gradient calls. INTEGER, INTENT(OUT) :: & iters, & ! Null step indicator. ! 0 - Null step. ! 1 - Serious step. iterm ! Cause of termination: ! 0 - Everything is ok. ! -3 - Failure in function or subgradient calculations. ! Local arrays REAL(KIND=prec), DIMENSION(n) :: & g ! Initial direction. ! Local Scalars REAL(KIND=prec) :: & tl,tu, & ! Lower and upper limits for t used in interpolation. fl,fu, & ! Values of the objective function for t=tl and t=tu. pl, & ! Previous directional derivative. epsa,epst,epslk,epsrk, & ! Line search parameters. thdnorm,epsawk,epstwk,epslwk,epsrwk ! Auxiliary scalars. INTEGER :: i,nin ! Number of interpolations. ! Parameters INTEGER, PARAMETER :: maxint = 200 ! Maximum number of interpolations. ! Intrinsic Functions INTRINSIC ABS,MAX ! Initialization nin = 0 epst = epsl+epsl epsa = half*(epsr - epst) thdnorm = theta*dnorm pl = -p tl = zero tu = t fl = fo fu = fo IF (theta < one) THEN epst = theta*epst epsa = theta*epsa epslk = epsl epsl = theta*epsl epsrk = epsr epsr = theta*epsr END IF epsawk = epsa*wk epslwk = epsl*wk epsrwk = epsr*wk epstwk = epst*wk ! Function and gradient evaluation at a new point DO i=1,n g(i)=d(i)/dnorm END DO lls_iteration: DO CALL scsum(n,theta*t,d,xo,x) CALL myf(n,x,f,iterm) nfe = nfe + 1 IF (iterm /= 0) RETURN ! Null/descent step test (ITERS=0/1) iters = 1 IF (f <= fo - t*epstwk) THEN tl = t fl = f ELSE tu = t fu = f END IF ! Additional interpolation IF (f > fo - t*epslwk .AND. t > one) THEN nin=nin+1 IF (tl == zero .AND. wk > zero) THEN t = qint(tu,fl,fu,pl,one-half/(one-epst)) ! If you are using cubic interpolation ELSE t = half*(tu+tl) END IF CYCLE lls_iteration END IF ! original additional interpolation IF (f > fo - t*epslwk .AND. tu-tl >= tmin*0.1_prec & .AND. nnk >= 1 .AND. nin < maxint) THEN nin=nin+1 IF (tl == zero .AND. wk > zero) THEN t = qint(tu,fl,fu,pl,one-half/(one-epst)) ! If you are using cubic interpolation ELSE t = half*(tu+tl) END IF CYCLE lls_iteration END IF CALL dgrad(n,x,sl,g,dg,f,ngrad,nfe,pwt,iterm) IF (iterm /= 0) RETURN p = theta*vdot(n,dg,d) alfn = MAX(ABS(fo-f+p*t),eta*(t*thdnorm)**2) ! Serious step IF (f <= fo - t*epslwk .AND. (t >= tmin .OR. alfn > epsawk)) EXIT lls_iteration ! Null step IF (p-alfn >= -epsrwk .OR. tu-tl < tmin*0.1_prec .OR. & nin >= maxint) THEN iters = 0 EXIT lls_iteration END IF ! Interpolation nin=nin+1 IF (tl == zero .AND. wk > zero) THEN t = cint(tu,fl,fu,pl,p,one-half/(one-epst)) ELSE t = half*(tu+tl) END IF END DO lls_iteration IF (theta /= one) THEN epsl = epslk epsr = epsrk END IF CONTAINS FUNCTION qint(tu,fl,fu,wk,kappa) RESULT(t) ! Quadratic interpolation ! USE param, ONLY : half,one ! given in host IMPLICIT NONE ! Scalar Arguments REAL(KIND=prec), INTENT(IN) :: & fl, & ! Value of the objective function. fu, & ! Value of the objective function for t=tu. wk, & ! Directional derivative. tu, & ! Upper value of the stepsize parameter. kappa ! Interpolation parameter. REAL(KIND=prec) :: & t ! Stepsize. ! Local Scalars REAL(KIND=prec) :: tmp1,tmp2 ! Intrinsic Functions INTRINSIC MAX tmp1 = (fu-fl)/ (-wk*tu) ! Quadratic interpolation with one directional derivative tmp2 = 2.0_prec * (one - tmp1) IF (tmp2 > one) THEN ! Interpolation accepted t = MAX(kappa*tu,tu/tmp2) RETURN END IF ! Bisection t = half*tu END FUNCTION qint FUNCTION cint(tu,fl,fu,wk,pu,kappa) RESULT(t) ! Cubic interpolation ! USE param, ONLY : half,one ! given in host IMPLICIT NONE ! Scalar Arguments REAL(KIND=prec), INTENT(IN) :: & fl, & ! Value of the objective function. fu, & ! Value of the objective function for t=tu. wk, & ! Directional derivative. pu, & ! Directional derivative at t=tu. tu, & ! Upper value of the stepsize parameter. kappa ! Interpolation parameter. REAL(KIND=prec) :: & t ! Stepsize. ! Local Scalars REAL(KIND=prec) :: tmp1,tmp2,tmp3,tmp4,tmp5,tmp6 ! Intrinsic Functions INTRINSIC MAX,SQRT ! Cubic interpolation tmp1 = (fu-fl)/ (-wk*tu) tmp3 = pu/ (-wk) tmp4 = tmp3 - 2.0_prec*tmp1 + one tmp5 = tmp3 - 3.0_prec*tmp1 + 2.0_prec tmp6 = tmp5*tmp5 - 3.0_prec*tmp4 IF (tmp6 < zero) THEN t = half*tu ! bisection RETURN END IF tmp2 = tmp5 + SQRT(tmp6) IF (tmp2 > one) THEN ! Interpolation accepted t = MAX(kappa*tu,tu/tmp2) RETURN END IF ! Bisection t = half*tu END FUNCTION cint END SUBROUTINE lls_cubic !************************************************************************ !* * !* * SUBROUTINE update * * !* * !* Matrix update. * !* * !************************************************************************ SUBROUTINE update(n,mcu,mcc,inew,s,u,ssmatrix,sumatrix) USE subpro, ONLY : & copy2 ! Copying of two vectors. IMPLICIT NONE ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & s, & ! Difference of current and previous variables. u ! Difference of current and previous subgradients. REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: & ssmatrix, & ! Matrix whose columns are stored corrections. sumatrix ! Matrix whose columns are stored subgradient differences. ! Scalar Arguments INTEGER, INTENT(IN) :: & n, & ! Number of variables mcu ! Declared number of stored corrections. INTEGER, INTENT(INOUT) :: & mcc, & ! Current number of stored corrections. inew ! Index for circular arrays. ! Local Arrays REAL(KIND=prec), DIMENSION(n) :: & ss,su ! Componentwise multiplication of s*s and s*u. ! Local Scalars INTEGER :: i ! Compute new values DO i = 1, n ss(i)=s(i)*s(i) su(i)=s(i)*u(i) END DO ! Update ssmatrix and sumatrix CALL copy2(n,ss,ssmatrix((inew-1)*n+1:),su,sumatrix((inew-1)*n+1:)) inew = inew + 1 IF (inew > mcu) inew = 1 IF (mcc < mcu) mcc = mcc + 1 END SUBROUTINE update !************************************************************************ !* * !* * SUBROUTINE aggre1 * * !* * !* Computation of aggregate values after first null step. * !* * !************************************************************************ SUBROUTINE aggre1(n,dg,dgp,dga,u,d,diag,alfn,alfv,gamma) USE param, ONLY : zero,half,one USE subpro, ONLY : & vxdiag, & ! Vector multiplied by a diagonal matrix. vdot ! Dot product. IMPLICIT NONE ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & diag, & ! Diagonal variable metric update. d, & ! Direction vector. dg, & ! Current (auxiliary) subgradient of the objective function. dgp, & ! Previous subgradient of the objective function. u ! Difference of trial and aggregate gradients. REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: & dga ! Next aggregate subgradient of the objective function. ! Scalar Arguments REAL(KIND=prec), INTENT(OUT) :: & alfv ! Aggregate locality measure. REAL(KIND=prec), INTENT(IN) :: & gamma, & ! Scaling parameter. alfn ! Locality measure. INTEGER, INTENT(IN) :: & n ! Number of variables ! Local arrays REAL(KIND=prec), DIMENSION(n) :: tmpn ! Local Scalars REAL(KIND=prec) :: & p, & ! p = trans(d)*u - alfn. q, & ! q = trans(u)*dm*u, where dm is the inverse approximation of ! the Hessian calculated by using the L-BFGS formula. lam, & ! Multiplier used to calculate aggregate values. w ! Correction. INTEGER :: i ! Intrinsic Functions INTRINSIC MAX,MIN,SIGN ! Computation of p=trans(d)*u - alfn and the product q=trans(u)*dm*u p = vdot(n,d,u) - alfn CALL vxdiag(n,diag,u,tmpn) q = vdot(n,u,tmpn) lam = half + SIGN(half,p) IF (q > zero) lam = MIN(one,MAX(zero,p/q)) ! Computation of the aggregate values p = one - lam DO i=1,n dga(i)=lam*dg(i) + p*dgp(i) END DO alfv = lam*alfn END SUBROUTINE aggre1 !************************************************************************ !* * !* * SUBROUTINE aggre2 * * !* * !* Computation of aggregate values after consequtive null steps. * !* * !************************************************************************ SUBROUTINE aggre2(n,dg,dgp,dga,d,diag,alfn,alfv,gamma) USE param, ONLY : zero,one,small USE subpro, ONLY : & vxdiag, & ! Vector multiplied by a diagonal matrix. vdot, & ! Dot product. xdiffy ! Difference of two vectors. IMPLICIT NONE ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & diag, & ! Diagonal variable metric update. d, & ! Direction vector. dg, & ! Current (auxiliary) subgradient of the objective function. dgp ! Previous subgradient of the objective function. REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: & dga ! Aggregate subgradient of the objective function. ! Scalar Arguments REAL(KIND=prec), INTENT(INOUT) :: & alfv ! Aggregate locality measure. REAL(KIND=prec), INTENT(IN) :: & gamma, & ! Scaling parameter. alfn ! Locality measure. INTEGER, INTENT(IN) :: & n ! Number of variables ! Local arrays REAL(KIND=prec), DIMENSION(n) :: tmpn2,tmpn3,tmpn4 ! Local Scalars REAL(KIND=prec) :: & pr, & ! pr = trans(dgp-dga) dm (dgp-dga), where dm ! presents the L-SR1- approximation of Hessian. rrp, & ! rrp = trans(dgp-dga) dm dga - alfv. prqr, & ! prqr = trans(dgp-dga) dm (dg-dga). rrq, & ! rrq = trans(dg-dga) dm dga - alfv + alfn. qr, & ! qr = trans(dg-dga) dm (dg-dga). pq, & ! pq = trans(dg-dgp) dm (dg-dgp). qqp, & ! qqp = trans(dg-dgp) dm dg + alfn. lam1, & ! Multiplier used to calculate aggregate values. lam2, & ! Multiplier used to calculate aggregate values. tmp1, & ! Auxiliary scalar. tmp2, & ! Auxiliary scalar. tmp3 ! Auxiliary scalar. INTEGER :: i ! Intrinsic Functions INTRINSIC MIN,MAX CALL xdiffy(n,dgp,dga,tmpn2) ! Calculation of tmpn3 = trans(dgp - dga)dm CALL vxdiag(n,diag,tmpn2,tmpn3) pr = vdot(n,tmpn3,tmpn2) rrp = vdot(n,tmpn3,dga) CALL xdiffy(n,dg,dga,tmpn4) prqr = vdot(n,tmpn3,tmpn4) rrq = -vdot(n,tmpn4,d) ! calculation of qr = trans(dg - dga) dm (dg - dga) CALL vxdiag(n,diag,tmpn4,tmpn3) qr = vdot(n,tmpn4,tmpn3) pq = qr - prqr - prqr + pr qqp = pq + prqr + rrq - pr - rrp + alfn rrp = rrp - alfv rrq = rrq + alfn - alfv ! computation of multipliers lam1 and lam2 IF (pr > zero .AND. qr > zero) THEN tmp1 = rrq/qr tmp2 = prqr/qr tmp3 = pr - prqr*tmp2 IF (tmp3 /= zero) THEN lam1 = (tmp1*prqr - rrp)/tmp3 lam2 = -tmp1 - lam1*tmp2 IF (lam1*(lam1 - one) < zero .AND. & lam2*(lam1 + lam2 - one) < zero) GO TO 200 END IF END IF ! Minimum on the boundary 100 continue lam1 = zero lam2 = zero IF (alfn <= alfv) lam2 = one IF (qr > zero) lam2 = MIN(one,MAX(zero,-rrq/qr)) tmp3 = (lam2*qr + rrq+rrq)*lam2 tmp1 = zero IF (alfv >= zero) tmp1 = one IF (pr > zero) tmp1 = MIN(one,MAX(zero,-rrp/pr)) tmp2 = (tmp1*pr + rrp+rrp)*tmp1 IF (tmp2 < tmp3) THEN tmp3 = tmp2 lam1 = tmp1 lam2 = zero END IF IF (qqp*(qqp - pq) < zero) THEN IF (qr + rrq + rrq - qqp*qqp/pq < tmp3) THEN lam1 = qqp/pq lam2 = one - lam1 END IF END IF 200 CONTINUE IF (lam1 == zero .AND. lam2*(lam2 - one) < zero & .AND. -rrp - lam2*prqr > zero .AND. pr > zero) & lam1 = MIN(one - lam2, (-rrp-lam2*prqr)/pr) ! Computation of the aggregate values tmp1 = one - lam1 - lam2 DO i=1,n dga(i)=lam1*dgp(i)+lam2*dg(i)+tmp1*dga(i) END DO alfv = lam2*alfn + tmp1*alfv END SUBROUTINE aggre2 !************************************************************************ !* * !* * SUBROUTINE dgrad * * !* * !* Calculation of discrete gradients. * !* * !************************************************************************ SUBROUTINE dgrad(n,x,sl,g,dg,f,ngrad,nfe,pwt,iterm) USE param, ONLY : zero,half,one,small USE obj_fun, ONLY : myf ! Computation of the value of the objective function. USE subpro, ONLY : vdot ! Dot product IMPLICIT NONE ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & x, & ! Current iteration point. g ! Direction of discrete gradient. REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: & dg ! Discrete gradient. ! Scalar Arguments REAL(KIND=prec), INTENT(IN) :: & f, & ! Value of the objective at x. sl, & ! Parameter for discrete gradient. (lambda) pwt ! Parameter: pwt = z(sl)*alpha^j*e_j. INTEGER, INTENT(IN) :: & n ! Number of variables INTEGER, INTENT(INOUT) :: & ngrad, & ! Number of discrete gradient calls. nfe ! Number of function evaluations. INTEGER, INTENT(OUT) :: & iterm ! Cause of termination: ! 0 - Everything is ok. ! -3 - Failure in function calculations. ! Local arrays REAL(KIND=prec), DIMENSION(n) :: x1 ! Local Scalars REAL(KIND=prec) :: & a1, & ! a1 = argmax{|g_i|} r2, & ! Function value at (x + sl*g) r3, & ! Function value at (x^(j-1)) r4, & ! Function value at (x^j) dsum, & ! trans(dg)*g - dg(imax)*g(imax) tmp1 ! Auxiliary scalar. INTEGER :: k, & imax ! Index of the largest component of g. ! Intrinsic Functions INTRINSIC ABS ! Compute a1 = argmax{|g_i|} and set imax imax = 1 a1 = zero DO k=1,n tmp1 = ABS(g(k)) IF (a1 < tmp1) THEN a1 = tmp1 imax = k END IF END DO DO k=1,n x1(k) = x(k)+sl*g(k) END DO CALL myf(n,x1,r2,iterm) IF (iterm /= 0) RETURN nfe = nfe+1 dsum = zero r4 = r2 DO k=1,n IF (k /= imax) THEN r3 = r4 x1(k) = x1(k) + pwt CALL myf(n,x1,r4,iterm) nfe = nfe+1 IF (iterm /= 0) RETURN dg(k) = (r4 - r3)/pwt dsum = dsum + dg(k)*g(k) END IF END DO dg(imax) = (r2 - f)/(sl*g(imax)) - dsum/g(imax) ngrad = ngrad+1 END SUBROUTINE dgrad !************************************************************************ !* * !* * SUBROUTINE wprint * * !* * !* Printout the warning and error messages. * !* * !************************************************************************ SUBROUTINE wprint(iterm,iprint,nout) IMPLICIT NONE ! Scalar Arguments INTEGER, INTENT(IN) :: & iprint, & ! Printout specification: ! -1 - No printout. ! 0 - Only the error messages. ! 1 - The final values of the objective ! function. ! 2 - The final values of the objective ! function and the most serious ! warning messages. ! 3 - The whole final solution. ! 4 - At each iteration values of the ! objective function. ! 5 - At each iteration the whole ! solution nout, & ! Auxilary printout specification. iterm ! Cause of termination: ! 1 - The problem has been solved with desired accuracy. ! 2 - Changes in function values < tolf in mtesf ! subsequent iterations. ! 3 - Changes in function value < tolf*small*MAX(|f_k|,|f_(k-1)|,1), ! where small is the smallest positive number such that ! 1.0 + small > 1.0. ! 4 - Number of function calls > mfe. ! 5 - Number of iterations > mit. ! 6 - Time limit exceeded. ! 7 - f < tolb. ! 8 - Failure in attaining the demanded accuracy. ! -1 - Two consecutive restarts. ! -2 - Number of restarts > maximum number of restarts. ! -3 - Failure in function calculations (assigned by the user). ! -4 ! -5 - Invalid input parameters. IF (iprint >= 0) THEN ! Initial error messages IF (iterm == -5) THEN IF (nout == 1) WRITE (6,FMT='(1X,''Error: '' & ''Number of variables (n) is too small, iterm='',I3)') iterm IF (nout == 2) WRITE (6,FMT='(1X,''Error: '' & ''The maximum number of stored corrections (mcu) '' & ''is too small, iterm='',I3)') iterm IF (nout == 4) WRITE (6,FMT='(1X,''Error: '' & ''Line search parameter epsl >= 0.25, iterm='',I3)') iterm RETURN END IF ! Warning messages IF (iprint >= 2) THEN IF (iterm == 0) THEN IF (nout == -1) WRITE (6,FMT='(1X,''Warning: '' & ''mc > mcu. Assigned mc = mcu.'')') IF (nout == -2) WRITE (6,FMT='(1X,''Warning: '' & ''A line search parameter epsr >= 0.5.'')') IF (nout == -3) WRITE (6,FMT='(1X,''Warning: '' & ''A nondescent search direction occured. Restart.'')') IF (nout == -4) WRITE (6,FMT='(1X,''Warning: '' & ''Does not converge.'')') IF (nout == -5) WRITE (6,FMT='(1X,''Warning: '' & ''tmax < tmin. Restart.'')') RETURN END IF ! Printout the final results IF (iterm == 6) WRITE (6,FMT='(1X,''Abnormal exit: Time is up.'')') IF (iterm == 7) WRITE (6,FMT='(1X,''Abnormal exit: f < tolb.'')') IF (iterm == 2) WRITE (6,FMT='(1X,''Abnormal exit: '' & ''Too many steps without significant progress.'')') IF (iterm == 3) WRITE (6,FMT='(1X,''Abnormal exit: '' & ''The value of the function does not change.'')') IF (iterm == 5) WRITE (6,FMT='(1X,''Abnormal exit: '' & ''Number of iterations > '',I5)') nout IF (iterm == 4) WRITE (6,FMT='(1X,''Abnormal exit: '' & ''Number of function evaluations > '',I5)') nout IF (iterm == -1) THEN IF (nout == -1) THEN WRITE (6,FMT='(1X,''Abnormal exit: Two consecutive restarts.'')') ELSE WRITE (6,FMT='(1X,''Abnormal exit: tmax < tmin in two subsequent iterations.'')') END IF END IF IF (iterm == -2) WRITE (6,FMT='(1X,''Abnormal exit: '' & ''Number of restarts > '',I5''.'')') nout IF (iterm == -3) WRITE (6,FMT='(1X,''Abnormal exit: '' & ''Failure in function or subgradient calculations.'')') IF (iterm == 8) WRITE (6,FMT='(1X,''Abnormal exit: '' & ''Failure in attaining the demanded accuracy.'')') END IF END IF END SUBROUTINE wprint !************************************************************************ !* * !* * SUBROUTINE rprint * * !* * !* Printout the (final) results. * !* * !************************************************************************ SUBROUTINE rprint(n,noutit,nit,nfe,ngrad,x,f,wk,qk,iterm,iprint) IMPLICIT NONE ! Scalar Arguments INTEGER, INTENT(IN) :: & n, & ! Number of variables noutit, & ! Number of used outer iterations. nit, & ! Number of used iterations. nfe, & ! Number of used function evaluations. ngrad, & ! Number of used discrete gradient calls. iprint, & ! Printout specification: ! -1 - No printout. ! 0 - Only the error messages. ! 1 - The final values of the objective ! function. ! 2 - The final values of the objective ! function and the most serious ! warning messages. ! 3 - The whole final solution. ! 4 - At each iteration values of the ! objective function. ! 5 - At each iteration the whole ! solution iterm ! Cause of termination: ! 1 - The problem has been solved with desired accuracy. ! 2 - Changes in function values < tolf in mtesf ! subsequent iterations. ! 3 - Changes in function value < tolf*small*MAX(|f_k|,|f_(k-1)|,1), ! where small is the smallest positive number such that ! 1.0 + small > 1.0. ! 4 - Number of function calls > mfe. ! 5 - Number of iterations > mit. ! 6 - Time limit exceeded. ! 7 - f < tolb. ! 8 - Failure in attaining the demanded accuracy. ! -1 - Two consecutive restarts. ! -2 - Number of restarts > maximum number of restarts. ! -3 - Failure in function or subgradient calculations ! (assigned by the user). ! -4 ! -5 - Invalid input parameters. REAL(KIND=prec), INTENT(IN) :: & f, & ! Value of the objective function. wk, & ! Value of the first stopping criterion. qk ! Value of the second stopping criterion. ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & x ! Vector of variables ! Local Scalars INTEGER :: i ! Intermediate results IF (iterm == 0) THEN IF (iprint > 3) WRITE (6,FMT='(1X,''nout='',I3,2X,''nit='',I5,2X, & ''nfe='',I7,2X,''ndg='',I7,2X,''f='',D15.8,2X,''wk='',D11.4,2X, & ''qk='',D11.4,2X)') noutit,nit,nfe,ngrad,f,wk,qk IF (iprint == 5) WRITE (6,FMT='(1X,''x='', & 5D15.7:/(4X,5D15.7))')(x(i),i=1,n) RETURN END IF ! Final results IF (iprint > 0) WRITE (6,FMT='(1X,''nout='',I3,2X,''nit='',I5,2X, & ''nfe='',I7,2X,''ndg='',I7,2X,''f='',D15.8,2X,''wk='',D11.4,2X, & ''qk='',D11.4,2X,''iterm='',I3)') noutit,nit,nfe,ngrad,f,wk,qk,iterm IF (IPRINT .EQ. 3 .OR. IPRINT .EQ. 5) & WRITE (6,FMT='(1X,''x='',5D15.7:/(4X,5D15.7))')(x(i),i=1,n) END SUBROUTINE rprint END MODULE ddgbundle_mod