* ************************************************************************* subroutine fdelhg( iel, elptr, elvar, xst, x, gptr, hptr, + fuval, negval, w1, inform ) * ************************************************************************* * Purpose : * --------- * This routine estimates the IELth element Hessian matrix * by forward differences in the element gradient values. * This IELth element Hessian has an elemental 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. * xst ( int ) * input : vector containing the status of the variables. * output : unmodified. * x ( dble ) * input : the current iterate vector. * output : unmodified. * 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. * 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. * 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, max, abs, sign. * gxfix, dcopy. * Programming : * ------------- * D. Tuyttens * ======================================================================== * Routine parameters integer iel, elptr(*), elvar(*), xst(*), gptr(*), + hptr(*), negval, inform double precision x(*), fuval(*), w1(*) * Internal variables integer i, j, k, ik, ii, ij, kk, ig, l, + debg, debh, eldim, istr, iend logical gxfix double precision sqrem, tempx, step * Save variables save i, j, k, ik, ii, ij, kk, ig, l, + debg, debh, eldim, istr, iend save sqrem, tempx, step * 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) eldim = elptr(iel+1) - elptr(iel) debg = gptr(iel) debh = hptr(iel) * call dcopy(eldim, fuval(debg), 1, w1, 1 ) ik = 0 istr = elptr(iel) iend = elptr(iel+1)-1 if ( iend .ge. istr ) then i = istr 10 continue k = elvar(i) ik = ik + 1 if( gxfix(xst(k)) ) then l = debh do 20 j = ik , eldim fuval(l) = zero l = l + 1 20 continue else tempx = x(k) step = sqrem * sign( max(abs(tempx),one) , tempx ) x(k) = tempx + step negval = negval + 1 inform = 4 return * * Element gradient evaluation. * 100 continue inform = 0 x(k) = tempx ig = debg do 30 ii = elptr(iel) , elptr(iel+1)-1 kk = elvar(ii) if( gxfix(xst(kk)) ) fuval(ig) = zero ig = ig + 1 30 continue l = debh do 40 j = 1 , ik-1 l = l - eldim + ik - j ij = ik - j fuval(l) = half*(fuval(l) + (fuval(debg+ij-1)-w1(ij))/step) 40 continue l = debh do 50 j = ik , eldim fuval(l) = (fuval(debg+j-1)-w1(j))/step l = l + 1 50 continue endif debh = debh + eldim - ik + 1 i = i + 1 if ( i .le. iend ) go to 10 endif call dcopy(eldim, w1, 1, fuval(debg), 1 ) * return end