MODULE smdb_mod ! Splitting Metrics D-Bundle Method USE r_precision, ONLY : prec ! Precision for reals. IMPLICIT NONE ! MODULE smdb_mod includes the following subroutines (S) and functions (F). PUBLIC :: & init_smdb, & ! S Initialization for SMDB subroutine. smdb ! S SMDB subroutine for nonsmooth large-scale ! optimization. Contains: ! S restar Initialization and reinitialization. ! S armijo Armijo line search. ! S nmarmijo Nonmonotone Armijo line search. ! S aggre1 Simplified subgradient aggregation. ! S aggre2 Subgradient aggregation. PRIVATE :: & wprint, & ! S Printout the error and warning messages. rprint ! S Printout the results. CONTAINS !*********************************************************************** !* * !* * SUBROUTINE init_smdb * * !* * !* Initialization for splitting metrics diagonal bundle method * !* for large-scale unconstrained nonsmooth optimization. * !* * !*********************************************************************** SUBROUTINE init_smdb(iterm) USE param, ONLY : zero,half,one,small,large USE initializat, ONLY : & n, & ! Number of variables. mcu, & ! Maximum number of stored corrections, mcu >= 1. iprint, & ! Printout specification (see initializat for more details). 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. 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,1) RETURN END IF IF (epsl >= 0.25_prec) THEN iterm = -5 CALL wprint(iterm,4) RETURN END IF IF (mcu <= 0) THEN iterm = -5 CALL wprint(iterm,2) RETURN END IF ! Default values IF (mit <= 0) mit = 10000 ! Maximum number of iterations. IF (mfe <= 0) mfe = n*mit ! Maximum number of function evaluations. IF (tolf <= zero) tolf = 1.0E-08_prec ! Tolerance for change of function values. IF (tolf2 == zero) tolf2 = 1.0E+04_prec ! Second tolerance for change of function values. IF (tolb == zero) tolb = -large + small ! Tolerance for the function value. IF (tolg <= zero) tolg = 1.0E-06_prec ! Tolerance for the termination criterion. IF (told <= zero) told = one ! Tolerance for the diagonal elements. IF (eta < zero) eta = 1.0E-03_prec ! 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. END SUBROUTINE init_smdb !*********************************************************************** !* * !* * SUBROUTINE smdb * * !* * !* Splitting metrics diagonal bundle subroutine * !* for nonconvex nonsmooth optimization. * !* * !*********************************************************************** SUBROUTINE smdb(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 mcu, & ! Maximum number of stored corrections, mcu >= 1. mnma, & ! Maximum number of function values used in nonmonotone line search. inma, & ! Selection of line search method: ! nma = 0, Armijo line search, ! nma = 1, nonmonotone Armijo line search. iprint, & ! Printout specification (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, & ! Upper bound for the diagonal elements of diag. mintold, & ! Lower bound for the diagonal elements of diag. 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 : & 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 objective function. ga, & ! Aggregate subgradient. s, & ! Difference of current and previous variables. u, & ! Difference of current and previous subgradients. d, & ! Direction vector. diag, & ! Diagonal matrix. cdiag, & ! "Convex" diagonal matrix. ncdiag, & ! "Concave" diagonal matrix. diagtmp ! Auxiliary diagonal matrix. REAL(KIND=prec), DIMENSION(n*mcu) :: & cssmatrix, & ! Matrix whose columns are stored differences of ! variables^2 when alfn >= 0. ncssmatrix, & ! Matrix whose columns are stored differences of ! variables^2 when alfn < 0. csumatrix, & ! Matrix whose columns are stored differences of ! variables and subgradients when alfn >= 0. ncsumatrix ! Matrix whose columns are stored differences of ! variables and subgradients when alfn < 0. REAL(KIND=prec), DIMENSION(mnma) :: & fold ! Old function values. ! Local Scalars REAL(KIND=prec) :: & alfn, & ! Linearization error. alfv, & ! Aggregate linearization error. 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. t, & ! Stepsize. fo, & ! Previous value of the objective. stu, & ! stu = trans(s)*u. utu, & ! utu = trans(u)*u. sts, & ! sts = trans(s)*s. ptmp, & ! Scalar for convex combination. a,b,c INTEGER :: i,j, & mccc, & ! Current number of stored "convex" corrections. mccnc, & ! Current number of stored "concave" corrections. inewc, & ! Index for the "convex" circular arrays. inewnc, & ! Index for the "concave" circular arrays. inewnma, & ! Index for the circular arrays containing function values. iters, & ! Null step indicator. ! 0 - Null step. ! 1 - Serious step. nnk, & ! Consecutive null steps counter. nser, & ! Serious steps counter. cnull, & ! Convex null steps counter. ncnull, & ! Concave null steps counter. neps, & ! Number of consecutive equal stopping criterions. ntesf, & ! Number of tests on function decrease. ncres, & ! Number of restarts. nres, & ! Number of consecutive restarts. nout ! Auxilary printout specification. ! Intrinsic Functions INTRINSIC ABS,MAX,SQRT,DOT_PRODUCT ! Computational Time REAL strtim,elapstim ! Parameters 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 SMD-Bundle:'')') ! Initialization inewnma = 1 nser = 0 cnull = 0 ncnull = 0 nout = 0 nit = 0 nfe = 0 nge = 0 ntesf = 0 nres = 1 ncres = -1 neps = 0 iterm = 0 iters = 1 nnk = 0 alfn = zero alfv = zero stu = one sts = one utu = one a = one fold = -large xnorm = large epsr = 0.25_prec+small IF (epsl+epsl >= epsr) THEN epsr = epsl+epsl + small IF (epsr >= half) THEN CALL wprint(iterm,-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,mccc,inewc,mccnc,inewnc,iters,gp,g,nnk,alfv,alfn,d,diag,cdiag,ncdiag,ncres) ! 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 = DOT_PRODUCT(g,g) dnorm = SQRT(DOT_PRODUCT(d,d)) p = DOT_PRODUCT(g,d) ELSE gnorm = DOT_PRODUCT(ga,ga) dnorm = SQRT(DOT_PRODUCT(d,d)) p = DOT_PRODUCT(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,-3) CALL restar(n,mccc,inewc,mccnc,inewnc,iters,gp,g,nnk,alfv,alfn,d,diag,cdiag,ncdiag,ncres) IF (ncres > maxnrs) THEN nout = maxnrs iterm = -2 EXIT iteration END IF CYCLE iteration END IF nout = -1 iterm = -1 EXIT iteration END IF ! Stopping criterion nit = nit + 1 pxnorm = xnorm IF (iters > 0) THEN CALL vxdiag(n,cdiag,g,diag) ! xnorm = DOT_PRODUCT(g,diag) + 2.0_prec*alfv ELSE CALL vxdiag(n,cdiag,ga,diag) ! xnorm = DOT_PRODUCT(ga,diag) + 2.0_prec*alfv END IF ! 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,-4) END IF CALL rprint(nit,nser,cnull,ncnull,nfe,nge,x,f,xnorm,iterm) ! Preparation of line search fo = f IF (iters > 0) THEN CALL copy2(n,x,xo,g,gp) END IF t = one ! Line search IF (inma == 0) THEN CALL armijo(x,g,d,xo,fo,f,p,alfn,xnorm,epsr,iters,nfe,nge,iterm) ELSE CALL nmarmijo(x,g,d,xo,fo,f,fold,p,alfn,xnorm,epsr,iters,nfe,nge,inewnma,iterm) END IF 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 ! Computation of variables and gradients differences CALL xdiffy(n,x,xo,s) CALL xdiffy(n,g,gp,u) ! Computation of aggregate values IF (iters == 0) THEN nnk = nnk + 1 IF (nnk == 1) THEN CALL aggre1(g,gp,ga,u,cdiag,alfn,alfv,dnorm) ELSE CALL aggre2(g,gp,ga,cdiag,alfn,alfv,dnorm) END IF CALL copy(n,xo,x) f = fo ELSE nnk = 0 END IF ! Direction finding IF (iters > 0) THEN ! Serious step update and "convex" direction determination nser = nser + 1 ! Serious step initialization alfv = zero ! Herskovits diagonal update IF (alfn >= zero) THEN CALL update(mccc,inewc,s,u,cssmatrix,csumatrix) ! convex update with sts and stu DO i = 1, n stu = zero sts = zero DO j = 1, mccc stu = stu + csumatrix((j-1)*n+i) sts = sts + cssmatrix((j-1)*n+i) END DO IF (stu > small) THEN cdiag(i) = sts/stu IF (cdiag(i) < mintold) cdiag(i) = mintold IF (cdiag(i) > told) cdiag(i) = told ELSE cdiag(i) = told END IF END DO ELSE CALL update(mccnc,inewnc,s,u,ncssmatrix,ncsumatrix) ! concave update DO i = 1, n stu = zero sts = zero DO j = 1, mccnc stu = stu + ncsumatrix((j-1)*n+i) sts = sts + ncssmatrix((j-1)*n+i) END DO IF (-stu > small) THEN ncdiag(i) = sts/stu IF (ncdiag(i) > -mintold) ncdiag(i) = -mintold IF (ncdiag(i) < -told) ncdiag(i) = -told ELSE ncdiag(i) = -told END IF END DO END IF ! Compute "convex" direction (we have a serious step, no need for "concave" direction). CALL vxdiag(n,cdiag,g,d) ! d = -one*d ELSE ! Update and direction determination for null steps. ! Herskovits's diagonal update IF (nnk == 1) THEN ! the first null step IF (alfn >= zero) THEN CALL update(mccc,inewc,s,u,cssmatrix,csumatrix) ! convex update with sts and stu DO i = 1, n stu = zero sts = zero DO j = 1, mccc stu = stu + csumatrix((j-1)*n+i) sts = sts + cssmatrix((j-1)*n+i) END DO IF (stu > small) THEN cdiag(i) = sts/stu IF (cdiag(i) < mintold) cdiag(i) = mintold IF (cdiag(i) > told) cdiag(i) = told ELSE cdiag(i) = told END IF END DO ELSE CALL update(mccnc,inewnc,s,u,ncssmatrix,ncsumatrix) ! concave update DO i = 1, n stu = zero sts = zero DO j = 1, mccnc stu = stu + ncsumatrix((j-1)*n+i) sts = sts + ncssmatrix((j-1)*n+i) END DO IF (-stu > small) THEN ncdiag(i) = sts/stu IF (ncdiag(i) > -mintold) ncdiag(i) = -mintold IF (ncdiag(i) < -told) ncdiag(i) = -told ELSE ncdiag(i) = -told END IF END DO END IF ELSE ! consecutive null step IF (alfn >= zero) THEN ! Skip the update or compute ! If you skip comment from here ... (preliminary results say skip always) ! DO i = 1, n ! stu = zero ! sts = zero ! DO j = 1, inewc-1 ! stu = stu + csumatrix((j-1)*n+i) ! sts = sts + cssmatrix((j-1)*n+i) ! END DO ! DO j = inewc+1, mccc ! stu = stu + csumatrix((j-1)*n+i) ! sts = sts + cssmatrix((j-1)*n+i) ! END DO ! stu = stu + s(i)*u(i) ! sts = sts + s(i)*s(i) ! ! IF (stu > small) THEN ! diagtmp(i) = sts/stu ! IF (diagtmp(i) < mintold) diagtmp(i) = mintold ! IF (diagtmp(i) > told) THEN ! diagtmp(i) = told ! END IF ! ELSE ! diagtmp(i) = told ! ! END IF ! ! END DO ! a=zero ! b=zero ! DO i=1,n ! c=ga(i)*ga(i) ! a=a+c*diagtmp(i) ! b=b+c*cdiag(i) ! END DO ! IF (a <= b) THEN ! CALL update(mccc,inewc,s,u,cssmatrix,csumatrix) ! convex update with sts and stu ! cdiag=diagtmp ! END IF ! ... to here ELSE CALL update(mccnc,inewnc,s,u,ncssmatrix,ncsumatrix) ! concave update DO i = 1, n stu = zero sts = zero DO j = 1, mccnc stu = stu + ncsumatrix((j-1)*n+i) sts = sts + ncssmatrix((j-1)*n+i) END DO IF (-stu > small) THEN ncdiag(i) = sts/stu IF (ncdiag(i) > -mintold) ncdiag(i) = -mintold IF (ncdiag(i) < -told) ncdiag(i) = -told ELSE ncdiag(i) = -told END IF END DO END IF END IF IF (alfn >= zero) THEN !convex null step cnull = cnull+1 CALL vxdiag(n,cdiag,ga,d) ! d = -one*d ELSE !concave null step ncnull = ncnull+1 ptmp = zero DO I=1,n IF (ABS(ncdiag(i)-cdiag(i)) > small) THEN ptmp = MAX(ncdiag(i)/(ncdiag(i)-cdiag(i)),ptmp) END IF END DO IF (ptmp > one) ptmp = one DO I=1,n diag(i) = ptmp*cdiag(i)+(one-ptmp)*ncdiag(i) END DO CALL vxdiag(n,diag,ga,d) d = -one*d END IF END IF END DO iteration 900 CONTINUE ! Printout the final results IF (iprint > 3) WRITE (6,FMT='(1X,''Exit from SMDB:'')') CALL wprint(iterm,nout) CALL rprint(nit,nser,cnull,ncnull,nfe,nge,x,f,xnorm,iterm) CONTAINS SUBROUTINE restar(n,mccc,inewc,mccnc,inewnc,iters,gp,g,nnk,alfv,alfn,d,diag,cdiag,ncdiag,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, & ! diag = I. cdiag, & ! cdiag = I. ncdiag ! ncdiag = -I. ! Scalar Arguments INTEGER, INTENT(IN) :: & n ! Number of variables (given in host) INTEGER, INTENT(OUT) :: & mccc, & ! mccnc, & ! inewc, & ! inewnc, & ! nnk ! Consecutive null steps counter. INTEGER, INTENT(INOUT) :: & iters, & ! Null step indicator. ! 0 - Null step. ! 1 - Serious step. ncres ! Number of restarts. REAL(KIND=prec), INTENT(OUT) :: & alfn, & ! Linearization error. alfv ! Aggregate linearization error. ! Restart mccc = 0 mccnc = 0 inewc = 1 inewnc = 1 ncres = ncres + 1 IF (iters == 0) THEN CALL copy(n,gp,g) iters = 1 nnk = 0 alfv=zero alfn=zero END IF d = -one*g diag = one cdiag = one ncdiag = -one END SUBROUTINE restar SUBROUTINE armijo(x,g,d,xo,fo,f,p,alfn,wk,epsr,iters,nfe,nge,iterm) ! Armijo line search ! USE param, ONLY : one,zero ! given in host USE initializat, ONLY : & ! n, & ! Number of variables. ! given in host ! epsl, & ! Descent parameter. ! given in host ! eta, & ! Distance measure parameter. ! given in host maxnin, & ! Maximum number of interpolations. sigma ! Tolerance for alfn. ! USE obj_fun, ONLY : & ! myf,myg ! Computation of the value and the subgradient of ! ! the objective function. ! given in host IMPLICIT NONE ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & d, & ! Direction vector. xo ! Previous vector of variables. REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: & x, & ! Vector of variables. g ! Subgradient of the objective function. ! Scalar Arguments REAL(KIND=prec), INTENT(INOUT) :: & epsr ! Linesearch parameter, used only in null step test to print an error message REAL(KIND=prec), INTENT(IN) :: & fo, & ! Previous value of the objective function. wk ! Stopping parameter. REAL(KIND=prec), INTENT(OUT) :: & f, & ! Value of the objective function. p, & ! Directional derivative. alfn ! Linearization error. 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 arrays REAL(KIND=prec), DIMENSION(n) :: & xu ! trial point ! Local Scalars REAL(KIND=prec) :: & fu, & ! Auxiliary values of the objective function a_orig, & ! Armijo parameter. inv_a_orig, & ! Inverse of Armijo parameter. alpha, & ! Armijo step size. epslwk,epsrwk ! Auxiliary scalars. INTEGER :: i,nin ! Number of interpolations. ! Intrinsic Functions INTRINSIC ABS,MAX,DOT_PRODUCT ! Initialization nin = 0 epslwk = epsl*wk epsrwk = epsr*wk alpha = 2.0_prec a_orig = alpha inv_a_orig = one/a_orig ! Function evaluation at a new point x = xo + d CALL myf(n,x,f,iterm) nfe = nfe + 1 IF (iterm /= 0) THEN f=fo RETURN END IF ! Serious step iters = 1 IF (f <= fo - epslwk) THEN ! descent condition xu = xo + alpha*d CALL myf(n,xu,fu,iterm) nfe = nfe + 1 IF (iterm /= 0) RETURN IF(fu > fo - alpha*epslwk .OR. maxnin==0) THEN ! Serious step with t=1 CALL myg(n,x,g,iterm) nge = nge + 1 IF (iterm /= 0) RETURN p = DOT_PRODUCT(g,d) alfn = fo-f+p ! Linearization error IF (alfn > -sigma .AND. alfn < zero) alfn=zero RETURN END IF inct: DO ! Increase t nin=nin+1 f=fu x=xu alpha = alpha * a_orig xu = xo + alpha*d CALL myf(n,xu,fu,iterm) nfe = nfe + 1 IF (iterm /= 0) RETURN IF (fu > fo - alpha*epslwk .OR. nin >= maxnin) THEN ! Serious step CALL myg(n,x,g,iterm) nge = nge + 1 IF (iterm /= 0) RETURN p = alpha * inv_a_orig * DOT_PRODUCT(g,d) alfn = fo-f+p ! Linearization error IF (alfn > -sigma .AND. alfn < zero) alfn=zero RETURN END IF END DO inct ELSE alpha = one dect: DO ! Decrease t nin=nin+1 IF (nin > maxnin) EXIT dect alpha = alpha * inv_a_orig ! x = xo + alpha * d CALL myf(n,x,f,iterm) nfe = nfe + 1 IF (iterm /= 0) RETURN IF (f <= fo - alpha*epslwk) THEN ! Serious step CALL myg(n,x,g,iterm) nge = nge + 1 IF (iterm /= 0) RETURN p = alpha * DOT_PRODUCT(g,d) alfn = fo-f+p ! Linearization error IF (alfn > -sigma .AND. alfn < zero) alfn=zero RETURN END IF END DO dect END IF iters = 0 ! Null step ! original d ! x = xo + d ! CALL myf(n,x,f,iterm) ! nfe = nfe + 1 CALL myg(n,x,g,iterm) ! x = xo + alpha*d nge = nge + 1 IF (iterm /= 0) RETURN ! p = DOT_PRODUCT(g,d) ! original d p = alpha * DOT_PRODUCT(g,d) alfn = fo-f+p ! Linearization error IF (alfn > -sigma .AND. alfn < zero) alfn=zero ! IF (p-alfn < -epsrwk) PRINT*,'not a good null step',p-alfn, -epsrwk ! IF (p-max(abs(alfn),eta*alpha*alpha*DOT_PRODUCT(d,d)) < -epsrwk) PRINT*,'not a good null step',p-max(abs(alfn),eta*alpha*alpha*DOT_PRODUCT(d,d)), -epsrwk END SUBROUTINE armijo SUBROUTINE nmarmijo(x,g,d,xo,fo,f,fold,p,alfn,wk,epsr,iters,nfe,nge,inewnma,iterm) ! Nonmonotone Armijo line search. ! USE param, ONLY : one,zero USE initializat, ONLY : & ! n, & ! Number of variables. ! given in host ! epsl, & ! Descent parameter. ! given in host ! eta, & ! Distance measure parameter. ! given in host ! mnma, & ! Maximum number of function values used in line search. ! given in host maxnin, & ! Maximum number of interpolations. sigma ! Tolerance for alfn. ! USE obj_fun, ONLY : & ! myf,myg ! Computation of the value and the subgradient of ! the objective function. ! given in host IMPLICIT NONE ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & d, & ! Direction vector. xo ! Previous vector of variables. REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: & x, & ! Vector of variables. g ! Subgradient of the objective function. REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: & fold ! Old function values. ! Scalar Arguments REAL(KIND=prec), INTENT(INOUT) :: & epsr ! Linesearch parameter, used only in null step test to print error REAL(KIND=prec), INTENT(IN) :: & fo, & ! Previous value of the objective function. wk ! Stopping parameter. REAL(KIND=prec), INTENT(OUT) :: & f, & ! Value of the objective function. p, & ! Directional derivative. alfn ! Linearization error. INTEGER, INTENT(INOUT) :: & nfe, & ! Number of function evaluations. nge, & ! Number of subgradient evaluations. iters, & ! Null step indicator. ! 0 - Null step. ! 1 - Serious step. inewnma ! Index for circular arrays. INTEGER, INTENT(OUT) :: & ! iters, & ! Null step indicator. ! 0 - Null step. ! 1 - Serious step. iterm ! Cause of termination: ! 0 - Everything is ok. ! -3 - Failure in function or subgradient calculations. ! Local arrays REAL(KIND=prec), DIMENSION(n) :: & xu ! trial point ! Local Scalars REAL(KIND=prec) :: & fu, & ! Auxiliary values of the objective function a_orig, & ! Armijo parameter. inv_a_orig, & ! Inverse of Armijo parameter. alpha, & ! Armijo step size. epslwk,epsrwk ! Auxiliary scalars. INTEGER :: i,nin ! Number of interpolations. ! Intrinsic Functions INTRINSIC ABS,MAX,MAXVAL,DOT_PRODUCT ! Initialization nin = 0 epslwk = epsl*wk epsrwk = epsr*wk alpha = 2.0_prec a_orig = alpha inv_a_orig = one/a_orig ! Updating circular array fold IF (iters==1) THEN fold(inewnma) = fo inewnma = inewnma + 1 IF (inewnma > mnma) inewnma = 1 END IF ! Function evaluation at a new point x = xo + d CALL myf(n,x,f,iterm) nfe = nfe + 1 IF (iterm /= 0) THEN f=fo RETURN END IF ! Serious step iters = 1 IF (f <= MAXVAL(fold) - epslwk) THEN ! Nonmonotone descent condition xu = xo + alpha*d CALL myf(n,xu,fu,iterm) nfe = nfe + 1 IF (iterm /= 0) RETURN IF(fu > fo - alpha*epslwk .OR. maxnin==0) THEN ! Serious step with t=1 CALL myg(n,x,g,iterm) nge = nge + 1 IF (iterm /= 0) RETURN p = DOT_PRODUCT(g,d) alfn = fo-f+p ! Linearization error IF (alfn > -sigma .AND. alfn < zero) alfn=zero RETURN END IF inct: DO ! Increase t nin=nin+1 f=fu x=xu alpha = alpha * a_orig xu = xo + alpha*d CALL myf(n,xu,fu,iterm) nfe = nfe + 1 IF (iterm /= 0) RETURN IF (fu > fo - alpha*epslwk .OR. nin >= maxnin) THEN ! Serious step CALL myg(n,x,g,iterm) nge = nge + 1 IF (iterm /= 0) RETURN p = alpha * inv_a_orig * DOT_PRODUCT(g,d) alfn = fo-f+p ! Linearization error IF (alfn > -sigma .AND. alfn < zero) alfn=zero RETURN END IF END DO inct ELSE alpha = one dect: DO ! Decrease t nin=nin+1 IF (nin > maxnin) EXIT dect alpha = alpha * inv_a_orig ! x = xo + alpha * d CALL myf(n,x,f,iterm) nfe = nfe + 1 IF (iterm /= 0) RETURN IF (f <= MAXVAL(fold) - alpha*epslwk) THEN ! Serious step CALL myg(n,x,g,iterm) nge = nge + 1 IF (iterm /= 0) RETURN p = alpha * DOT_PRODUCT(g,d) alfn = fo-f+p ! Linearization error IF (alfn > -sigma .AND. alfn < zero) alfn=zero RETURN END IF END DO dect END IF iters = 0 ! Null step ! original d ! x = xo + d ! CALL myf(n,x,f,iterm) ! nfe = nfe + 1 CALL myg(n,x,g,iterm) ! x = xo + alpha*d nge = nge + 1 IF (iterm /= 0) RETURN ! p = DOT_PRODUCT(g,d) ! original d p = alpha * DOT_PRODUCT(g,d) alfn = fo-f+p ! Linearization error IF (alfn > -sigma .AND. alfn < zero) alfn=zero ! IF (p-alfn < -epsrwk) PRINT*,'not a good null step',p-alfn, -epsrwk ! IF (p-max(abs(alfn),eta*alpha*alpha*DOT_PRODUCT(d,d)) < -epsrwk) PRINT*,'not a good null step',p-max(abs(alfn),eta*alpha*alpha*DOT_PRODUCT(d,d)), -epsrwk END SUBROUTINE nmarmijo SUBROUTINE update(mcc,inew,s,u,ssmatrix,sumatrix) ! Matrix update for SMDB ! USE subpro, ONLY : & ! copy2 ! Copying of two vectors. ! Given in host ! USE initializat, ONLY : & ! n, & ! Number of variables. ! Given in host ! mcu ! Upper limit for maximum number of stored corrections. ! Given in host 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(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(g,gp,ga,u,diag,alfn,alfv,dnorm) ! Computation of aggregate values after first null step. ! USE param, ONLY : zero,half,one ! Given in host ! USE subpro, ONLY : & ! Given in host ! vxdiag ! Vector multiplied by a diagonal matrix. ! USE initializat, ONLY : & ! Given in host ! n, & ! Number of variables. ! eta ! Distance measure parameter, eta >= 0. IMPLICIT NONE ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & diag, & ! Diagonal variable metric update. 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) :: & dnorm ! sqrt(norm(d)). REAL(KIND=prec), INTENT(INOUT) :: & alfn ! Locality measure. ! Local arrays REAL(KIND=prec), DIMENSION(n) :: tmpn ! Local Scalars REAL(KIND=prec) :: & tmp, & ! temporary storing of alfn. p, & ! trans(gp)*diag*u - alfn. q, & ! q = trans(u)*diag*u, where diag is the diagonal inverse ! approximation of the Hessian. lam ! Multiplier used to calculate aggregate values. INTEGER :: i ! Intrinsic Functions INTRINSIC ABS,MAX,MIN,SIGN,DOT_PRODUCT tmp=alfn alfn=MAX(ABS(alfn), eta*dnorm*dnorm) ! Computation of p=trans(gp)*diag*u - alfn and the product q=trans(u)*diag*u CALL vxdiag(n,diag,u,tmpn) p = - DOT_PRODUCT(tmpn,gp) - alfn q = DOT_PRODUCT(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 alfn=tmp END SUBROUTINE aggre1 SUBROUTINE aggre2(g,gp,ga,diag,alfn,alfv,dnorm) ! Computation of aggregate values after consequtive null steps. ! USE param, ONLY : zero,one,small ! Given in host ! USE subpro, ONLY : & ! Given in host ! vxdiag, & ! Vector multiplied by a diagonal matrix. ! xdiffy ! Difference of two vectors. ! USE initializat, ONLY : & ! Given in host ! n, & ! Number of variables. ! eta ! Distance measure parameter, eta >= 0. IMPLICIT NONE ! Array Arguments REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: & diag, & ! Diagonal variable metric update. 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) :: & alfn, & ! Locality measure. alfv ! Aggregate locality measure. REAL(KIND=prec), INTENT(IN) :: & dnorm ! sqrt(norm(d)). ! Local arrays REAL(KIND=prec), DIMENSION(n) :: tmpn2,tmpn3,tmpn4 ! Local Scalars REAL(KIND=prec) :: & tmp, & ! temporary storage for alfn. pr, & ! pr = trans(gp-ga) diag (gp-ga). rrp, & ! rrp = trans(gp-ga) diag ga - alfv. prqr, & ! prqr = trans(gp-ga) diag (g-ga). rrq, & ! rrq = trans(g-ga) diag ga - alfv + alfn. qr, & ! qr = trans(g-ga) diag (g-ga). pq, & ! pq = trans(g-gp) diag (g-gp). qqp, & ! qqp = trans(g-gp) diag 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 ABS,MIN,MAX,DOT_PRODUCT tmp=alfn alfn=MAX(ABS(alfn), eta*dnorm*dnorm) CALL xdiffy(n,gp,ga,tmpn2) ! Calculation of tmpn3 = trans(gp - ga)dm CALL vxdiag(n,diag,tmpn2,tmpn3) pr = DOT_PRODUCT(tmpn3,tmpn2) rrp = DOT_PRODUCT(tmpn3,ga) CALL xdiffy(n,g,ga,tmpn4) prqr = DOT_PRODUCT(tmpn3,tmpn4) ! calculation of qr = trans(g - ga) diag (g - ga) and rrq = trans(g-ga) diag ga + alfn - alfv. CALL vxdiag(n,diag,tmpn4,tmpn3) rrq = DOT_PRODUCT(tmpn3,ga) qr = DOT_PRODUCT(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) THEN 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 alfn=tmp RETURN END IF END IF END IF ! Minimum on the boundary 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 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 alfn=tmp END SUBROUTINE aggre2 END SUBROUTINE smdb !************************************************************************ !* !* * SUBROUTINE wprint * !* !* Printout the warning and error messages. !* !************************************************************************ SUBROUTINE wprint(iterm,nout) USE initializat, ONLY : & 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 IMPLICIT NONE ! Scalar Arguments INTEGER, INTENT(IN) :: & 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 == 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 == -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.'')') RETURN END IF ! Printout the final results 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 == 4) WRITE (6,FMT='(1X,''Abnormal exit: '' & ''Number of function evaluations > '',I5)') nout IF (iterm == 5) WRITE (6,FMT='(1X,''Abnormal exit: '' & ''Number of iterations > '',I5)') nout 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 == 8) WRITE (6,FMT='(1X,''Abnormal exit: '' & ''Failure in attaining the demanded accuracy.'')') IF (iterm == -1) WRITE (6,FMT='(1X,''Abnormal exit: '' & ''Two consecutive restarts.'')') 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.'')') END IF END IF END SUBROUTINE wprint !************************************************************************ !* !* * SUBROUTINE rprint * !* !* Printout the (final) results. !* !************************************************************************ SUBROUTINE rprint(nit,nser,cnull,ncnull,nfe,nge,x,f,wk,iterm) USE initializat, ONLY : & n, & ! Number of variables 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 IMPLICIT NONE ! Scalar Arguments INTEGER, INTENT(IN) :: & nit, & ! Number of used iterations. nser, & ! Number of serious steps. cnull, & ! Number of convex null steps. ncnull, & ! Number of nonconvex null steps. nfe, & ! Number of used function evaluations. nge, & ! Number of used subgradient evaluations. 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 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 & )') nit,nfe,nge,f,wk 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,''nser='',I5,2X,''cnull='',I5,2X, & ''ncnull='',I5,2X,''nfe='',I5,2X,''nge='',I5,2X,''f='',D15.8,2X,''wk='',D11.4,2X, & ''iterm='',I3)') nit,nser,cnull,ncnull,nfe,nge,f,wk,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 smdb_mod