MODULE dbundle_mod ! D-Bundle method USE r_precision, ONLY : prec ! Precision for reals. IMPLICIT NONE ! MODULE dbundle_mod includes the following subroutines (S) and functions (F). PUBLIC :: & init_db, & ! S Initialization for D-bundle subroutine. dbundle ! S D-bundle subroutine for nonsmooth large-scale ! optimization. Contains: ! S restar Initialization and reinitialization. ! S dobun Bundle construction. ! F sclpar Selection of scaling parameter. PRIVATE :: & 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. lls, & ! S Line search. Contains: ! F qint Quadratic interpolation. aggre1, & ! S Simplified subgradient aggregation. aggre2, & ! S Subgradient aggregation. wprint, & ! S Printout the error and warning messages. rprint ! S Printout the results. CONTAINS !*********************************************************************** !* * !* * SUBROUTINE init_db * * !* * !* Initialization for diagonal bundle subroutine for * !* large-scale unconstrained nonsmooth optimization. * !* * !*********************************************************************** SUBROUTINE init_db(iterm) USE param, ONLY : zero,half,one,small,large USE initializat, ONLY : & n, & ! Number of variables. na, & ! Maximum bundle dimension, na >= 2. mcu, & ! Maximum number of stored corrections, mcu >= 1. 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. 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 (na < 2) THEN iterm = -5 CALL wprint(iterm,iprint,3) 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 ldgbm_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 first termination criterion. IF (told <= zero) told = one ! 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 scaling. END SUBROUTINE init_db !*********************************************************************** !* * !* * SUBROUTINE dbundle * * !* * !* Diagonal bundle subroutine for nonsmooth optimization. * !* * !*********************************************************************** SUBROUTINE dbundle(f,nit,nfe,nge,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 na, & ! Maximum bundle dimension, na >= 2. mcu, & ! Maximum number of stored corrections, mcu >= 1. 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. 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. myg ! Computation of the subgradient 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) :: & nit,nfe,nge ! Number of iterations, and function and subgradient 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 or subgradient 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. g, & ! Subgradient of the objective function. gp, & ! Previous subgradient of the ohjective function. ga, & ! Aggregate subgradient. s, & ! Difference of current and previous variables. u, & ! Difference of current and previous subgradients. 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 subgradients. REAL(KIND=prec), DIMENSION(n*na) :: & ax,ag ! Matrix whose columns are bundle points and subgradients. REAL(KIND=prec), DIMENSION(na) :: & af ! Vector of bundle values. ! Local Scalars REAL(KIND=prec) :: & 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. stu, & ! stu = trans(s)*u. utu, & ! utu = trans(u)*u. sts, & ! sts = trans(s)*s. gamma_old, & ! Previous scaling parameter. gamma ! Scaling parameter. INTEGER :: i,j, & 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. nress, & ! Number of consecutive restarts in case of tmax < tmin. 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. tmin = 1.0E-12_prec, & ! Minimum stepsize. lengthd = 1.0E+20_prec ! Direction vector length. 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 D-Bundle:'')') ! Initialization nout = 0 nit = 0 nfe = 0 nge = 0 ntesf = 0 nres = 1 ncres = -1 nress = 0 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 ! Computation of the value and the subgradient of the objective ! function and the search direction for the first iteration CALL myf(n,x,f,iterm) CALL myg(n,x,g,iterm) nfe = nfe + 1 nge = nge + 1 IF (iterm /= 0) GO TO 900 CALL restar(n,mcc,inew,ibun,iters,gp,g,nnk,alfv,alfn,gamma,d,diag,mal,ncres) CALL dobun(n,na,mal,x,g,f,ax,ag,af,iters,ibun) ! Start of the iteration iteration: DO ! Computational time IF (time > 0.0E+00) THEN CALL getime(elapstim) IF (elapstim-strtim > time) THEN iterm = 6 EXIT iteration END IF END IF ! Computation of norms IF (iters > 0) THEN gnorm = vdot(n,g,g) dnorm = SQRT(vdot(n,d,d)) p = vdot(n,g,d) ELSE gnorm = vdot(n,ga,ga) dnorm = SQRT(vdot(n,d,d)) p = vdot(n,ga,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 wprint(iterm,iprint,-3) CALL restar(n,mcc,inew,ibun,iters,gp,g,nnk,alfv,alfn,gamma,d,diag,mal,ncres) IF (ncres > maxnrs) THEN nout = maxnrs iterm = -2 EXIT iteration END IF CALL dobun(n,na,mal,x,g,f,ax,ag,af,iters,ibun) CYCLE iteration END IF nout = -1 iterm = -1 EXIT iteration END IF ! Stopping criterion nit = nit + 1 pxnorm = xnorm xnorm = -p + 2.0_prec*alfv ! Tests for termination IF(xnorm <= tolg) THEN ! desired accuracy iterm = 1 EXIT iteration END IF IF (nfe >= mfe) THEN ! too many function calls nout = mfe iterm = 4 EXIT iteration END IF IF (nit >= mit) THEN ! too many iterations nout = mit iterm = 5 EXIT iteration END IF IF (f <= tolb) THEN ! too small f iterm = 7 EXIT iteration END IF IF (iters == 0) THEN IF (ABS(xnorm - pxnorm) <= small) THEN neps = neps + 1 IF (neps > maxeps) THEN iterm = 8 EXIT iteration 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,nit,nfe,nge,x,f,xnorm,half*gnorm+alfv,iterm,iprint) ! Preparation of line search fo = f IF (iters > 0) THEN CALL copy2(n,x,xo,g,gp) END IF IF (dnorm > zero) tmax = xmax/dnorm IF (tmax > tmin) THEN nress = 0 ELSE nress = nress + 1 IF (nress == 1) THEN CALL wprint(iterm,iprint,-5) CALL restar(n,mcc,inew,ibun,iters,gp,g,nnk,alfv,alfn,gamma,d,diag,mal,ncres) IF (ncres > maxnrs) THEN nout = maxnrs iterm = -2 EXIT iteration END IF CALL dobun(n,na,mal,x,g,f,ax,ag,af,iters,ibun) CYCLE iteration END IF iterm = -1 EXIT iteration END IF ! 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 iteration ! Line search with directional derivatives which allows null steps theta = one IF (dnorm > lengthd) THEN theta=lengthd/dnorm END IF CALL lls(n,x,g,d,xo,t,fo,f,p,alfn,tmin,dnorm,xnorm, & theta,epsl,epsr,eta,iters,nfe,nge,nnk,iterm) IF (iterm /= 0) EXIT iteration IF (tolf2 >= 0) THEN IF (ABS(fo-f) <= tolf2*small*MAX(ABS(f),ABS(fo),one) & .AND. iters == 1) THEN iterm = 3 EXIT iteration END IF END IF IF (ABS(fo-f) <= tolf) THEN ntesf = ntesf + 1 if (ntesf >= mtesf .AND. iters == 1) THEN iterm = 2 EXIT iteration END IF ELSE ntesf = 0 END IF ! Bundle updating CALL dobun(n,na,mal,x,g,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 xdiffy(n,g,gp,u) CALL aggre1(n,g,gp,ga,u,d,diag,alfn,alfv,gamma) ELSE CALL xdiffy(n,g,gp,u) CALL aggre2(n,g,gp,ga,d,diag,alfn,alfv,gamma) END IF CALL copy(n,xo,x) f = fo ELSE nnk = 0 CALL xdiffy(n,g,gp,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(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*g ! 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*g ! diag = gamma*I 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-6_prec) diag(i) = 1.0E-6_prec IF (diag(i) > told) diag(i) = told ELSE diag(i) = told END IF END DO CALL vxdiag(n,diag,g,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 ELSE diag(i) = told END IF END DO CALL vxdiag(n,diag,g,d) ! d = -one*d gamma = one CASE DEFAULT ! no scaling: diag = I ! no need to update diag or gamma d = -one*g ! END SELECT ELSE ! Update and direction determination for null steps. SELECT CASE(iscale) CASE(0,2) ! diag = 1 at null steps ! diag = one ! d = -one*ga ! diag = diag at null steps (skipping the update) ! no need to update diag or gamma d = -gamma*ga ! 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*ga ! diag = gamma CASE(1,3) ! diag = 1 at null steps ! diag = one ! d = -one*ga ! diag = diag at null steps (skipping the update) ! no need to update diag or gamma d = -gamma*ga ! 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*ga ! diag = gamma CASE(4,5) ! Herskovits's diagonal update CALL vxdiag(n,diag,ga,d) ! diag = diag at null steps. d = -one*d ! diag = one ! diag = 1 at null steps. ! d = -one*ga ! 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,ga,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*ga ! diag = gamma ! ELSE ! CALL vxdiag(n,diag,ga,d) ! diag = diag ! d = -one*d ! gamma = gamma_old ! END IF ! ELSE ! stu <= 0 ! CALL vxdiag(n,diag,ga,d) ! diag = diag ! d = -one*d ! END IF CASE DEFAULT ! no scaling: diag = I ! no need to update diag or gamma d = -one*ga ! END SELECT END IF END DO iteration 900 CONTINUE ! Printout the final results IF (iprint > 3) WRITE (6,FMT='(1X,''Exit from D-Bundle:'')') CALL wprint(iterm,iprint,nout) CALL rprint(n,nit,nfe,nge,x,f,xnorm,half*gnorm+alfv,iterm,iprint) CONTAINS SUBROUTINE restar(n,mcc,inew,ibun,iters,gp,g,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) :: & 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. 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,gp,g) iters = 1 nnk = 0 alfv=zero alfn=zero END IF ! gamma = large gamma = one d = -one*g diag = one END SUBROUTINE restar FUNCTION sclpar(iscale,sts,stu,utu) RESULT(spar) ! Calculation of the scaling parameter for D-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. ! 4 - Preliminary scaling with STU/UTU. ! tämmöisiä ei varmaan tarvi ! 5 - Preliminary scaling with STS/STU. ! tämmöisiä ei varmaan tarvi ! 6 - No scaling. !ok for Dbundle ! 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 < 1.0E-11_prec) spar = 1.0E-11_prec ! no need for correction ! IF (spar < 1.0E-5_prec) spar = 1.0E-5_prec ! tolg 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 SUBROUTINE dobun(n,ma,mal,x,g,f,ay,ag,af,iters,ibun) ! Bundle construction ! USE subpro, 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 dbundle !************************************************************************ !* * !* * 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 ! USE param, ONLY : zero,one ! these are the one needed in tinit itself ! 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. ! 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 lls * * !* * !* Special line search for limited memory bundle method * !* * !************************************************************************ SUBROUTINE lls(n,x,g,d,xo,t,fo,f,p,alfn,tmin,dnorm,wk,theta,epsl,epsr,& eta,iters,nfe,nge,nnk,iterm) USE param, ONLY : zero,half,one USE obj_fun, ONLY : & myf,myg ! Computation of the value and the subgradient 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) :: & g ! Subgradient 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) :: & 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. nge ! Number of subgradient evaluations. 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 Scalars INTEGER :: nin ! Number of interpolations. 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. ! 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 evaluation at a new point 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 myg(n,x,g,iterm) nge = nge + 1 IF (iterm /= 0) RETURN p = theta*vdot(n,g,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 lls !************************************************************************ !* * !* * SUBROUTINE update * * !* * !* Matrix update for D-bundle * !* * !************************************************************************ 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,g,gp,ga,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. g, & ! Current (auxiliary) subgradient of the objective function. gp, & ! Previous subgradient of the objective function. u ! Difference of trial and aggregate gradients. 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. 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 ga(i)=lam*g(i) + p*gp(i) END DO alfv = lam*alfn END SUBROUTINE aggre1 !************************************************************************ !* !* * SUBROUTINE aggre2 * !* !* Computation of aggregate values after consequtive null steps. !* !************************************************************************ SUBROUTINE aggre2(n,g,gp,ga,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. g, & ! Current (auxiliary) subgradient of the objective function. gp ! Previous subgradient of the objective function. 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. 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(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. tmp1, & ! Auxiliary scalar. tmp2, & ! Auxiliary scalar. tmp3 ! Auxiliary scalar. INTEGER :: i ! Intrinsic Functions INTRINSIC MIN,MAX CALL xdiffy(n,gp,ga,tmpn2) ! Calculation of tmpn3 = trans(gp - ga)dm CALL vxdiag(n,diag,tmpn2,tmpn3) 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) ! calculation of qr = trans(g - ga) dm (g - ga) 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 ga(i)=lam1*gp(i)+lam2*g(i)+tmp1*ga(i) END DO alfv = lam2*alfn + tmp1*alfv END SUBROUTINE aggre2 !************************************************************************ !* !* * 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 or subgradient 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 == 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 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,nit,nfe,nge,x,f,wk,qk,iterm,iprint) IMPLICIT NONE ! Scalar Arguments INTEGER, INTENT(IN) :: & n, & ! Number of variables nit, & ! Number of used iterations. nfe, & ! Number of used function evaluations. nge, & ! Number of used subgradient evaluations. 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,''nit='',I5,2X, & ''nfe='',I5,2X,''nge='',I5,2X,''f='',D15.8,2X,''wk='',D11.4,2X, & ''qk='',D11.4,2X)') nit,nfe,nge,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,''nit='',I5,2X, & ''nfe='',I5,2X,''nge='',I5,2X,''f='',D15.8,2X,''wk='',D11.4,2X, & ''qk='',D11.4,2X,''iterm='',I3)') nit,nfe,nge,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 dbundle_mod