* ************************************************************************* subroutine fdinhf( iel, elptr, elvar, x, w1, w2, w3, w4, + hptr, fuval, nefval, inform ) * ************************************************************************* * Purpose : * --------- * This routine estimates the IELth element Hessian matrix * by forward differences in the element function values. * This IELth element Hessian has an internal representation. * Every time information on the element function 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. * hptr ( int ) * input : array whose ith value is the position of the * first component of the ith element Hessian * in FUVAL. * output : unmodified. * 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 function values. * nefval ( int ) * input : the number of element function evaluations. * output : this number is increased by one each time * one element function 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 function value that was needed * to proceed the calculation. * output : If it is not equal to zero, an element function * value is needed to proceed the calculation. A return * is then made to the main program. * Otherwise, no further information is needed. * Routines used : * --------------- * range, dcopcd, dsetvl, dxpycd. * Programming : * ------------- * D. Tuyttens * ======================================================================== * Routine parameters integer iel, elptr(*), elvar(*), hptr(*), nefval, + inform double precision x(*), w1(*), w2(*), w3(*), w4(*), fuval(*) * Internal variables integer i, j, debh, eldim, intdim double precision thrdem, step2, tempf * Save variables save i, j, debh, eldim, intdim save thrdem, step2, tempf * 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.eq.11 ) goto 100 if( inform.eq.12 ) goto 200 if( inform.eq.13 ) goto 300 * * Estimate of the IELth element Hessian matrix by * forward differences in the element function values. * eldim = elptr(iel+1) - elptr(iel) thrdem = epsmch**(one/three) step2 = thrdem**2 call range(iel, eldim, intdim, w1, w2, 0 ) call dcopcd(eldim, elvar(elptr(iel)), x, w4 ) call dsetvl( intdim, w1, 1, zero ) * tempf = fuval(iel) i = 1 10 continue w1(i) = thrdem call range( iel, eldim, intdim, w1, w2, 3 ) w1(i) = zero call dxpycd( eldim, elvar(elptr(iel)), x, w2 ) nefval = nefval + 1 inform = 11 return * * Element function evaluation. * 100 continue w3(i) = fuval(iel) inform = 0 call dcopcd(eldim, elvar(elptr(iel)), w4, x ) i = i + 1 if ( i .le. intdim ) go to 10 * debh = hptr(iel) i = 1 20 continue w1(i) = thrdem j = i+1 if( j .le. intdim ) then 30 continue w1(j) = thrdem call range( iel, eldim, intdim, w1, w2, 3 ) w1(j) = 0.0d0 call dxpycd( eldim, elvar(elptr(iel)), x, w2 ) nefval = nefval + 1 inform = 12 return * * Element function evaluation. * 200 continue inform = 0 call dcopcd(eldim, elvar(elptr(iel)), w4, x ) fuval(debh+j-i) = ((tempf-w3(j))+(fuval(iel)-w3(i))) / step2 j = j + 1 if( j .le. intdim ) go to 30 endif * w1(i) = two * thrdem call range( iel, eldim, intdim, w1, w2, 3 ) w1(i) = zero call dxpycd( eldim, elvar(elptr(iel)), x, w2 ) nefval = nefval + 1 inform = 13 return * * Element function evaluation. * 300 continue inform = 0 call dcopcd(eldim, elvar(elptr(iel)), w4, x ) fuval(debh) = ((tempf-w3(i))+(fuval(iel)-w3(i))) / step2 i = i + 1 if ( i .le. intdim ) go to 20 fuval(iel) = tempf * return end