MODULE ldgbm_mod ! Limited memory discrete gradient bundle method USE r_precision, ONLY : prec ! Precision for reals. IMPLICIT NONE ! MODULE ldgbm_mod includes the following subroutines (S) and functions (F). PUBLIC :: & init_ldgbm, & ! S Initialization for LDGBM subroutine. ldgbm ! S LDGBM subroutine for nonsmooth moderate and large- ! scale optimization. Contains: ! S restar Initialization and reinitialization. ! S dobun Bundle construction. PRIVATE :: & indic1, & ! S Initialization of indices. tinit, & ! S Calculation of initial step size. Contains: ! S destep Stepsize selection using polyhedral ! approximation for descent steps. ! S nulstep Stepsize selection using polyhedral ! approximation for null steps. ldgbm_lls, & ! S Line search. Contains: ! F qint Quadratic interpolation. ldgbm_dlbfgs, & ! S Computing the search direction by limited memory ! BFGS update. Contains: ! F ldgbm_sclpar Selection of scaling parameter. dlsr1, & ! S Computing the search direction by limited ! memory SR1 update. agbfgs, & ! S Simplified subgradient aggregation. aggsr1, & ! S Subgradient aggregation. ldgbm_wprint, & ! S Printout the error and warning messages. ldgbm_rprint ! S Printout the results. CONTAINS !*********************************************************************** !* * !* * SUBROUTINE init_ldgbm * * !* * !* Initialization for limited memory discrete gradient bundle * !* subroutine for unconstrained nonsmooth optimization. * !* * !*********************************************************************** SUBROUTINE init_ldgbm(mc,iterm) USE param, ONLY : zero,half,small,large USE initializat, ONLY : & n, & ! Number of variables. na, & ! Maximum bundle dimension, na >= 2. mcu, & ! Upper limit for maximum number of stored corrections, mcu >= 3. iprint, & ! Printout specification (see initializat for more details). iscale, & ! Selection of the scaling (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. tolg2, & ! Tolerance for the second termination criterion. 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(INOUT) :: mc ! Initial maximum number of stored corrections. 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 ldgbm_wprint(iterm,iprint,1) RETURN END IF IF (na < 2) THEN iterm = -5 CALL ldgbm_wprint(iterm,iprint,3) RETURN END IF IF (epsl >= 0.25_prec) THEN iterm = -5 CALL ldgbm_wprint(iterm,iprint,4) RETURN END IF IF (mcu <= 3) THEN iterm = -5 CALL ldgbm_wprint(iterm,iprint,2) RETURN END IF IF (mc > mcu) THEN mc = mcu CALL ldgbm_wprint(iterm,iprint,-1) END IF ! Default values IF (mc <= 0) mc = 3 ! Initial maximum number of corrections. IF (mit <= 0) mit = 10000 ! Maximum number of iterations. IF (mfe <= 0) mfe = mit*n ! 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 first termination criterion. IF (tolg2 <= zero) tolg2 = tolg ! Tolerance for the second termination criterion. 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 = 0 ! Selection of the scaling. END SUBROUTINE init_ldgbm !*********************************************************************** !* * !* * SUBROUTINE ldgbm * * !* * !* Limited memory discrete gradient bundle subroutine * !* for nonsmooth optimization. * !* * !*********************************************************************** SUBROUTINE ldgbm(mc,f,noutit,nit,nfe,iterm,strtim) USE param, ONLY : small,large,zero,half,one USE exe_time, ONLY : getime ! Execution time. USE initializat, ONLY : & n, & ! Number of variables. x, & ! Vector of variables na, & ! Maximum bundle dimension, na >= 2. mcu, & ! Upper limit for maximum number of stored corrections, mcu >= 3. iprint, & ! Printout specification (see initializat for more details). iscale, & ! Selection of the scaling (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. tolg2, & ! Tolerance for the second termination criterion. 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 lmbm_sub, ONLY : & vdot, & ! Dot product. vneg, & ! Copying of a vector with change of the sign. 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(INOUT) :: & mc ! Maximum number of stored corrections. 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. ! -1 - Two consecutive restarts. ! -2 - Number of restarts > maximum number of restarts. ! -3 - Failure in function or subgradient calculations ! (assigned by the user). ! -4 - Failure in attaining the demanded accuracy. ! -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. s, & ! Difference of current and previous variables. u, & ! Difference of current and previous discrete gradients. d, & ! Direction vector. tmpn1 ! Auxiliary array. REAL(KIND=prec), DIMENSION(n*(mcu+1)) :: & sm,um ! Matrises whose columns are stored differences of ! variables (sm) and discrete gradients (um). REAL(KIND=prec), DIMENSION((mcu+2)*(mcu+1)/2) :: & rm, & ! pper triangular matrix stored columnwise in the ! one-dimensional array. umtum ! Matrix whose columns are stored discrete gradient differences. REAL(KIND=prec), DIMENSION(mcu+1) :: & cm, & ! Diagonal matrix. smtgp, & ! smtgp = trans(sm)*dgp. umtgp, & ! umtgp = trans(um)*dgp. tmpmc1, & ! Auxiliary array. tmpmc2 ! Auxiliary array. REAL(KIND=prec), DIMENSION(n*na) :: & ax,ag ! Matrix whose columns are bundle points and discrete gradients. REAL(KIND=prec), DIMENSION(na) :: & af ! Vector of bundle values. ! Local Scalars REAL(KIND=prec) :: & delta, & ! Outer iteration parameter. sl, & ! Size parameter for discrete gradient. pwt, & ! Parameter for discrete gradient. sigma, & ! Acceleration multiplier for outer iteration parameter. alfn, & ! Locality measure. alfv, & ! Aggregate locality measure. epsr, & ! Line search parameter. dnorm, & ! Euclidean norm of the direction vector. gnorm, & ! Euclidean norm of the aggregate subgradient. 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. gamma, & ! Scaling parameter. slout, & ! tmp ! Auxiliary value. INTEGER :: i, & mcinit, & ! Initial maximum number of stored corrections. mcc, & ! Current number of stored corrections. inew, & ! Index for the circular arrays. ibfgs, & ! Index of the type of BFGS update. isr1, & ! Index of the type of SR1 update. 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. ic, & ! Correction indicator. icn, & ! Correction indicator for null steps. iflag, & ! Index for adaptive version. neps, & ! Number of consecutive equal stopping criterions. ntesf, & ! Number of tests on function decrease. ncres, & ! Number of restarts. nres, & ! Number of consecutive restarts. nress, & ! Number of consecutive restarts in case of tmax < tmin. nout, & ! Auxilary printout specification. ngrad ! Number of discrete gradient calls. ! Intrinsic Functions INTRINSIC ABS,MAX,SQRT ! Computational Time REAL strtim,elapstim ! Parameters REAL(KIND=prec), PARAMETER :: & fmin = -large, & ! Smallest acceptable value of the function. tmin = 1.0E-12_prec, & ! Minimum stepsize. slmin = 1.0E-6_prec, & ! Minimum size parameter for the discrete gradient. lengthd = 1.0E+20_prec, & ! Direction vector length. rho = 1.0E-12_prec ! Correction parameter. INTEGER, PARAMETER :: & maxeps = 10, & ! Maximum number of consecutive equal stopping criterions. maxnrs = 2000 ! Maximum number of restarts. IF (iprint > 3) WRITE (6,FMT='(1X,''Entry to LDGBM:'')') ! Initialization nout = 0 nit = 0 nfe = 0 ntesf = 0 nres = 1 ncres = -1 nress = 0 neps = 0 iterm = 0 iters = 1 nnk = 0 isr1 = 0 alfn = zero alfv = zero mcinit = mc tmax = xmax xnorm = large epsr = 0.25_prec+small IF (epsl+epsl >= epsr) THEN epsr = epsl+epsl + small IF (epsr >= half) THEN CALL ldgbm_wprint(iterm,iprint,-2) END IF END IF ! Outer iteration tmp = one/SQRT(REAL(n, kind=prec)) ! ngrad = 0 noutit = 0 pwt = 1.0E-08_prec ! slout = 1.25E+0_prec sigma = 1.0E-1_prec ! ! sigma = 2.0E-1_prec ! delta = 1.0E+6_prec ! ! delta = n*n 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 nress = 0 ! First direction for the discrete gradient DO i=1,n tmpn1(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,tmpn1,dg,f,ngrad,nfe,pwt,iterm) IF (iterm /= 0) EXIT outer_iteration CALL restar(n,mc,mcc,mcinit,inew,ibun,ibfgs,iters,dgp,dg,nnk, & alfv,alfn,gamma,d,ic,icn,mal,ncres,iflag) CALL dobun(n,na,mal,x,dg,f,ax,ag,af,iters,ibun) ! 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 gnorm = vdot(n,dg,dg) dnorm = SQRT(vdot(n,d,d)) p = vdot(n,dg,d) ELSE gnorm = 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(gnorm)*dnorm <= zero) THEN nres = 0 ELSE nres = nres + 1 IF (nres == 1) THEN CALL ldgbm_wprint(iterm,iprint,-3) CALL restar(n,mc,mcc,mcinit,inew,ibun,ibfgs,iters,dgp,dg,nnk, & alfv,alfn,gamma,d,ic,icn,mal,ncres,iflag) IF (ncres > maxnrs) THEN nout = maxnrs iterm = -2 EXIT outer_iteration ! real exit END IF CALL dobun(n,na,mal,x,dg,f,ax,ag,af,iters,ibun) 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 + alfv + alfv ! Tests for termination IF (xnorm <= 1.0E+03_prec*delta .AND. (mcc > 0 .OR. ibfgs == 2)) THEN IF(xnorm <= delta) THEN ! IF(half*gnorm + alfv <= delta .AND. xnorm <= delta) THEN ! More accurate results iterm = 1 EXIT inner_iteration END IF IF (mc < mcu .AND. iflag == 0) THEN mc = mc+1 iflag = 1 END IF END IF IF (nfe >= mfe) THEN nout = mfe iterm = 4 EXIT outer_iteration ! real exit END IF IF (nit >= mit) THEN nout = mit iterm = 5 EXIT outer_iteration ! Real exit END IF IF (f <= tolb) THEN iterm = 7 EXIT outer_iteration ! Unbounded minima, real exit END IF IF (iters == 0) THEN IF (ABS(xnorm - pxnorm) <= small) THEN neps = neps + 1 IF (neps > maxeps) THEN iterm = -4 EXIT inner_iteration ! try again with different delta END IF ELSE neps = 0 END IF ELSE neps = 0 END IF ! Correction IF (-p < rho*gnorm .OR. icn == 1) THEN xnorm = xnorm + rho*gnorm dnorm = SQRT(dnorm*dnorm-2.0_prec*rho*p+rho*rho*gnorm) IF (iters > 0) THEN DO i=1,n d(i) = d(i)-rho*dg(i) END DO ELSE DO i=1,n d(i) = d(i)-rho*dga(i) END DO icn = 1 END IF ic = 1 ELSE ic = 0 END IF IF (pxnorm < xnorm .AND. nnk > 2) THEN CALL ldgbm_wprint(iterm,iprint,-4) END IF CALL ldgbm_rprint(n,noutit,nit,nfe,ngrad,x,f,xnorm,half*gnorm+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 CALL tinit(n,na,mal,x,af,ag,ax,ibun,d,f,p,t,tmax,tmin, & eta,iters,iterm) IF (iterm /= 0) EXIT outer_iteration ! Real exit ! Line search with directional derivatives which allows null steps theta = one IF (dnorm > lengthd) THEN theta=lengthd/dnorm END IF CALL ldgbm_lls(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 ! Bundle updating CALL dobun(n,na,mal,x,dg,f,ax,ag,af,iters,ibun) ! 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 copy(n,dgp,tmpn1) CALL xdiffy(n,dg,dgp,u) CALL agbfgs(n,mc,mcc,inew,ibfgs,iflag,dg,dgp,dga,u,d,sm,um, & rm,cm,umtum,alfn,alfv,gamma,ic,rho) ELSE CALL copy(n,dga,tmpn1) CALL aggsr1(n,mc,mcc,inew,iflag,dg,dgp,dga,d,alfn,alfv & ,umtum,rm,gamma,smtgp,umtgp,tmpmc1,tmpmc2,sm & ,um,icn,rho) CALL xdiffy(n,dg,dgp,u) END IF CALL copy(n,xo,x) f = fo ELSE IF (nnk /= 0) THEN CALL copy(n,dga,tmpn1) ELSE CALL copy(n,dgp,tmpn1) END IF nnk = 0 CALL xdiffy(n,dg,dgp,u) END IF ! Serious step initialization IF (iters > 0) THEN sl = MAX(slmin,half*sl) ! icn = 0 alfn = zero alfv = zero END IF ! Direction finding IF (iters > 0) THEN ! BFGS update and direction determination CALL ldgbm_dlbfgs(n,mc,mcc,inew,ibfgs,iflag,d,dg,dgp,s,u,sm,um,rm, & umtum,cm,smtgp,umtgp,gamma,tmpn1,iscale) ELSE ! SR1 update and direction determination CALL dlsr1(n,mc,mcc,inew,isr1,iflag,d,dgp,dga,s,u,sm,um,rm, & umtum,cm,smtgp,umtgp,gamma,tmpmc1,tmpmc2,tmpn1,nnk,iprint) ibfgs=0 END IF END DO inner_iteration delta = MIN(delta*sigma,xnorm) END DO outer_iteration ! Printout the final results IF (iprint > 3) THEN WRITE (6,FMT='(1X,''Exit from LDGBM:'')') END IF CALL ldgbm_wprint(iterm,iprint,nout) CALL ldgbm_rprint(n,noutit,nit,nfe,ngrad,x,f,xnorm,half*gnorm+alfv,iterm,iprint) CONTAINS SUBROUTINE restar(n,mc,mcc,mcinit,inew,ibun,ibfgs,iters,gp,g,nnk, & alfv,alfn,gamma,d,ic,icn,mal,ncres,iflag) ! Initialization ! USE param, ONLY : zero,one ! given in host ! USE lmbm_sub, ONLY : vneg,copy ! given in host IMPLICIT NONE ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & gp ! Basic subgradient of the objective function. REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: & g ! Current (auxiliary) subgradient of the objective function. REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: & d ! Search direction. ! Scalar Arguments INTEGER, INTENT(IN) :: & n, & ! Number of variables mcinit ! Initial maximum number of stored corrections. INTEGER, INTENT(OUT) :: & mc, & ! Current maximum number of stored corrections. mcc, & ! Current number of stored corrections. inew, & ! Index for the circular arrays. ibun, & ! Index for the circular arrays in bundle updating. ibfgs, & ! Index of the type of BFGS update. nnk, & ! Consecutive null steps counter. ic, & ! Correction indicator. icn, & ! Correction indicator for null steps. mal, & ! Current size of the bundle. iflag ! Index for adaptive version. 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 mc = mcinit mcc = 0 inew = 1 ibun = 1 ibfgs = 0 ic = 0 icn = 0 mal = 0 ncres = ncres + 1 iflag = 0 IF (iters == 0) THEN CALL copy(n,gp,g) iters = 1 nnk = 0 alfv=zero alfn=zero END IF gamma = one CALL vneg(n,g,d) END SUBROUTINE restar SUBROUTINE dobun(n,ma,mal,x,g,f,ay,ag,af,iters,ibun) ! Bundle construction ! USE lmbm_sub, ONLY : copy2 ! given in host IMPLICIT NONE ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & x, & ! Vector of variables g ! Subgradient of the objective function. REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: & ay, & ! Matrix whose columns are bundle points ! (stored in one-dimensional n*ma -array). ag, & ! Matrix whose columns are bundle subgradients ! (stored in one-dimensional n*ma -array). af ! Vector of values of bundle functions ! (stored in one-dimensional ma -array). ! Scalar Arguments INTEGER, INTENT(IN) :: & n, & ! Number of variables iters, & ! Null step indicator. ! 0 - Null step. ! 1 - Serious step. ma ! Maximum size of the bundle. INTEGER, INTENT(INOUT) :: & ibun, & ! Index for the circular arrays in bundle updating. mal ! Current size of the bundle. REAL(KIND=prec), INTENT(IN) :: & f ! Value of the objective function. ! Local Scalars INTEGER :: i,j IF (iters == 1) THEN ! Serious step af(ibun) = f i = (ibun-1)*n + 1 CALL copy2(n,g,ag(i:),x,ay(i:)) ELSE ! Null step IF (mal < ma) THEN af(ibun) = af(mal) af(mal) = f i = mal*n + 1 CALL copy2(n,ag(i-n:),ag(i:),ay(i-n:),ay(i:)) CALL copy2(n,g,ag(i-n:),x,ay(i-n:)) ELSE i = ibun-1 IF (i < 1) i = mal af(ibun) = af(i) af(i) = f i = (ibun-2)*n + 1 IF (i < 1) i = (mal-1)*n + 1 j = (ibun-1)*n + 1 CALL copy2(n,ag(i:),ag(j:),ay(i:),ay(j:)) CALL copy2(n,g,ag(i:),x,ay(i:)) END IF END IF mal = mal + 1 IF (mal > ma) mal = ma ibun = ibun + 1 IF (ibun > ma) ibun = 1 END SUBROUTINE dobun END SUBROUTINE ldgbm !************************************************************************ !* * !* * SUBROUTINE tinit * * !* * !* Initial stepsize selection for limited memory bundle method * !* * !************************************************************************ SUBROUTINE tinit(n,na,mal,x,af,ag,ay,ibun,d,f,p,t,tmax,tmin,eta,iters,iterm) USE param, ONLY : zero,half,one,large ! half and large are used in destep and nulstep IMPLICIT NONE ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & x, & ! Vector of variables (n array). d, & ! Direction vector (n array). ay, & ! Matrix whose columns are bundle points (stored in one-dimensional n*ma -array). ag ! Matrix whose columns are bundle subgradients (stored in one-dimensional n*ma -array). REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: & af ! Vector of values of bundle functions (stored in one-dimensional ma -array). ! Scalar Arguments REAL(KIND=prec), INTENT(OUT) :: & t ! Initial stepsize REAL(KIND=prec), INTENT(IN) :: & p, & ! Directional derivative. eta, & ! Distance measure parameter. f, & ! Value of the objective function. tmax, & ! Upper limit for stepsize parameter. tmin ! Lower limit for stepsize parameter. INTEGER, INTENT(IN) :: & n, & ! Number of variables na, & ! Maximum size of the bundle. mal, & ! Current size of the bundle. ibun, & ! Index for the circular arrays in bundle updating. iters ! Null step indicator. ! 0 - Null step. ! 1 - Serious step. INTEGER, INTENT(OUT) :: & iterm ! Cause of termination: ! 0 - Everything is ok. ! -6 - Error. ! INTEGER :: koe = 1 ! Intrinsic Functions INTRINSIC MAX,MIN t = MIN(one,tmax) IF (p == zero) RETURN IF (iters == 1) THEN CALL destep(n,na,mal,x,af,ag,ay,ibun,d,f,p,t,eta,iterm) ELSE CALL nulstep(n,na,mal,x,af,ag,ay,ibun,d,f,p,t,eta,iterm) END IF t = MIN(MAX(t,tmin),tmax) CONTAINS SUBROUTINE destep(n,ma,mal,x,af,ag,ay,ibun,d,f,df,t,eta,iterm) ! Stepsize selection ! USE param, ONLY : zero,half,one,large ! given in host IMPLICIT NONE ! Scalar Arguments REAL(KIND=prec), INTENT(INOUT) :: & t ! Initial stepsize REAL(KIND=prec), INTENT(IN) :: & df, & ! Directional derivative. eta, & ! Distance measure parameter. f ! Value of the objective function. INTEGER, INTENT(IN) :: & n, & ! Number of variables ma, & ! Maximum size of the bundle. mal, & ! Current size of the bundle. ibun ! Index for the circular arrays in bundle updating. INTEGER, INTENT(OUT) :: & iterm ! Cause of termination: ! 0 - Everything is ok. ! -6 - Error. ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & x, & ! Vector of variables (n array). d, & ! Direction vector (n array). ay, & ! Matrix whose columns are bundle points (stored in one-dimensional n*ma -array). ag, & ! Matrix whose columns are bundle subgradients (stored in one-dimensional n*ma -array). af ! Vector of values of bundle functions (stored in one-dimensional ma -array). ! Local Arrays REAL(KIND=prec), DIMENSION(2*ma) :: & tmparray ! Auxiliary array. ! Local Scalars REAL(KIND=prec) :: alf,alfl,alfr,bet,betl,betr,dx,q,r,w INTEGER :: i,j,jn,k,l,lq,ib ! Intrinsic Functions INTRINSIC ABS,REAL,MAX,MIN,SQRT iterm = 0 alfl = zero betl = zero w = df*t* (one - t*half) ! Initial choice of possibly active lines k = 0 l = -1 jn = (ibun-1)*n betr = - large DO j=1,mal-1 ib = ibun - 1 + j IF (ib > mal) ib = ib - mal IF (jn >= mal*n) jn = jn - mal*n r = zero bet = zero alfl = af(ib) - f DO i=1,n dx = x(i) - ay(jn+i) q = ag(jn+i) r = r + dx*dx alfl = alfl + dx*q bet = bet + d(i)*q END DO alf = MAX(ABS(alfl),eta*r) r = one - bet/df IF (r*r + (alf+alf)/df > 1.0E-6_prec) THEN k = k + 1 tmparray(k) = alf tmparray(ma+k) = bet r = t*bet - alf IF (r > w) THEN w = r l = k END IF END IF betr = MAX(betr,bet-alf) jn = jn + n END DO lq = -1 IF (betr <= df*half) RETURN lq = 1 IF (l <= 0) RETURN betr = tmparray(ma+l) IF (betr <= zero) THEN IF (t < one .OR. betr == zero) RETURN lq = 2 END IF alfr = tmparray(l) ! Iteration loop ds_iteration: DO IF (lq >= 1) THEN q = one - betr/df r = q + SQRT(q*q + (alfr+alfr)/df) IF (betr >= zero) r = - (alfr+alfr)/ (df*r) r = MIN(1.95_prec,MAX(zero,r)) ELSE IF (ABS(betr-betl)+ABS(alfr-alfl)< -1.0E-4_prec*df) RETURN IF (betr-betl == zero) THEN iterm = -6 RETURN END IF r = (alfr-alfl)/ (betr-betl) END IF IF (ABS(t-r) < 1.0E-4_prec) RETURN t = r tmparray(l) = - one w = t*betr - alfr l = -1 DO j = 1,k alf = tmparray(j) IF (alf < zero) EXIT bet = tmparray(ma+j) r = t*bet - alf IF (r > w) THEN w = r l = j END IF END DO IF (l < 0) RETURN bet = tmparray(ma+l) IF (bet == zero) RETURN ! New interval selection alf = tmparray(l) IF (bet < zero) THEN IF (lq == 2) THEN alfr = alf betr = bet ELSE alfl = alf betl = bet lq = 0 END IF ELSE IF (lq == 2) THEN alfl = alfr betl = betr lq = 0 END IF alfr = alf betr = bet END IF END DO ds_iteration END SUBROUTINE destep SUBROUTINE nulstep(n,ma,mal,x,af,ag,ay,ibun,d,f,df,t,eta,iterm) ! Stepsize selection ! USE param, ONLY : zero,half,one,large ! given in host IMPLICIT NONE ! Scalar Arguments REAL(KIND=prec), INTENT(INOUT) :: & t ! Initial stepsize REAL(KIND=prec), INTENT(IN) :: & df, & ! Directional derivative. eta, & ! Distance measure parameter. f ! Value of the objective function. INTEGER, INTENT(IN) :: & n, & ! Number of variables ma, & ! Maximum size of the bundle. mal, & ! Current size of the bundle. ibun ! Index for the circular arrays in bundle updating. INTEGER, INTENT(OUT) :: & iterm ! Cause of termination: ! 0 - Everything is ok. ! -6 - Error. ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & x, & ! Vector of variables (n array). d, & ! Direction vector (n array). ay, & ! Matrix whose columns are bundle points (stored in one-dimensional n*ma -array). ag, & ! Matrix whose columns are bundle subgradients (stored in one-dimensional n*ma -array). af ! Vector of values of bundle functions (stored in one-dimensional 4*ma -array). ! Local Arrays REAL(KIND=prec), DIMENSION(2*ma) :: & tmparray ! Auxiliary array. ! Local Scalars REAL(KIND=prec) :: alf,alfl,alfr,bet,betl,betr,dx,q,r,w INTEGER :: i,j,jn,k,l,ib ! Intrinsic Functions INTRINSIC ABS,REAL,MAX,MIN,SQRT iterm = 0 w = df*t ! Initial choice of possibly active parabolas k = 0 l = -1 jn = (ibun-1)*n betr = - large DO j = 1,mal - 1 ib = ibun - 1 + j IF (ib > mal) ib = ib - mal IF (jn >= mal*n) jn = jn - mal*n bet = zero r = zero alfl = af(ib) - f DO i = 1,n dx = x(i) - ay(jn+i) r = r + dx*dx q = ag(jn+i) alfl = alfl + dx*q bet = bet + d(i)*q END DO alf = MAX(ABS(alfl),eta*r) betr = MAX(betr,bet-alf) IF (alf < bet-df) THEN k = k + 1 r = t*bet - alf tmparray(k) = alf tmparray(ma+k) = bet IF (r > w) THEN w = r l = k END IF END IF jn = jn + n END DO IF (l <= 0) RETURN betr = tmparray(ma+l) alfr = tmparray(l) alf = alfr bet = betr alfl = zero betl = df ! Iteration loop ns_iteration: DO w = bet/df IF (ABS(betr-betl)+ABS(alfr-alfl)< -1.0E-4_prec*df) RETURN IF (betr-betl == zero) THEN iterm = -6 RETURN END IF r = (alfr-alfl)/ (betr-betl) IF (ABS(t-w) < ABS(t-r)) r = w q = t t = r IF (ABS(t-q) < 1.0E-3_prec) RETURN tmparray(l) = - one w = t*bet - alf l = -1 DO j=1,k alf = tmparray(j) IF (alf < zero) EXIT bet = tmparray(ma+j) r = t*bet - alf IF (r > w) THEN w = r l = j END IF END DO IF (l <= 0) RETURN bet = tmparray(ma+l) q = bet - t*df IF (Q == zero) RETURN ! New interval selection alf = tmparray(l) IF (q < zero) THEN alfl = alf betl = bet ELSE alfr = alf betr = bet END IF END DO ns_iteration END SUBROUTINE nulstep END SUBROUTINE tinit !************************************************************************ !* * !* * SUBROUTINE ldgbm_lls * * !* * !* Special line search for limited memory discrete gradient * !* bundle method. * !* * !************************************************************************ SUBROUTINE ldgbm_lls(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 lmbm_sub, 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. t, & ! Stepsize p ! Directional derivative. 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) :: & iters, & ! Null step indicator. ! 0 - Null step. ! 1 - Serious step. nfe,ngrad ! Number of function and discrete gradient evaluations. INTEGER, INTENT(OUT) :: & 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. 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 tl = zero tu = t fl = 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 .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,wk,one-half/(one-epst)) 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 = qint(tu,fl,fu,wk,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 END SUBROUTINE ldgbm_lls !************************************************************************ !* !* * SUBROUTINE ldgbm_dlbfgs * !* !* Matrix update and computation of the search direction d = -dm*g !* by the limited memory BFGS update. !* !************************************************************************ SUBROUTINE ldgbm_dlbfgs(n,mc,mcc,inew,ibfgs,iflag,d,g,gp,s,u,sm,um,rm, & umtum,cm,smtgp,umtgp,gamma,tmpn1,iscale) USE param, ONLY : zero,small,one,half ! half is used at sclpar USE lmbm_sub, ONLY : & vdot, & ! Dot product. vneg, & ! Copying of a vector with change of the sign. xdiffy, & ! Difference of two vectors. xsumy, & ! Sum of two vectors. scdiff, & ! Difference of the scaled vector and a vector. scsum, & ! Sum of a vector and the scaled vector. vxdiag, & ! Multiplication of a vector and a diagonal matrix. symax, & ! Multiplication of a dense symmetric matrix by a vector. cwmaxv, & ! Multiplication of a vector by a dense rectangular matrix. rwaxv2, & ! Multiplication of two rowwise stored dense rectangular ! matrices A and B by vectors x and y. trlieq, & ! Solving x from linear equation L*x=y or trans(L)*x=y. copy2 ! Copying of two vectors. IMPLICIT NONE ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & g, & ! Current subgradient of the objective function. gp ! Previous subgradient of the objective function. REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: & sm, & ! Matrix whose columns are stored corrections. um, & ! Matrix whose columns are stored subgradient differences. rm, & ! Upper triangular matrix. cm, & ! Diagonal matrix. umtum, & ! Matrix umtum = trans(um) * um. smtgp, & ! Vector smtgp = trans(sm)*gp. umtgp, & ! vector umtgp = trans(um)*gp. s, & ! Difference of current and previous variables. u, & ! Difference of current and previous subgradients. tmpn1 ! Auxiliary array. On input: previous aggregate subgradient. REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: & d ! Direction vector. ! Scalar Arguments REAL(KIND=prec), INTENT(INOUT) :: & gamma ! Scaling parameter. INTEGER, INTENT(IN) :: & n, & ! Number of variables mc, & ! Declared number of stored corrections. 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. ! 4 - Preliminary scaling with STU/UTU. ! 5 - Preliminary scaling with STS/STU. ! 6 - No scaling. INTEGER, INTENT(INOUT) :: & mcc, & ! Current number of stored corrections. inew, & ! Index for circular arrays. iflag ! Index for adaptive version: ! 0 - Maximum number of stored corrections ! has not been changed at this iteration. ! 1 - Maximum number of stored corrections ! has been changed at this iteration. INTEGER, INTENT(OUT) :: & ibfgs ! Index of the type of BFGS update: ! 1 - BFGS update: the corrections are stored. ! 2 - BFGS update: the corrections are not stored. ! 3 - BFGS update is skipped. ! Local Arrays REAL(KIND=prec), DIMENSION(mcc+1) :: & tmpmc1,tmpmc2,tmpmc3,tmpmc4 ! Local Scalars REAL(KIND=prec) :: & stu, & ! stu = trans(s)*u. sts ! sts = trans(s)*s. INTEGER :: i,j,k, & mcnew, & ! Current size of vectors. iold, & ! Index of the oldest corrections. iflag2, & ! Index for adaptive version. ierr ! Error indicator ! Intrinsic Functions INTRINSIC SQRT,MIN,MAX ierr = 0 ibfgs = 0 iflag2 = 0 stu = vdot(n,s,u) sts = vdot(n,s,s) ! Positive definiteness IF (stu > zero) THEN IF (-vdot(n,d,u)-vdot(n,tmpn1,s) < -small) THEN ! Update matrices ibfgs = 1 ! Initialization of indices. CALL indic1(mc,mcc,mcnew,inew,iold,iflag,iflag2,ibfgs) ! Update sm and um CALL copy2(n,s,sm((inew-1)*n+1:),u,um((inew-1)*n+1:)) ! Computation of trans(sm)*g and trans(um)*g IF (inew >= mcnew) THEN CALL rwaxv2(n,mcnew,sm((inew-mcnew)*n+1:),& um((inew-mcnew)*n+1:),g,g,tmpmc1(iold:),tmpmc2(iold:)) ELSE CALL rwaxv2(n,inew,sm,um,g,g,tmpmc1,tmpmc2) CALL rwaxv2(n,mcnew-inew,sm((iold-1)*n+1:),& um((iold-1)*n+1:),g,g,tmpmc1(iold:),tmpmc2(iold:)) END IF ! Computation of trans(sm)*u and trans(um)*u IF (inew >= mcnew) THEN DO i=iold,inew-1 tmpmc3(i) = tmpmc1(i) - smtgp(i) smtgp(i) = tmpmc1(i) tmpmc4(i) = tmpmc2(i) - umtgp(i) umtgp(i) = tmpmc2(i) END DO ELSE DO i=1,inew-1 tmpmc3(i) = tmpmc1(i) - smtgp(i) smtgp(i) = tmpmc1(i) tmpmc4(i) = tmpmc2(i) - umtgp(i) umtgp(i) = tmpmc2(i) END DO DO i=iold,mcnew+1 tmpmc3(i) = tmpmc1(i) - smtgp(i) smtgp(i) = tmpmc1(i) tmpmc4(i) = tmpmc2(i) - umtgp(i) umtgp(i) = tmpmc2(i) END DO END IF tmpmc3(inew) = tmpmc1(inew) - vdot(n,s,gp) smtgp(inew) = tmpmc1(inew) tmpmc4(inew) = tmpmc2(inew) - vdot(n,u,gp) umtgp(inew) = tmpmc2(inew) ! Update rm and umtum IF (mcc >= mc .AND. iflag2 /= 1) THEN DO i=1,mcnew-1 j=(i-1)*i/2+1 k=i*(i+1)/2+2 CALL copy2(i,rm(k:),rm(j:),umtum(k:),umtum(j:)) END DO END IF IF (inew >= mcnew) THEN CALL copy2(mcnew,tmpmc3(iold:),rm((mcnew-1)*mcnew/2+1:),& tmpmc4(iold:),umtum((mcnew-1)*mcnew/2+1:)) ELSE CALL copy2(mcnew-inew,tmpmc3(iold:),rm((mcnew-1)*mcnew/2+1:)& ,tmpmc4(iold:),umtum((mcnew-1)*mcnew/2+1:)) CALL copy2(inew,tmpmc3,rm((mcnew-1)*mcnew/2+mcnew-inew+1:),& tmpmc4,umtum((mcnew-1)*mcnew/2+mcnew-inew+1:)) END IF ! Update CM cm(inew) = stu ! Computation of gamma gamma = ldgbm_sclpar(mcc,iscale,sts,stu,tmpmc4(inew)) inew = inew + 1 IF (inew > mc + 1) inew = 1 IF (iflag == 0 .AND. mcc < mc + 1) mcc = mcc + 1 ELSE ! BFGS update, corrections are not saved. ibfgs = 2 ! Initialization of indices. CALL indic1(mc,mcc,mcnew,inew,iold,iflag,iflag2,ibfgs) ! Update sm and um CALL copy2(n,s,sm((inew-1)*n+1:),u,um((inew-1)*n+1:)) ! Computation of trans(sm)*g and trans(um)*g CALL rwaxv2(n,mcnew,sm,um,g,g,tmpmc1,tmpmc2) ! Computation of trans(sm)*u and trans(um)*u IF (iold /= 1) THEN DO i=1,inew-1 tmpmc3(i) = tmpmc1(i) - smtgp(i) smtgp(i) = tmpmc1(i) tmpmc4(i) = tmpmc2(i) - umtgp(i) umtgp(i) = tmpmc2(i) END DO DO i=iold,mcnew tmpmc3(i) = tmpmc1(i) - smtgp(i) smtgp(i) = tmpmc1(i) tmpmc4(i) = tmpmc2(i) - umtgp(i) umtgp(i) = tmpmc2(i) END DO ELSE DO i=1,mcnew-1 tmpmc3(i) = tmpmc1(i) - smtgp(i) smtgp(i) = tmpmc1(i) tmpmc4(i) = tmpmc2(i) - umtgp(i) umtgp(i) = tmpmc2(i) END DO END IF tmpmc3(inew) = tmpmc1(inew) - vdot(n,s,gp) smtgp(inew) = tmpmc1(inew) tmpmc4(inew) = tmpmc2(inew) - vdot(n,u,gp) umtgp(inew) = tmpmc2(inew) ! Update rm and umtum IF (iold /= 1) THEN CALL copy2(mcnew-inew,tmpmc3(iold:),& rm((mcnew-1)*mcnew/2+1:),tmpmc4(iold:),& umtum((mcnew-1)*mcnew/2+1:)) CALL copy2(inew,tmpmc3,rm((mcnew-1)*mcnew/2+mcnew-inew+1:),tmpmc4,& umtum((mcnew-1)*mcnew/2+mcnew-inew+1:)) ELSE CALL copy2(mcnew,tmpmc3,rm((mcnew-1)*mcnew/2+1:)& ,tmpmc4,umtum((mcnew-1)*mcnew/2+1:)) END IF ! Update cm cm(inew) = stu ! Computation of gamma gamma = ldgbm_sclpar(mcc,iscale,sts,stu,tmpmc4(inew)) END IF ELSE ! BFGS update is skipped ibfgs = 3 IF (mcc == 0) THEN iflag = 0 CALL vneg(n,g,d) RETURN END IF ! Initialization of indices. CALL indic1(mc,mcc,mcnew,inew,iold,iflag,iflag2,ibfgs) ! Computation of gamma IF (iscale >= 4) gamma = one ! Computation of trans(sm)*g and trans(um)*g and the two intermediate values IF (iold <= 2) THEN CALL rwaxv2(n,mcnew,sm((iold-1)*n+1:),um((iold-1)*n+1:),g,g,& smtgp(iold:),umtgp(iold:)) ELSE CALL rwaxv2(n,inew-1,sm,um,g,g,smtgp,umtgp) CALL rwaxv2(n,mcnew-inew+1,sm((iold-1)*n+1:),& um((iold-1)*n+1:),g,g,smtgp(iold:),umtgp(iold:)) END IF END IF ! Computation of two intermediate values tmpmc1 and tmpmc2 IF (iold == 1 .OR. ibfgs == 2) THEN CALL trlieq(mcnew,mcnew,iold,rm,tmpmc1,smtgp,1,ierr) CALL symax(mcnew,mcnew,iold,umtum,tmpmc1,tmpmc3) CALL vxdiag(mcnew,cm,tmpmc1,tmpmc2) CALL scsum(mcnew,gamma,tmpmc3,tmpmc2,tmpmc2) CALL scsum(mcnew,-gamma,umtgp,tmpmc2,tmpmc3) CALL trlieq(mcnew,mcnew,iold,rm,tmpmc2,tmpmc3,0,ierr) ELSE IF (iflag == 0) THEN CALL trlieq(mcnew,mc+1,iold,rm,tmpmc1,smtgp,1,ierr) CALL symax(mcnew,mc+1,iold,umtum,tmpmc1,tmpmc3) CALL vxdiag(mc+1,cm,tmpmc1,tmpmc2) CALL scsum(mc+1,gamma,tmpmc3,tmpmc2,tmpmc2) CALL scsum(mc+1,-gamma,umtgp,tmpmc2,tmpmc3) CALL trlieq(mcnew,mc+1,iold,rm,tmpmc2,tmpmc3,0,ierr) ELSE CALL trlieq(mcnew,mc,iold,rm,tmpmc1,smtgp,1,ierr) CALL symax(mcnew,mc,iold,umtum,tmpmc1,tmpmc3) CALL vxdiag(mc,cm,tmpmc1,tmpmc2) CALL scsum(mc,gamma,tmpmc3,tmpmc2,tmpmc2) CALL scsum(mc,-gamma,umtgp,tmpmc2,tmpmc3) CALL trlieq(mcnew,mc,iold,rm,tmpmc2,tmpmc3,0,ierr) END IF ! Computation of the search direction d IF (iold == 1 .OR. ibfgs == 2) THEN CALL cwmaxv(n,mcnew,um,tmpmc1,d) CALL xdiffy(n,d,g,d) CALL cwmaxv(n,mcnew,sm,tmpmc2,tmpn1) CALL scdiff(n,gamma,d,tmpn1,d) ELSE CALL cwmaxv(n,inew-1,um,tmpmc1,d) CALL cwmaxv(n,mcnew-inew+1,um((iold-1)*n+1:),tmpmc1(iold:),tmpn1) CALL xsumy(n,d,tmpn1,d) CALL xdiffy(n,d,g,d) CALL cwmaxv(n,inew-1,sm,tmpmc2,tmpn1) CALL scdiff(n,gamma,d,tmpn1,d) CALL cwmaxv(n,mcnew-inew+1,sm((iold-1)*n+1:),tmpmc2(iold:),tmpn1) CALL xdiffy(n,d,tmpn1,d) END IF CONTAINS FUNCTION ldgbm_sclpar(mcc,iscale,sts,stu,utu) RESULT(spar) ! Calculation of the scaling parameter. ! USE param, ONLY : small,one,half ! given in host 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) :: & mcc, & ! Current number of stored corrections. 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. ! 4 - Preliminary scaling with STU/UTU. ! 5 - Preliminary 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,4) IF (utu < SQRT(small)) THEN spar = one RETURN ELSE spar = stu/utu END IF ! Scaling parameter = STS/STU CASE(1,3,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 IF (MCC == 0) THEN IF (spar < 0.01_prec) spar=0.01_prec IF (spar > 100.0_prec) spar=100.0_prec ELSE SELECT CASE(iscale) ! Interval scaling CASE(2) IF (spar < 0.6_prec .OR. spar > 6.0_prec) spar = one CASE(3) IF (spar < half .OR. spar > 5.0_prec) spar = one ! Preliminary scaling CASE(4,5) spar = one ! Scaling at every iteration CASE DEFAULT CONTINUE END SELECT END IF IF (spar < 1.0E+03_prec*small) spar = 1.0E+03_prec*small END FUNCTION ldgbm_sclpar END SUBROUTINE ldgbm_dlbfgs !************************************************************************ !* !* * SUBROUTINE dlsr1 * !* !* Matrix update and computation of the search direction d = -dm*ga !* by the limited memory SR1 update. !* !************************************************************************ SUBROUTINE dlsr1(n,mc,mcc,inew,isr1,iflag,d,gp,ga,s,u,sm,um,rm,& umtum,cm,smtgp,umtgp,gamma,tmpmc1,tmpmc2,tmpn1,nnk,iprint) USE param, ONLY : zero,small,one USE lmbm_sub, ONLY : & vdot, & ! Dot product. vneg, & ! Copying of a vector with change of the sign. scalex, & ! Scaling a vector. xdiffy, & ! Difference of two vectors. scdiff, & ! Difference of the scaled vector and a vector. xsumy, & ! Sum of two vectors. cwmaxv, & ! Multiplication of a vector by a dense rectangular matrix. rwaxv2, & ! Multiplication of two rowwise stored dense rectangular ! matrices A and B by vectors x and y. calq, & ! Solving x from linear equation A*x=y. copy, & ! Copying of a vector. copy2 ! Copying of two vectors. IMPLICIT NONE ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & ga, & ! Current aggregate subgradient of the objective function. gp ! Basic subgradient of the objective function. REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: & sm, & ! Matrix whose columns are stored corrections. um, & ! Matrix whose columns are stored subgradient differences. rm, & ! Upper triangular matrix. cm, & ! Diagonal matrix. umtum, & ! Matrix umtum = trans(um) * um. smtgp, & ! Vector smtgp = trans(sm)*gp. umtgp, & ! vector umtgp = trans(um)*gp. s, & ! Difference of current and previous variables. u, & ! Difference of current and previous subgradients. tmpn1 ! Auxiliary array. On input: previous aggregate subgradient. REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: & tmpmc1, & ! Auxiliary array. On output: trans(sm)*ga. tmpmc2, & ! Auxiliary array. On output: trans(um)*ga. d ! Direction vector. ! Scalar Arguments REAL(KIND=prec), INTENT(INOUT) :: & gamma ! Scaling parameter. INTEGER, INTENT(IN) :: & n, & ! Number of variables mc, & ! Declared number of stored corrections. nnk, & ! Consecutive null steps counter. iprint ! Printout specification. INTEGER, INTENT(INOUT) :: & mcc, & ! Current number of stored corrections. inew, & ! Index for circular arrays. iflag ! Index for adaptive version: ! 0 - Maximum number of stored corrections ! has not been changed at this iteration. ! 1 - Maximum number of stored corrections ! has been changed at this iteration. INTEGER, INTENT(OUT) :: & isr1 ! Index of the type of L-SR1 update: ! 1 - SR1 update: the corrections are stored. ! 3 - SR1 update is skipped. ! Local Arrays REAL(KIND=prec), DIMENSION(n) :: tmpn2 REAL(KIND=prec), DIMENSION((mcc+1)*(mcc+2)/2) :: tmpmat REAL(KIND=prec), DIMENSION(mcc+1) :: tmpmc3,tmpmc4,tmpmc5,tmpmc6 ! Local Scalars REAL(KIND=prec) :: & stu, & ! stu = trans(s)*u. a, & ! a = trans(ga) dm_(k-1) ga. b ! b = trans(ga) dm_k ga. INTEGER :: i,j,k, & mcnew, & ! Current size of vectors. iold, & ! Index of the oldest corrections. iflag2 ! Index for adaptive version. iflag2 = 0 isr1 = 0 ! Initialization of indices CALL indic1(mc,mcc,mcnew,inew,iold,iflag,iflag2,3) ! Computation of GAMMA gamma = one ! Computation of trans(sm)*ga and trans(um)*ga IF (iold <= 2) THEN CALL rwaxv2(n,mcnew,sm((iold-1)*n+1:),um((iold-1)*n+1:),ga,ga,& tmpmc1(iold:),tmpmc2(iold:)) ELSE CALL rwaxv2(n,inew-1,sm,um,ga,ga,tmpmc1,tmpmc2) CALL rwaxv2(n,mcnew-inew+1,sm((iold-1)*n+1:),um((iold-1)*n+1:),& ga,ga,tmpmc1(iold:),tmpmc2(iold:)) END IF ! Positive definiteness IF (-vdot(n,d,u) - vdot(n,tmpn1,s) >= -small) THEN ! SR1 update is skipped isr1 = 3 IF (mcc == 0) THEN iflag = 0 CALL vneg(n,ga,d) RETURN END IF ELSE stu = vdot(n,s,u) tmpmc1(inew) = vdot(n,s,ga) tmpmc2(inew) = vdot(n,u,ga) ! Convergence conditions IF ((nnk == 1 .OR. mcc < mc) .OR. & (iflag == 1 .AND. (inew == 1 .OR. inew == mc))) THEN ! SR1 Update isr1 = 1 ! Initialization of indices CALL indic1(mc,mcc,mcnew,inew,iold,iflag,iflag2,1) if (iflag2 == 1 .and. iold == 2) then tmpmc1(inew) = tmpmc1(1) tmpmc2(inew) = tmpmc2(1) end if ! Update sm and um CALL copy2(n,s,sm((inew-1)*n+1:),u,um((inew-1)*n+1:)) ! Update trans(sm)*gp and trans(um)*gp smtgp(inew) = vdot(n,s,gp) umtgp(inew) = vdot(n,u,gp) ! Computation of trans(sm)*u and trans(um)*u IF (iold <= 2) THEN CALL rwaxv2(n,mcnew-1,sm((iold-1)*n+1:),um((iold-1)*n+1:),& u,u,tmpmc3(iold:),tmpmc4(iold:)) ELSE CALL rwaxv2(n,inew-1,sm,um,u,u,tmpmc3,tmpmc4) CALL rwaxv2(n,mcnew-inew,sm((iold-1)*n+1:),um((iold-1)*n+1:),& u,u,tmpmc3(iold:),tmpmc4(iold:)) END IF tmpmc3(inew) = stu tmpmc4(inew) = vdot(n,u,u) ! Update rm and umtum IF (mcc >= mc .AND. iflag2 /= 1) THEN DO i=1,mcnew-1 j=(i-1)*i/2+1 k=i*(i+1)/2+2 CALL copy2(i,rm(k:),rm(j:),umtum(k:),umtum(j:)) END DO END IF IF (inew >= mcnew) THEN CALL copy2(mcnew,tmpmc3(iold:),rm((mcnew-1)*mcnew/2+1:),& tmpmc4(iold:),umtum((mcnew-1)*mcnew/2+1:)) ELSE CALL copy2(mcnew-inew,tmpmc3(iold:),& rm((mcnew-1)*mcnew/2+1:),tmpmc4(iold:),& umtum((mcnew-1)*mcnew/2+1:)) CALL COPY2(inew,tmpmc3,rm((mcnew-1)*mcnew/2+mcnew-inew+1:),tmpmc4,& umtum((mcnew-1)*mcnew/2+mcnew-inew+1:)) END IF ! Update cm cm(inew) = stu inew = inew + 1 IF (inew > mc + 1) inew = 1 IF (iflag == 0 .AND. mcc < mc + 1) mcc = mcc + 1 ELSE ! Calculation of matrix (umtum-rm-trans(rm)+cm) from previous iteration DO i=1,mcnew*(mcnew+1)/2 tmpmat(i)= gamma * umtum(i) - rm(i) END DO ! Computation of tmpmat*tmpmc4 = gamma*trans(um)*ga-trans(sm)*ga IF (iold == 1) THEN CALL scdiff(mcnew,gamma,tmpmc2,tmpmc1,tmpmc5) CALL calq(mcnew,mcnew,iold,tmpmat,tmpmc4,tmpmc5,0) ELSE IF (iflag == 0) THEN CALL scdiff(mc+1,gamma,tmpmc2,tmpmc1,tmpmc5) CALL calq(mcnew,mc+1,iold,tmpmat,tmpmc4,tmpmc5,0) ELSE CALL scdiff(mc,gamma,tmpmc2,tmpmc1,tmpmc5) CALL calq(mcnew,mc,iold,tmpmat,tmpmc4,tmpmc5,0) END IF ! Computation of a = -trans(ga)*dm_(k-1)*ga IF (iold <= 2) THEN CALL scalex(mcnew,gamma,tmpmc4(iold:),tmpmc3(iold:)) CALL cwmaxv(n,mcnew,sm((iold-1)*n+1:),tmpmc4(iold:),tmpn1) CALL scdiff(n,-gamma,ga,tmpn1,tmpn2) CALL cwmaxv(n,mcnew,um((iold-1)*n+1:),tmpmc3(iold:),tmpn1) CALL xsumy(n,tmpn2,tmpn1,tmpn2) ELSE CALL scalex(mcc,gamma,tmpmc4,tmpmc3) CALL cwmaxv(n,inew-1,sm,tmpmc4,tmpn1) CALL scdiff(n,-gamma,ga,tmpn1,tmpn2) CALL cwmaxv(n,mcnew-inew+1,sm((iold-1)*n+1:),tmpmc4(iold:),& tmpn1) CALL xdiffy(n,tmpn2,tmpn1,tmpn2) CALL cwmaxv(n,inew-1,um,tmpmc3,tmpn1) CALL xsumy(n,tmpn2,tmpn1,tmpn2) CALL cwmaxv(n,mcnew-inew+1,um((iold-1)*n+1:),tmpmc3(iold:),& tmpn1) CALL xsumy(n,tmpn2,tmpn1,tmpn2) END IF a = vdot(n,ga,tmpn2) IF (iflag == 0) THEN mcnew = mc iold = inew + 2 IF (iold > mc+1) iold = iold - mc - 1 ELSE mcnew = mc - 1 iold = inew + 2 IF (iold > mc) iold = iold - mc END IF ! Calculation of the new canditate for search direction ! Updates are not necessarily saved ! Update sm and um CALL copy2(n,s,sm((inew-1)*n+1:),u,um((inew-1)*n+1:)) ! Computation of trans(sm)*u and trans(um)*u IF (iold == 1 .OR. iold == 2) THEN CALL rwaxv2(n,mcnew-1,sm((iold-1)*n+1:),um((iold-1)*n+1:),u,u,& tmpmc3(iold:),tmpmc4(iold:)) ELSE CALL rwaxv2(n,inew-1,sm,um,u,u,tmpmc3,tmpmc4) CALL rwaxv2(n,mcnew-inew,sm((iold-1)*n+1:),um((iold-1)*n+1:),u,u,& tmpmc3(iold:),tmpmc4(iold:)) END IF tmpmc3(inew) = stu tmpmc4(inew) = vdot(n,u,u) ! Calculation of matrix (umtum-rm-trans(rm)+cm) without updating ! matrices rm, umtum and cm DO i=1,mcnew*(mcnew+1)/2 tmpmat(i)= gamma * umtum(i) - rm(i) END DO DO i=1,mcnew-1 j=(i-1)*i/2+1 k=i*(i+1)/2+2 CALL copy(i,tmpmat(k:),tmpmat(j:)) END DO CALL scdiff(mcnew+1,gamma,tmpmc4,tmpmc3,tmpmc5) IF (inew >= mcnew) THEN CALL copy(mcnew,tmpmc5(iold:),tmpmat((mcnew-1)*mcnew/2+1:)) ELSE CALL copy(mcnew-inew,tmpmc5(iold:),tmpmat((mcnew-1)*mcnew/2+1:)) CALL copy(inew,tmpmc5,tmpmat((mcnew-1)*mcnew/2+mcnew-inew+1:)) END IF IF (iflag == 0) THEN CALL scdiff(mc+1,gamma,tmpmc2,tmpmc1,tmpmc5) CALL calq(mcnew,mc+1,iold,tmpmat,tmpmc5,tmpmc5,iprint) ELSE CALL scdiff(mc,gamma,tmpmc2,tmpmc1,tmpmc5) CALL calq(mcnew,mc,iold,tmpmat,tmpmc5,tmpmc5,iprint) END IF ! Calculation of the new canditate for search direction d = -dm_k*ga ! and computation of b = -trans(ga)*dm_k*ga IF (iold <= 2) THEN CALL scalex(mcnew,gamma,tmpmc5(iold:),tmpmc6(iold:)) CALL cwmaxv(n,mcnew,sm((iold-1)*n+1:),tmpmc5(iold:),tmpn1) CALL scdiff(n,-gamma,ga,tmpn1,d) CALL cwmaxv(n,mcnew,um((iold-1)*n+1:),tmpmc6(iold:),tmpn1) CALL xsumy(n,d,tmpn1,d) ELSE CALL scalex(mcnew+1,gamma,tmpmc5,tmpmc6) CALL cwmaxv(n,inew,sm,tmpmc5,tmpn1) CALL scdiff(n,-gamma,ga,tmpn1,d) CALL cwmaxv(n,mcnew-inew,sm((iold-1)*n+1:),tmpmc5(iold:),& tmpn1) CALL xdiffy(n,d,tmpn1,d) CALL cwmaxv(n,inew,um,tmpmc6,tmpn1) CALL xsumy(n,d,tmpn1,d) CALL cwmaxv(n,mcnew-inew,um((iold-1)*n+1:),tmpmc6(iold:),& tmpn1) CALL xsumy(n,d,tmpn1,d) END IF b = vdot(n,ga,d) ! Checking the convergence conditions IF (b - a < zero) THEN isr1 = 3 CALL copy(n,tmpn2,d) ELSE isr1 = 1 ! Update trans(sm)*gp and trans(um)*gp smtgp(inew) = vdot(n,s,gp) umtgp(inew) = vdot(n,u,gp) ! Update rm and umtum DO i=1,mcnew-1 j=(i-1)*i/2+1 k=i*(i+1)/2+2 CALL copy2(i,rm(k:),rm(j:),umtum(k:),umtum(j:)) END DO IF (inew >= mcnew) THEN CALL copy2(mcnew,tmpmc3(iold:),rm((mcnew-1)*mcnew/2+1:),tmpmc4(iold:),& umtum((mcnew-1)*mcnew/2+1:)) ELSE CALL copy2(mcnew-inew,tmpmc3(iold:),rm((mcnew-1)*mcnew/2+1:),tmpmc4(iold:),& umtum((mcnew-1)*mcnew/2+1:)) CALL copy2(inew,tmpmc3,rm((mcnew-1)*mcnew/2+mcnew-inew+1:),tmpmc4,& umtum((mcnew-1)*mcnew/2+mcnew-inew+1:)) END IF ! Update cm cm(inew) = stu inew = inew + 1 IF (inew > mc + 1) inew = 1 IF (iflag == 0 .AND. mcc < mc + 1) mcc = mcc + 1 END IF RETURN END IF END IF DO i=1,mcnew*(mcnew+1)/2 tmpmat(i)= gamma * umtum(i) - rm(I) END DO ! Computation of tmpmat*tmpmc4 = gamma*trans(um)*ga-trans(sm)*ga IF (iold == 1) THEN CALL scdiff(mcnew,gamma,tmpmc2,tmpmc1,tmpmc4) CALL calq(mcnew,mcnew,iold,tmpmat,tmpmc4,tmpmc4,iprint) ELSE IF (iflag == 0) THEN CALL scdiff(mc+1,gamma,tmpmc2,tmpmc1,tmpmc4) CALL calq(mcnew,mc+1,iold,tmpmat,tmpmc4,tmpmc4,iprint) ELSE CALL scdiff(mc,gamma,tmpmc2,tmpmc1,tmpmc4) CALL calq(mcnew,mc,iold,tmpmat,tmpmc4,tmpmc4,iprint) END IF ! Computation of the search direction d IF (iold <= 2) THEN CALL scalex(mcnew,gamma,tmpmc4(iold:),tmpmc3(iold:)) CALL cwmaxv(n,mcnew,sm((iold-1)*n+1:),tmpmc4(iold:),tmpn1) CALL scdiff(n,-gamma,ga,tmpn1,d) CALL cwmaxv(n,mcnew,um((iold-1)*n+1:),tmpmc3(iold:),tmpn1) CALL xsumy(n,d,tmpn1,d) ELSE CALL scalex(mcc,gamma,tmpmc4,tmpmc3) CALL cwmaxv(n,inew-1,sm,tmpmc4,tmpn1) CALL scdiff(n,-gamma,ga,tmpn1,d) CALL cwmaxv(n,mcnew-inew+1,sm((iold-1)*n+1:),tmpmc4(iold:),& tmpn1) CALL xdiffy(n,d,tmpn1,d) CALL cwmaxv(n,inew-1,um,tmpmc3,tmpn1) CALL xsumy(n,d,tmpn1,d) CALL cwmaxv(n,mcnew-inew+1,um((iold-1)*n+1:),tmpmc3(iold:),& tmpn1) CALL xsumy(n,d,tmpn1,d) END IF END SUBROUTINE dlsr1 !************************************************************************ !* * !* * SUBROUTINE indic1 * * !* * !* Initialization of indices. * !* * !************************************************************************ SUBROUTINE indic1(mc,mcc,mcnew,inew,iold,iflag,iflag2,itype) IMPLICIT NONE ! Scalar Arguments INTEGER, INTENT(IN) :: & mc, & ! Declared number of stored corrections. mcc, & ! Current number of stored corrections. itype ! Type of Initialization: ! 1 - corrections are stored, ! 2 - corrections are not stored, ! 3 - update is skipped. INTEGER, INTENT(INOUT) :: & inew, & ! Index for circular arrays. iflag, & ! Index for adaptive version: ! 0 - Maximum number of stored corrections ! has not been changed at this iteration. ! 1 - Maximum number of stored corrections ! has been changed at this iteration. iflag2 ! Index for adaptive version. ! 0 - iflag has not been changed. ! 1 - iflag has been changed. INTEGER, INTENT(OUT) :: & mcnew, & ! Current size of vectors. iold ! Index of the oldest corrections. IF (itype == 1) THEN IF (mcc < mc) THEN mcnew = mcc + 1 iold = 1 iflag = 0 ELSE IF (iflag == 0) THEN mcnew = mc iold = inew + 2 IF (iold > mc+1) iold = iold - mc - 1 ELSE IF (inew == 1) THEN inew = mc + 1 mcnew = mc iold = 2 iflag = 0 iflag2 = 1 ELSE IF (inew == mc) THEN mcnew = mc iold = 1 iflag = 0 iflag2 = 1 ELSE mcnew = mc - 1 iold = inew + 2 IF (iold > mc) iold = iold - mc END IF END IF END IF ELSE IF (itype == 2) THEN IF (mcc < mc) THEN mcnew = mcc + 1 iold = 1 iflag = 0 ELSE IF (iflag == 0) THEN mcnew = mc + 1 iold = inew + 1 IF (iold > mc + 1) iold = 1 ELSE mcnew = mc iold = inew + 1 IF (iold > mc) iold = 1 END IF END IF ELSE IF (mcc < mc) THEN mcnew = mcc iold = 1 iflag = 0 ELSE IF (iflag == 0) THEN mcnew = mc iold = inew + 1 IF (iold > mc + 1) iold = 1 ELSE mcnew = mc - 1 iold = inew + 1 IF (iold > mc) iold = 1 END IF END IF END IF END SUBROUTINE indic1 !************************************************************************ !* !* * SUBROUTINE agbfgs * !* !* Computation of aggregate values by the limited memory BFGS update. !* !************************************************************************ SUBROUTINE agbfgs(n,mc,mcc,inew,ibfgs,iflag,g,gp,ga,u,d,sm,um, & rm,cm,umtum,alfn,alfv,gamma,ic,rho) USE param, ONLY : zero,half,one USE lmbm_sub, ONLY : & symax, & ! Multiplication of a dense symmetric matrix by a vector. rwaxv2, & ! Multiplication of two rowwise stored dense rectangular ! matrices A and B by vectors x and y. trlieq, & ! Solving x from linear equation L*x=y or trans(L)*x=y. vdot ! Dot product. IMPLICIT NONE ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & d, & ! Direction vector. g, & ! Current (auxiliary) subgradient of the objective function. gp, & ! Previous subgradient of the objective function. u, & ! Difference of trial and aggregate gradients. sm, & ! Matrix whose columns are stored corrections. um, & ! Matrix whose columns are stored subgradient differences. rm, & ! Upper triangular matrix. umtum, & ! Matrix umtum = trans(um) * um. cm ! Diagonal matrix. REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: & ga ! 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. rho ! Correction parameter. INTEGER, INTENT(IN) :: & n, & ! Number of variables mc, & ! Declared number of stored corrections. mcc, & ! Current number of stored corrections. inew, & ! Index for circular arrays. ibfgs, & ! Index of the type of BFGS update. ic ! Correction indicator. INTEGER, INTENT(INOUT) :: & iflag ! Index for adaptive version: ! 0 - Maximum number of stored corrections ! has not been changed at this iteration. ! 1 - Maximum number of stored corrections ! has been changed at this iteration. ! Local arrays REAL(KIND=prec), DIMENSION(mcc+1) :: tmpmc1, tmpmc2 ! 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, & mcnew, & ! Current size of vectors. iold, & ! Index of the oldest corrections. ierr ! Error indicador. ! Intrinsic Functions INTRINSIC MAX,MIN,SIGN ierr = 0 IF (mcc < mc) THEN IF (ibfgs == 2) THEN mcnew = mcc + 1 ELSE mcnew = mcc END IF iold = 1 ELSE IF (iflag == 0) THEN IF (ibfgs == 2) THEN mcnew = mc + 1 ELSE mcnew = mc END IF iold = inew + 1 IF (iold > mc+1) iold = 1 ELSE IF (ibfgs == 2) THEN mcnew = mc ELSE mcnew = mc - 1 END IF iold = inew + 1 IF (iold > mc) iold = 1 END IF END IF ! Computation of trans(d) * u - alfn p = vdot(n,d,u) - alfn q = vdot(n,u,u) IF (ic == 1) THEN w = rho * q ELSE w = zero END IF ! Computation of the product trans(u)*dm*u IF (mcc > 0 .OR. ibfgs == 2) THEN IF (iold == 1 .OR. ibfgs == 2) THEN CALL rwaxv2(n,mcnew,sm,um,u,u,tmpmc1,tmpmc2) CALL trlieq(mcnew,mcnew,iold,rm,tmpmc1,tmpmc1,1,ierr) q = q - 2.0_prec*vdot(mcnew,tmpmc2,tmpmc1) q = gamma*q DO i=1,mcnew tmpmc2(i) = cm(i)*tmpmc1(i) END DO q = q + vdot(mcnew,tmpmc1,tmpmc2) CALL symax(mcnew,mcnew,iold,umtum,tmpmc1,tmpmc2) q = q + gamma*vdot(mcnew,tmpmc1,tmpmc2) ELSE CALL rwaxv2(n,inew-1,sm,um,u,u,tmpmc1,tmpmc2) CALL rwaxv2(n,mcc-inew,sm((iold-1)*n+1:),um((iold-1)*n+1:), & u,u,tmpmc1(iold:),tmpmc2(iold:)) CALL trlieq(mcnew,mcc,iold,rm,tmpmc1,tmpmc1,1,ierr) q = q - 2.0_prec*(vdot(mcc-inew,tmpmc2(iold:),tmpmc1(iold:)) + & vdot(inew-1,tmpmc2,tmpmc1)) q = gamma*q DO i=1,mcc tmpmc2(i) = cm(i)*tmpmc1(i) END DO q = q + vdot(mcc-inew,tmpmc1(iold:),tmpmc2(iold:)) + & vdot(inew-1,tmpmc1,tmpmc2) CALL symax(mcnew,mcc,iold,umtum,tmpmc1,tmpmc2) q = q + gamma*(vdot(mcc-inew,tmpmc1(iold:),tmpmc2(iold:)) + & vdot(inew-1,tmpmc1,tmpmc2)) END IF END IF q = q + w 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 ga(i)=lam*g(i) + p*gp(i) END DO alfv = lam*alfn END SUBROUTINE agbfgs !************************************************************************ !* !* * SUBROUTINE aggsr1 * !* !* Computation of aggregate values by the limited memory SR1 update. !* !************************************************************************ SUBROUTINE aggsr1(n,mc,mcc,inew,iflag,g,gp,ga,d,alfn,alfv, & umtum,rm,gamma,smtgp,umtgp,smtga,umtga,sm,um,icn,rho) USE param, ONLY : zero,one,small USE lmbm_sub, ONLY : & vdot, & ! Dot product. scalex, & ! Scaling a vector. xsumy, & ! Sum of two vectors. xdiffy, & ! Difference of two vectors. scsum, & ! Sum of a vector and the scaled vector. scdiff, & ! Difference of the scaled vector and a vector. rwaxv2, & ! Multiplication of two rowwise stored dense rectangular ! matrices A and B by vectors x and y. cwmaxv, & ! Multiplication of a vector by a dense rectangular matrix. lineq, & ! Solver for linear equation. calq ! Solving x from linear equation A*x=y. IMPLICIT NONE ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & d, & ! Direction vector. g, & ! Current (auxiliary) subgradient of the objective function. gp, & ! Previous subgradient of the objective function. sm, & ! Matrix whose columns are stored corrections. um, & ! Matrix whose columns are stored subgradient differences. rm, & ! Upper triangular matrix. umtum, & ! Matrix umtum = trans(um) * um. smtgp, & ! Vector smtgp = trans(sm)*gp. umtgp, & ! vector umtgp = trans(um)*gp. smtga, & ! vector smtga = trans(sm)*ga. umtga ! vector umtga = trans(um)*ga. REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: & ga ! 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. rho ! Correction parameter. INTEGER, INTENT(IN) :: & n, & ! Number of variables mc, & ! Declared number of stored corrections. mcc, & ! Current number of stored corrections. inew, & ! Index for circular arrays. icn ! Correction indicator. INTEGER, INTENT(INOUT) :: & iflag ! Index for adaptive version: ! 0 - Maximum number of stored corrections ! has not been changed at this iteration. ! 1 - Maximum number of stored corrections ! has been changed at this iteration. ! Local arrays REAL(KIND=prec), DIMENSION(n) :: tmpn2,tmpn3,tmpn4 REAL(KIND=prec), DIMENSION((mcc+1)*(mcc)/2) :: tmpmat REAL(KIND=prec), DIMENSION(mcc+1) :: tmpmc3, tmpmc4 ! Local Scalars REAL(KIND=prec) :: & pr, & ! pr = trans(gp-ga) dm (gp-ga), where dm ! presents the L-SR1- approximation of Hessian. rrp, & ! rrp = trans(gp-ga) dm ga - alfv. prqr, & ! prqr = trans(gp-ga) dm (g-ga). rrq, & ! rrq = trans(g-ga) dm ga - alfv + alfn. qr, & ! qr = trans(g-ga) dm (g-ga). pq, & ! pq = trans(g-gp) dm (g-gp). qqp, & ! qqp = trans(g-gp) dm g + alfn. lam1, & ! Multiplier used to calculate aggregate values. lam2, & ! Multiplier used to calculate aggregate values. w, & ! Correction. tmp1, & ! Auxiliary scalar. tmp2 ! Auxiliary scalar. INTEGER :: i, & mcnew, & ! Current size of vectors. iold, & ! Index of the oldest corrections. ierr ! Error indicador. ! Intrinsic Functions INTRINSIC MIN,MAX ierr = 0 IF (mcc < mc) THEN iold = 1 mcnew = mcc ELSE IF (iflag == 0) THEN mcnew = mc iold = inew + 1 IF (iold > mc+1) iold = 1 ELSE mcnew = mc - 1 iold = inew + 1 IF (iold > mc) iold = 1 END IF CALL xdiffy(n,gp,ga,tmpn2) ! Calculation of tmpn3 = trans(gp - ga)dm IF (mcc > 0) THEN DO i=1,mcnew*(mcnew+1)/2 tmpmat(i)= gamma * umtum(i) - rm(i) END DO IF (iold == 1) THEN CALL xdiffy(mcnew,umtgp,umtga,tmpmc4) CALL scdiff(mcnew,gamma,tmpmc4,smtgp,tmpmc4) CALL xsumy(mcnew,tmpmc4,smtga,tmpmc4) CALL calq(mcnew,mcnew,iold,tmpmat,tmpmc3,tmpmc4,0) CALL scalex(mcnew,gamma,tmpmc3,tmpmc4) CALL cwmaxv(n,mcnew,sm,tmpmc3,tmpn4) CALL scsum(n,gamma,tmpn2,tmpn4,tmpn3) CALL cwmaxv(n,mcnew,um,tmpmc4,tmpn4) CALL xdiffy(n,tmpn3,tmpn4,tmpn3) ELSE CALL xdiffy(mcc,umtgp,umtga,tmpmc4) CALL scdiff(mcc,gamma,tmpmc4,smtgp,tmpmc4) CALL xsumy(mcc,tmpmc4,smtga,tmpmc4) CALL calq(mcnew,mcc,iold,tmpmat,tmpmc3,tmpmc4,0) CALL scalex(mcc,gamma,tmpmc3,tmpmc4) CALL cwmaxv(n,inew-1,sm,tmpmc3,tmpn4) CALL scsum(n,gamma,tmpn2,tmpn4,tmpn3) CALL cwmaxv(n,mcnew-inew+1,sm((iold-1)*n+1:),tmpmc3(iold:)& ,tmpn4) CALL xsumy(n,tmpn3,tmpn4,tmpn3) CALL cwmaxv(n,inew-1,um,tmpmc4,tmpn4) CALL xdiffy(n,tmpn3,tmpn4,tmpn3) CALL cwmaxv(n,mcnew-inew+1,um((iold-1)*n+1:),tmpmc4(iold:)& ,tmpn4) CALL xdiffy(n,tmpn3,tmpn4,tmpn3) END IF IF (icn == 1) THEN CALL scsum(n,rho,tmpn2,tmpn3,tmpn3) END IF pr = vdot(n,tmpn3,tmpn2) rrp = vdot(n,tmpn3,ga) CALL xdiffy(n,g,ga,tmpn4) prqr = vdot(n,tmpn3,tmpn4) rrq = -vdot(n,tmpn4,d) ELSE pr = vdot(n,tmpn2,tmpn2) rrp = vdot(n,tmpn2,ga) CALL xdiffy(n,g,ga,tmpn4) prqr = vdot(n,tmpn2,tmpn4) rrq = -vdot(n,tmpn4,d) END IF ! calculation of qr = trans(g - ga) dm (g - ga) qr = vdot(n,tmpn4,tmpn4) IF (icn == 1) THEN w = rho*qr ELSE w = zero END IF IF (mcc > 0) THEN qr = gamma*qr IF (iold == 1) THEN CALL rwaxv2(n,mcnew,sm,um,tmpn4,tmpn4,tmpmc4,tmpmc3) CALL scsum(mcnew,-gamma,tmpmc3,tmpmc4,tmpmc4) CALL lineq(mcnew,mcnew,iold,tmpmat,tmpmc3,tmpmc4,ierr) qr = qr - vdot(mcnew,tmpmc4,tmpmc3) + w ELSE CALL rwaxv2(n,inew-1,sm,um,tmpn4,tmpn4,tmpmc4,tmpmc3) CALL rwaxv2(n,mcnew-inew+1,sm((iold-1)*n+1:),& um((iold-1)*n+1:),tmpn4,tmpn4,tmpmc4(iold:),tmpmc3(iold:)) CALL scsum(mcc,-gamma,tmpmc3,tmpmc4,tmpmc4) CALL lineq(mcnew,mcc,iold,tmpmat,tmpmc3,tmpmc4,ierr) qr = qr - vdot(mcc-inew,tmpmc4(iold:),tmpmc3(iold:)) -& vdot(inew-1,tmpmc4,tmpmc3) + w END IF END IF 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 w = pr - prqr*tmp2 IF (w /= zero) THEN lam1 = (tmp1*prqr - rrp)/w 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)) w = (lam2*qr + rrq+rrq)*lam2 ! w = (lam2*qr + 2.0_prec*rrq)*lam2 tmp1 = zero IF (alfv >= zero) tmp1 = one IF (pr > zero) tmp1 = MIN(one,MAX(zero,-rrp/pr)) ! tmp2 = (tmp1*pr + 2.0_prec*rrp)*tmp1 tmp2 = (tmp1*pr + rrp+rrp)*tmp1 IF (tmp2 < w) THEN w = tmp2 lam1 = tmp1 lam2 = zero END IF IF (qqp*(qqp - pq) < zero) THEN IF (qr + rrq + rrq - qqp*qqp/pq < W) 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 ga(i)=lam1*gp(i)+lam2*g(i)+tmp1*ga(i) END DO alfv = lam2*alfn + tmp1*alfv END SUBROUTINE aggsr1 !************************************************************************ !* * !* * 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 lmbm_sub, 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 ! Disgrete 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 ldgbm_wprint * !* !* Printout the warning and error messages. !* !************************************************************************ SUBROUTINE ldgbm_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. ! -1 - Two consecutive restarts. ! -2 - Number of restarts > maximum number of restarts. ! -3 - Failure in function or subgradient calculations ! (assigned by the user). ! -4 - Failure in attaining the demanded accuracy. ! -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 == 3) WRITE (6,FMT='(1X,''Error: '' & ''The size of the bundle (na) 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 calculations.'')') IF (iterm == -4) WRITE (6,FMT='(1X,''Abnormal exit: '' & ''Failure in attaining the demanded accuracy.'')') END IF END IF END SUBROUTINE ldgbm_wprint !************************************************************************ !* !* * SUBROUTINE ldgbm_rprint * !* !* Printout the (final) results. !* !************************************************************************ SUBROUTINE ldgbm_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. ! -1 - Two consecutive restarts. ! -2 - Number of restarts > maximum number of restarts. ! -3 - Failure in function or subgradient calculations ! (assigned by the user). ! -4 - Failure in attaining the demanded accuracy. ! -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 ldgbm_rprint END MODULE ldgbm_mod