* ************************************************************************* subroutine update( inform, x, xst, elst, elptr, elvar, + fuval, fcalc, lfc, olgra, gptr, + lg, hptr, nefval, negval, nehval, + dx, s, y, w4, recalc, prec, dpr, + nbfgs, nrk1, noupd ) * ************************************************************************* * Purpose : * --------- * This routine updates the element gradient vectors and the * element Hessians matrices that are stored in array FUVAL. * A considerable flexibility in the computation of the * derivatives is provided in LSNNO. If the element gradient * and/or the element Hessian are not supplied by the user, * the software contains an assortment of routines for * computing suitable approximations of these quantities. * The way these quantities will be computed is indicated * in array ELST that contains the status of the elemnt * functions. * Every time informations on some element function values or * on some element gradient values or on some element Hessian * values are needed, a return is made to the main program. * This routine is then re-entered with the information available * and the calculation proceeds. * Parameters : * ------------ * inform ( int ) * input : must be set to zero the first time this routine * is called. * If it is not equal to zero, the routine is then * re-entered with some information that was needed * to proceed the calculation. * output : If it is not equal to zero, some element function * values or some element gradient values or some * element Hessian values are needed to proceed the * calculation. A return is then made to the main * program. * Otherwise, no further information is needed. * x ( dble ) * input : the current iterate vector. * output : unmodified. * xst ( int ) * input : vector containing the status of the variables. * output : unmodified. * elst ( int ) * input : vector containing the status of the element * functions. * output : the status of the element functions updated * by a quasi-Newtom formula may have been * changed. * elptr ( int ) * input : array whose kth value is the position of * the first variable of element k, in the * list ELVAR. * output : unmodified. * elvar ( int ) * input : array containing the indices of the varaibles * in the first element, followed by those in the * second element, etc. * output : unmodified. * fuval ( dble ) * input : array used to store the function and derivative * values of the element functions. * output : the derivative values of the element functions * are updated. * fcalc ( int ) * input : meaningless. * output : this vector contains the indices of the element * functions for which information is required. * lfc ( int ) * input : meaningless. * output : the number of element functions for which * information is required. * olgra ( dble ) * input : this array contains the element gradient * vectors evaluated at the previous iterate. * output : it contains the differences between the * element gradient vectors evaluated at point * X and the element gradient vectors evaluated * at the previous iterate. * gptr ( int ) * input : array whose ith value is the position of the * first component of the ith element gradient * in FUVAL. * output : unmodified. * lg ( int ) * input : the total lenght of the element gradient vectors. * output : unmodified. * hptr ( int ) * input : array whose ith value is the position of the * first component of the ith element Hessian * in FUVAL. * output : unmodified. * nefval ( int ) * input : the number of element function evaluations. * output : this number is increased by one each time * an element function needs to be evaluated. * negval ( int ) * input : the number of element gradient evaluations. * output : this number is increased by one each time * an element gradient needs to be evaluated. * nehval ( int ) * input : the number of element Hessian evaluations. * output : this number is increased by one each time * an element Hessian needs to be evaluated. * dx ( dble ) * input : array used as workspace. * output : it contains the differences between the * current iterate X and the previous iterate. * s ( dble ) * input : array used as workspace. * output : it contains for each element function the * differences between the current iterate X * and the previous iterate. * y ( dble ) * input : array used as workspace. * output : it contains for each element function the * differences between the element gradient * vector evaluated at the current iterate X * and the element gradient vector evaluated * at the previous iterate. * w4 ( dble ) * input : array used as workspace. * output : meaningless. * recalc ( int ) * input : array containing the indices of the elements whose * derivatives need updating. * output : unmodified. * prec ( log ) * input : .true. if the user wants to use a preconditioned * conjugate gradient scheme to solve the linear * systems of equations which arise at each * iteration of the minimization procedure. * .false. if the user does not want to use a * preconditioner. * output : unmodified. * dpr ( dble ) * input : array containing the information necessary to build * the preconditioner associated with Dembo's strategy. * output : the updated information. * nbfgs ( int ) * input : meaningless. * output : the number of element Hessians that have been * updated by the BFGS formula. * nrk1 ( int ) * input : meaningless. * output : the number of element Hessians that have been * updated by the Rank-One formula. * noupd ( int ) * input : meaningless. * output : the number of element Hessians for which * no update was needed. * Routines used : * --------------- * mod, sqrt, abs, max. * dcopy, ddot, analg, elint, fdingf, fdelgf, * fdinhg, fdelhg, fdinhf, fdelhf, range, gathr0, * hessid, selhqn, elcvx, scbfgs, scrk1. * Programming : * ------------- * D. Tuyttens * ======================================================================== * Routine parameters integer inform, xst(*), elst(*), elptr(*), + elvar(*), gptr(*), lg, hptr(*), nefval, + negval, nehval, fcalc(*), lfc, nbfgs, + nrk1, noupd, recalc(*) logical prec double precision x(*), fuval(*), dx(*), s(*), y(*), + w4(*), olgra(*), dpr(*) * Internal variables integer iel, els, dim, eldim, intdim, i, ij, k, + debx, debg, debh, ifc, nrec logical analg, elint, elhqn, elcvx, newp, + upd, negcur, qnd double precision sbs, yts, alpha, sqrem, ddot * Save variables save iel, els, dim, eldim, intdim, i, ij, k, + debx, debg, debh, ifc, nrec save newp, upd, negcur, qnd save sbs, yts, alpha, sqrem * Common specifications integer arcs, nodes, elem double precision epsmch, huge, tiny, tol common / prbdim / arcs, nodes, elem common / prbmch / epsmch, huge, tiny, tol double precision zero, one, two, three, half, tenm1, tenm2, tenm4 common / prbcst / zero, one, two, three, half, tenm1, tenm2, tenm4 * * Reverse communication test. * if( inform.ne.0 ) then inform = inform - 14 go to( 200, 910, 920, 500, 500, + 300, 300, 700, 700, 700, + 700, 700, 700 ) inform endif * * We obtain from the informations stored in array * ELST, the way the element gradients and element * Hessians will be estimated. * nrec = lfc call dcopy( lg, fuval(gptr(1)), 1, olgra, 1 ) * * * The gradient and Hessian are provided by the user. * qnd = .false. lfc = 0 do 100 ifc = 1, nrec iel = recalc( ifc ) els = mod( elst(iel), 8 ) if( analg(elst(iel)) .and. els.eq.0 ) then lfc = lfc + 1 fcalc(lfc) = iel endif qnd = qnd .or. els .ge. 3 100 continue * * Element gradient and Hessian evaluations. * if( lfc .gt. 0 ) then negval = negval + lfc nehval = nehval + lfc inform = 2 return endif 910 continue * * Only the gradient is provided by the user. * lfc = 0 do 210 ifc = 1, nrec iel = recalc( ifc ) els = mod( elst(iel), 8 ) if( analg(elst(iel)) .and. els.ne.0 ) then lfc = lfc + 1 fcalc(lfc) = iel endif 210 continue * * Element gradient evaluations. * if( lfc .gt. 0 ) then negval = negval + lfc inform = 3 return endif 920 continue * * Only the Hessian is provided by the user. * lfc = 0 do 120 ifc = 1, nrec iel = recalc( ifc ) els = mod( elst(iel), 8 ) if( .not. analg(elst(iel)) .and. els.eq.0 ) then lfc = lfc + 1 fcalc(lfc) = iel endif 120 continue * * Element Hessian evaluations. * if( lfc .gt. 0 ) then nehval = nehval + lfc inform = 1 return endif 200 continue inform = 0 * * The gradient vector is estimated by * differences in the function values. * lfc = 1 if( nrec .ge. 1 ) then ifc = 1 10 continue iel = recalc(ifc) if( .not. analg(elst(iel)) ) then fcalc(1) = iel 300 continue if( elint(elst(iel)) ) then * * The element gradient has an internal representation. * call fdingf( iel, elptr, elvar, x, w4(arcs+1), s, y, + gptr, fuval, nefval, inform ) else * * The element gradient has an elemental representation. * call fdelgf( iel, elptr, elvar, xst, x, + gptr, fuval, nefval, inform ) endif if( inform.ne.0 ) return * * Element function evaluations. * endif ifc = ifc + 1 if ( ifc .le. nrec ) go to 10 endif * * The Hessian matrix is estimated by * differences in the gradient values. * lfc = 1 if( nrec .ge. 1 ) then ifc = 1 40 continue iel = recalc(ifc) if( mod( elst(iel), 8 ) .eq. 1 ) then fcalc(1) = iel 500 continue if( elint(elst(iel)) ) then * * The element Hessian has an internal representation. * call fdinhg( iel, elptr, elvar, x, w4(arcs+1), s, y, w4, + fuval, gptr, hptr, negval, inform ) else * * The element Hessian has an elemental representation. * call fdelhg( iel, elptr, elvar, xst, x, gptr, hptr, + fuval, negval, w4(arcs+1), inform ) endif if( inform.ne.0 ) return * * Element gradient evaluations. * endif ifc = ifc + 1 if ( ifc .le. nrec ) go to 40 endif * * The Hessian matrix is estimated by * differences in the function values. * lfc = 1 if( nrec .ge. 1 ) then ifc = 1 50 continue iel = recalc(ifc) if( mod( elst(iel), 8 ) .eq. 2 ) then fcalc(1) = iel 700 continue if( elint(elst(iel)) ) then * * The element Hessian has an internal representation. * call fdinhf( iel, elptr, elvar, x, w4(arcs+1), s, y, + w4, hptr, fuval, nefval, inform ) else * * The element Hessian has an elemental representation. * call fdelhf( iel, elptr, elvar, xst, x, w4(arcs+1), + s, hptr, fuval, nefval, inform ) endif if( inform.ne.0 ) return * * Element function evaluations. * endif ifc = ifc + 1 if ( ifc .le. nrec ) go to 50 endif * * The element Hessian matrices are updated by a quasi-Newton formula * or Dembo's preconditioner needs computing * if( qnd ) then sqrem = sqrt(epsmch) nbfgs = 0 nrk1 = 0 noupd = 0 * * We compute the differences between the * iterate X and the previous iterate. * do 70 k = 1 , arcs dx(k) = x(k) - dx(k) 70 continue * * We compute the differences between the gradient * evaluated at point X and the gradient evaluated * at the previous iterate. * ij = gptr(1) do 75 i = 1 , lg olgra(i) = fuval(ij) - olgra(i) ij = ij + 1 75 continue * * We loop on the element functions * that need to be updated. * do 80 ifc = 1 , nrec iel = recalc(ifc) els = mod( elst(iel), 8 ) * * Dembo's strategy without preconditioning. * In this case, nothing needs updating. * if( els.eq.7 .and. .not.prec ) goto 80 * eldim = elptr(iel+1) - elptr(iel) debx = elptr(iel) debg = gptr(iel) - elem debh = hptr(iel) * * We obtain the dimension of the * element Hessian matrix. * if( elint(elst(iel)) ) then call range( iel, eldim, intdim, s, y , 0 ) dim = intdim else dim = eldim endif * * We compute the vectors S and Y. * if( dim .lt. eldim ) then if( els.eq.7 .and. prec ) then call dcopy( dim, olgra(debg), 1, s, 1 ) call range( iel, eldim, intdim, s, y, 2 ) call gathr0( eldim, s, 1, dx, elvar(debx) ) else call gathr0( eldim, y, 1, dx, elvar(debx) ) call range( iel, eldim, intdim, y, s, 1 ) call dcopy( dim, olgra(debg), 1, y, 1 ) endif else call gathr0( dim, s, 1, dx, elvar(debx) ) call dcopy( dim, olgra(debg), 1, y, 1 ) endif * * We test if this element Hessian matrix needs * really to be updated. * newp = .false. do 90 i = 1 , dim newp = newp .or. abs(s(i)).gt.epsmch 90 continue * if( newp ) then * * The element Hessian matrix needs * to be updated. We compute its * curvate estimate. * yts = ddot( dim, s, 1, y, 1 ) * * The Hessian matrix is set to a scaling of * the identity matrix before its is first * updated by a quasi-Newton formula. * if( els.eq.6 .and. .not.elhqn(elst(iel)) ) then sbs = ddot( dim, s, 1, s, 1 ) if( sbs.gt.sqrem ) then alpha = max( sqrem , abs(yts / sbs) ) else alpha = one endif call hessid( dim, fuval(debh), alpha ) call selhqn( elst(iel) ) endif * * When Dembo's preconditioner is used, * we update the diagonal of the Hessian * matrix and we store it in array FUVAL. * if( els.eq.7 .and. prec ) then debh = arcs + lg + elptr(iel) do 95 ij = 1 , eldim if( abs(s(ij)).gt.sqrem ) then dpr(debh) = y(ij) / s(ij) else dpr(debh) = one endif debh = debh + 1 95 continue else if( dim.eq.1 ) then * * The element function has * only one variable. * fuval(debh) = y(1) / s(1) * else if( elcvx(elst(iel)) ) then * * We test if the element function is convex. * if( yts.le.zero ) go to 110 * * The BFGS update is applied. * call scbfgs( dim, fuval(debh), s, y, epsmch, + epsmch, negcur, w4 ) * * We test if a negative curvate has been observed. * if( negcur ) then * * A negative curvate has been detected. * We will used the Rank-One update. * go to 110 else nbfgs = nbfgs + 1 endif else * * The Rank-One update is applied. * 110 continue call scrk1( dim, fuval(debh), s, y, epsmch, + epsmch, upd, w4 ) * if( upd ) then * * We update the status of the element * function to correspond to the status * of an element function that is non * convex. * call snelcv( elst(iel) ) nrk1 = nrk1 + 1 else * * No update has been done. * noupd = noupd + 1 endif endif endif 80 continue endif * 900 continue inform = 0 * return end