* ************************************************************************* subroutine prod( x, fuval, gptr, lg, negval, hptr, fcalc, lfc, + inform, curv, ds, gds, dw1, dw2, elvar, + elptr, elst, ba, nba, su, nsu, bep, beptr, + ftstem, kdist, bgx ) * ************************************************************************* * Purpose : * --------- * This routine computes the product of the reduced Hessian * matrix times the conjugate gradient direction ds and store * it in array gds. * t * gds = Z G Z ds * This matrix vector product can be obtained by calling a * special procedure that computes the product of a symmetric * partially separable matrix times a vector, or by computing * differences in the element gradient values. * Parameters : * ------------ * x ( dble ) * input : the current iterate vector. * output : unmodified. * fuval ( dble ) * input : array used to store the function and derivative * values of the element functions. * 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. * lg ( int ) * input : the total lenght of the element gradient vectors. * output : unmodified. * 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. * hptr ( int ) * input : array whose ith value is the position of the * first component of the ith element Hessian * in FUVAL. * output : unmodified. * fcalc ( int ) * input : meaningless. * output : this vector contains the indices of the element * Hessian matrices involved in the matrix vector * product. * lfc ( int ) * input : meaningless. * output : the number of element Hessian matrices involved * in the matrix vector product. * 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 gradient * values are needed to proceed the calculation. A * return is then made to the main program. * Otherwise, no further information is needed. * curv ( dble ) * input : meaningless. * output : the curvature of the conjugate gradient * direction ds. * ds ( dble ) * input : the superbasic components of the conjugate * gradient direction. * output : the basic corresponding components of the * conjugate gradient direction are computed. * gds ( dble ) * input : meaningless * output : it contains the product of the reduced Hessian * matrix times the conjugate gradient direction ds. * dw1 ( dble ) * input : array used as workspace. * output : meaningless. * dw2 ( dble ) * input : array used as workspace. * output : meaningless. * elvar ( int ) * input : array containing the indices of the varaibles * in the first element, followed by those in the * second element, etc. * 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. * elst ( int ) * input : vector containing the status of the element * functions. * output : unmodified. * ba ( int ) * input : vector containing the indices of the basic * variables. * output : unmodified. * nba ( int ) * input : the number of basic variables. * output : unmodified. * su ( int ) * input : vector containing the indices of the superbasic * variables. * output : unmodified. * nsu ( int ) * input : the number of superbasic variables. * output : unmodified. * bep ( int ) * input : array containing the flow augmenting paths of * the superbasic variables. * output : unmodified. * beptr ( int ) * input : array whose kth value is the position of the * first element of the flow augmenting path of * superbasic arc k, in array BEP. * output : unmodified. * ftstem ( int ) * input : vector containing for each superbasic variable, * the length from its origine node to the stem * node, and the length of its flow augmenting path. * output : unmodified. * kdist ( int ) * input : vector containing for each superbasic variable * the length of its flow augmenting path. * output : unmodified. * bgx ( dble ) * input : array with the gradient at x. * output : unmodified (only used for Dembo's strategy). * Routines used : * --------------- * mod. * dsetcd, dsetvl, dembo, psprdt, ddotcd. * Programming : * ------------- * D. Tuyttens * ======================================================================== * Routine parameters integer elvar(*), elptr(*), gptr(*), lg, negval, + hptr(*), elst(*), ba(*), nba, su(*), nsu, + bep(*), beptr(*), ftstem(*), inform, + fcalc(*), lfc, kdist(*) double precision x(*), fuval(*), curv, ds(*), gds(*), dw1(*), + dw2(*), bgx(*) * Internal variables integer ik, k, kk, kdis, i, j, iel, ndmb double precision ddotcd * Save variables save ik, k, kk, kdis, i, j, iel, ndmb * Common specifications integer arcs, nodes, elem common / prbdim / arcs, nodes, elem 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 * * On input, vector ds contains only the superbasic * components of the conjugate gradient direction. * We now compute the corresponding basic components * of the conjugate gradient direction and store them * in vector ds (i.e. ds = Z ds ). * call dsetcd( nba, ba, ds, zero ) do 10 ik = 1 , nsu k = su(ik) kdis = kdist(ik) * * We compute the corresponding basic components * of the conjugate gradient direction. * do 20 i = beptr(k) , beptr(k)+kdis-1 kk = bep(i) if( kk.lt.0 ) then ds(-kk) = ds(-kk) - ds(k) else ds(kk) = ds(kk) + ds(k) endif 20 continue 10 continue * * We search the indices of the element Hessians * that are involved in the matrix vector product. * lfc = 0 ndmb = 0 do 30 iel = 1 , elem do 40 j = elptr(iel) , elptr(iel+1)-1 if( ds(elvar(j)).ne.zero ) then lfc = lfc + 1 if( mod(elst(iel),8) .eq. 7 ) then ndmb = ndmb + 1 fcalc(lfc) = -iel else fcalc(lfc) = iel endif go to 30 endif 40 continue 30 continue * * We compute the product G ds = G (Z ds), * and we store it in array gds. * call dsetvl( arcs, gds, 1, zero ) * * The matrix vector product is obtained by calling a special * procedure that computes the product of a symmetric partially * separable matrix times a vector. * if ( ndmb .lt. lfc ) then call psprdt( fcalc, lfc, fuval, hptr, ds, gds, elvar, + elptr, elst, dw1, dw2 ) endif * * The matrix vector product is obtained by computing * differences in the element gradient values. * if ( ndmb .gt. 0 ) then 100 continue call dembo( inform, elptr, elvar, elst, gptr, lg, negval, + fcalc, lfc, fuval, x, dw1(arcs+1), ds, gds, + dw1, dw2, bgx(arcs+1) ) if( inform.ne.0 ) return endif * * t t * We compute the product Z G Z ds = Z gds * and we store it in array gds. * do 50 ik = 1 , nsu k = su(ik) kdis = kdist(ik) do 60 i = beptr(k) , beptr(k)+kdis-1 kk = bep(i) if( kk.lt.0 ) then gds(k) = gds(k) - gds(-kk) else gds(k) = gds(k) + gds(kk) endif 60 continue 50 continue * * We compute the curvate of the conjugate * gradient direction ds. * curv = ddotcd( nsu, su, gds, ds ) * return end