* ************************************************************************* subroutine prcond( fuval, lg, hptr, diag, dpr, w1, w2, elptr, + elvar, elst, iss, ip, su, nsu, nqa, + ftstem, beptr, bep, kdist ) * ************************************************************************* * Purpose : * --------- * This routine computes the diagonal preconditioner that is * used in the preconditioned conjugate gradient method. This * preconditioner is stored in array dpr and is computed as * follows : * t * dpr = | diag( Z diag(G) Z ) | * * where diag(G) is the diagonal of matrix G. * When Dembo's method is used, the Hessian matrix is not * stored at all and the diagonal of this Hessian matrix * is estimated as follows : * (g(x) - g(x_))[i] * (diag(G))[i] = ----------------- * ( x - x_ )[i] * Parameters : * ------------ * fuval ( dble ) * input : array used to store the function and derivative * values of the element functions. * 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. * diag ( dble ) * input : meaningless. * output : it contains the diagonal elements of the * Hessian matrix, for the variables that are * included in the independent set IP. * dpr ( dble ) * input : meaningless. * output : it contains the diagonal preconditioner * corresponding to independent set IP. * w1 ( dble ) * input : array used as workspace. * output : meaningless. * w2 ( dble ) * input : array used as workspace. * output : meaningless. * 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. * elst ( int ) * input : vector containing the status of the element * functions. * output : unmodified. * iss ( int ) * input : vector containing for each variable, the indice * of the independent set in which it is included. * It is equal to zero when the variable does not * belongs to any independent set. * output : unmodified. * ip ( int ) * input : the indice of the current independent set. * 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. * nqa ( int ) * input : the number of superbasic variables that are * quasi active. * 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. * 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. * bep ( int ) * input : array containing the flow augmenting paths of * the superbasic variables. * output : unmodified. * kdist ( int ) * input : vector containing for each superbasic variable * the length of its flow augmenting path. * output : unmodified. * Routines used : * --------------- * mod, sqrt, abs. * dsetvl, elint, range, ltcmvp. * Programming : * ------------- * D. Tuyttens * ======================================================================== * Routine parameters integer lg, hptr(*), elptr(*), elvar(*), elst(*), + iss(*), ip, su(*), nsu, nqa, ftstem(*), + beptr(*), bep(*), kdist(*) double precision fuval(*), diag(*), dpr(*), w1(*), w2(*) * Internal variables integer i, j, jj, jp, ik, k, kk, iel, deb, intdim, + eldim, kdis, els logical elint double precision 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 * * We initialize the diagonal values to zero. * call dsetvl( arcs, diag, 1, zero ) * * We loop on the elements functions. * do 10 iel = 1 , elem * * We get some informations about * the element function IEL. * els = mod(elst(iel),8) if( els.eq.7 ) then deb = arcs + lg else deb = hptr(iel) endif eldim = elptr(iel+1) - elptr(iel) i = 0 * * We loop on the variables that belong * to the element function IEL. * do 20 j = elptr(iel), elptr(iel+1)-1 i = i + 1 jj = elvar(j) * * We test if the variable JJ is included * in the current independent set IP. * if( iss(jj).eq.ip ) then * * The variable JJ is included * in independent set IP. * if( els.eq.7 ) then * * Diagonal estimation of the Hessian matrix * when Dembo's method is applied. * diag(jj) = diag(jj) + dpr(deb+j) * else if( elint(elst(iel)) ) then * * We estimate the JJ-th component of * the diagonal of the Hessian matrix. * The element Hessian matrix has an * internal representation. * call range( iel, eldim, intdim, w1, w2, 0 ) call dsetvl( eldim, w1, 1, zero ) w1(i) = one call range( iel, eldim, intdim, w1, w2, 1 ) w1(i) = zero call ltcmvp( intdim, fuval(deb), w2, w1,.true. ) call range( iel, eldim, intdim, w1, w2, 2 ) diag(jj) = diag(jj) + w2(i) * else * * We estimate the JJ-th component of * the diagonal of the Hessian matrix. * The element Hessian matrix has an * internal representation. * if( eldim.eq.1 .or. i.eq.1 ) then diag(jj) = diag(jj) + fuval(deb) else jp = (i-1)*eldim - (( i**2 - 3*i + 2 ) / 2 ) diag(jj) = diag(jj) + fuval(deb+jp) endif endif endif 20 continue 10 continue * * We obtain the diagonal preconditioner * from the diagonal of the Hessian matrix. * sqrem = sqrt(epsmch) do 30 ik = 1 , nsu-nqa k = su(ik) kdis = kdist(ik) dpr(ik) = diag(k) * * We go over the flow augmenting path * of the superbasic variable K. * do 40 j = beptr(k) , beptr(k)+kdis-1 kk = abs(bep(j)) dpr(ik) = dpr(ik) + diag(kk) 40 continue * * The components of the diagonal preconditioner * must be positive and not too small. Otherwise, * they are set to one. * dpr(ik) = abs(dpr(ik)) if( dpr(ik).le.sqrem ) dpr(ik) = one 30 continue * return end