* ************************************************************************* subroutine fdinhg( iel, elptr, elvar, x, w1, w2, w3, w4, + fuval, gptr, hptr, negval, inform ) * ************************************************************************* * Purpose : * --------- * This routine estimates the IELth element Hessian matrix * by forward differences in the element gradient values. * This IELth element Hessian has an internal representation. * Every time information on the element gradient value is needed, * a return is made to the main program. This routine is then * re-entered with the information available and the calculation * proceeds. * Parameters : * ------------ * iel ( int ) * input : the indice of the element Hessian matrix * that will be estimated. * output : unmodified. * 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. * x ( dble ) * input : the current iterate vector. * output : unmodified. * w1 ( dble ) * input : array used as workspace. * output : meaningless. * w2 ( dble ) * input : array used as workspace. * output : meaningless. * w3 ( dble ) * input : array used as workspace. * output : meaningless. * w4 ( dble ) * input : array used as workspace. * output : meaningless. * fuval ( dble ) * input : array used to store the function and derivative * values of the element functions. * output : the components of indices HPTR(IEL) up to * HPTR(IEL+1)-1 will be set to the estimate * of the element Hessian matrix by forward * differences in the element gradient values. * gptr ( int ) * input : array whose ith value is the position of the * first component of the ith element gradient * in FUVAL. * 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. * negval ( int ) * input : the number of element gradient evaluations. * output : this number is increased by one each time * one element gradient needs to be evaluated. * 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 re-entered * with an element gradient value that was needed * to proceed the calculation. * output : If it is not equal to zero, an element gradient * value is needed to proceed the calculation. A return * is then made to the main program. * Otherwise, no further information is needed. * Routines used : * --------------- * sqrt. * range, dcopcd, dcopy, dsetvl, dxpycd. * Programming : * ------------- * D. Tuyttens * ======================================================================== * Routine parameters integer iel, elptr(*), elvar(*), gptr(*), + hptr(*), negval, inform double precision x(*), w1(*), w2(*), w3(*), w4(*), + fuval(*) * Internal variables integer i, j, ij, l, eldim, intdim, debg, + debh double precision sqrem * Save variables save i, j, ij, l, eldim, intdim, debg save debh save sqrem * Common specifications double precision epsmch, huge, tiny, tol 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 ) goto 100 * * Estimate of the IELth element Hessian matrix by * forward differences in the element gradient values. * sqrem = sqrt(epsmch) debg = gptr(iel) debh = hptr(iel) eldim = elptr(iel+1) - elptr(iel) call range( iel, eldim, intdim, w1, w1, 0 ) call dcopcd( eldim, elvar(elptr(iel)), x, w3 ) call dcopy( intdim, fuval(debg), 1, w4, 1 ) call dsetvl( intdim, w1, 1, zero ) * i = 1 10 continue w1(i) = sqrem call range( iel, eldim, intdim, w1, w2, 3 ) call dxpycd( eldim, elvar(elptr(iel)), x, w2 ) negval = negval + 1 inform = 5 return * * Element gradient evaluation. * 100 continue inform = 0 call dcopcd( eldim, elvar(elptr(iel)), w3, x ) l = debh do 20 j = 1 , i-1 l = l - intdim + i - j ij = i - j fuval(l) = half*(fuval(l) + (fuval(debg+ij-1)-w4(ij))/sqrem) 20 continue l = debh do 30 j = i , intdim fuval(l) = (fuval(debg+j-1)-w4(j))/sqrem l = l + 1 30 continue debh = debh + intdim - i + 1 i = i + 1 if( i .le. intdim ) go to 10 call dcopy( intdim, w4, 1, fuval(debg), 1 ) * return end