* ************************************************************************* subroutine pcg( itcg, itcgmx, dpr, rgra, dir, ds, gds, res, dw1, + pres, x, fuval, gptr, lg, negval, hptr, fcalc, + lfc, nqa, nsu, su, nba, ba, elst, elvar, elptr, + bep, beptr, ftstem, kdist, inform, ipdevc, + what, tprt ) * ************************************************************************* * Purpose : * --------- * This routine solves the following linear system of equations * t t * Z G Z ds = - Z g(x) * using a preconditioned truncated conjugate gradient method. * The preconditioner used is a diagonal matrix. The solution * of this system is the superbasic stepdirection ds. * Parameters : * ------------ * itcg ( int ) * input : the total number of conjugate gradient iterations * executed before solving this linear system. * output : the total number of conjugate gradient iterations * executed after solving this linear system. * itcgmx ( int ) * input : the maximal number of conjugate gradient * iterations needed if all the linear system * were solved exactly before solving this * linear system. * We consider here that a NxN linear system * requires N conjugate gradient iterations to * be solved exactly. * This quantity is used to estimate the average * size of the linear systems solved during the * minimization. * output : the maximal number of conjugate gradient * iterations needed if all the linear system * were solved exactly after solving this * linear system. * dpr ( dble ) * input : vector containing the diagonal preconditioner. * output : unmodified. * rgra ( dble ) * input : the reduced gradient vector. * output : unmodified. * dir ( dble ) * input : meaningless. * output : the superbasic stepdirection computed by the * truncated conjugated gradient method. * ds ( dble ) * input : array used as workspace. * output : it contains the direction of the last * conjugate gradient iteration. * gds ( dble ) * input : array used as workspace. * output : it contains the product of the reduced Hessian * matrix times the direction of the last conjugate * gradient iteration. * res ( dble ) * input : array used as workspace. * output : it is the residual vector of the linear system. * dw1 ( dble ) * input : array used as workspace. * output : meaningless. * pres ( dble ) * input : array used as workspace. * output : it is the residual vector of the linear system * scaled by the diagonal preconditioner. * 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. * nqa ( int ) * input : the number of superbasic variables that are * quasi active. * output : unmodified. * nsu ( int ) * input : the number of superbasic variables. * output : unmodified. * su ( int ) * input : vector containing the indices of the superbasic * variables. * output : unmodified. * nba ( int ) * input : the number of basic variables. * output : unmodified. * ba ( int ) * input : vector containing the indices of the basic * variables. * output : unmodified. * elst ( int ) * input : vector containing the status of the element * functions. * 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. * elptr ( int ) * input : array whose kth value is the position of * the first variable of element k, in the * list ELVAR. * 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. * 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. * ipdevc ( int ) * input : output device unit number for printing messages. * output : unmodified. * what ( int ) * input : If >= 3 and TPRT = .true, then the conjugate * gradient iterations results are printed out. * Otherwise, no output is produced. * output : unmodified. * tprt ( log ) * input : .true. if output is generated this iteration. * .false. if no output is generated this iteration. * output : unmodified. * Routines used : * --------------- * sqrt, mqx, min. * dsetvl, ddot, second, prod, rg. * Programming : * ------------- * D. Tuyttens * ======================================================================== * Routine parameters integer itcg, itcgmx, nqa, nsu, su(*), nba, ba(*), + elst(*), elvar(*), elptr(*), bep(*), beptr(*), + ftstem(*), gptr(*), lg, negval, hptr(*), + inform, ipdevc, what, fcalc(*), lfc, kdist(*) logical tprt double precision x(*), dpr(*), rgra(*), dir(*), ds(*), gds(*), + res(*), pres(*), dw1(*), fuval(*) * Internal variables integer ns, it, ik, k double precision etatn, cond, rayl, oresn, nresn, sqrem, + crmax, crmin, rgnorm, dsnor2, alpha, beta, + curv, ddot real time1, time2 * Save variables save ns, it, ik, k save etatn, cond, rayl, oresn, nresn, sqrem, + crmax, crmin, rgnorm, dsnor2, alpha, beta, + curv save time1, time2 * Common specifications integer arcs, nodes, elem double precision epsmch, huge, tiny, tol real xocpu, totcpu, cgcpu, prdcpu, lnscpu, updcpu common / prbdim / arcs, nodes, elem common / prbmch / epsmch, huge, tiny, tol common / prbcpu / xocpu, totcpu, cgcpu, prdcpu, lnscpu, updcpu 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 200 * * Some initializations. * ns = nsu - nqa beta = zero crmin = huge crmax = -huge sqrem = sqrt(epsmch) * * We define the initial direction and * residual vectors. * call dsetvl( arcs, ds, 1, zero ) do 10 ik = 1 , ns k = su(ik) res(ik) = rgra(k) pres(ik) = res(ik) / dpr(ik) ds(k) = -pres(ik) 10 continue * * We define the stopping criteria ETATN for the * preconditioned truncated conjugate gradient. * oresn = ddot( ns, res, 1, pres, 1 ) rgnorm = sqrt(oresn) etatn = max( sqrem , min( tenm2 , rgnorm ) )**2 * * Preconditioned conjugate gradient iterations loop. * it = 1 20 continue * * We update the minimization direction DS * and its norm DSNOR2. * do 30 ik = 1 , ns k = su(ik) ds(k) = -pres(ik) + beta * ds(k) 30 continue dsnor2 = oresn + dsnor2 * beta**2 * * We compute the product of the reduced Hessian * matrix times the vector DS. * call second(time1) 200 continue call prod( x, fuval, gptr, lg, negval, hptr, fcalc, lfc, + inform, curv, ds, gds, dw1, pres, elvar, + elptr, elst, ba, nba, su, ns, bep, beptr, + ftstem, kdist, dpr ) if( inform.ne.0 ) return call second(time2) prdcpu = prdcpu + (time2-time1) * rayl = curv / dsnor2 * * We obtain informations about the conditioning of the problem * and the positive definiteness of the reduced Hessian matrix. * if( rayl .le. sqrem ) then * * The reduced Hessian matrix is not positive definite. * if( tprt .and. what.ge.3 ) then write(ipdevc,1000) 1000 format(/' Reduced Hessian is not positive definite !!') endif * * We stop the conjugate gradient iterations. * go to 99 else * * We update the estimate COND of the conditioning * of the problem. * alpha = oresn / curv crmax = max( crmax , rayl ) crmin = max( sqrem , min( crmin , rayl ) ) cond = crmax / crmin endif * * We test the conditioning of the problem. * if( sqrem*cond .gt. one ) then * * The problem is ill-conditioned. * if( tprt .and. what.ge.3 ) then write(ipdevc,1001) cond 1001 format(/' The problem is ill-conditioned !!! ',/, + ' Cond(Hessian) >>> ',d12.5) endif * * We stop the conjugate gradient iterations. * go to 99 endif * * The preconditioned conjugate gradient step. * do 40 ik = 1 , ns k = su(ik) dir(k) = dir(k) + alpha * ds(k) res(ik) = res(ik) + alpha * gds(k) pres(ik) = res(ik) / dpr(ik) 40 continue * * We update the scaled residual norm. * nresn = ddot( ns, res, 1, pres, 1 ) * * We test if the stopping criteria is verified. * if( nresn/rgnorm**2 .gt. etatn ) then * * The stopping criteria is not verified. * A new conjugate gradient iteration is * executed. * beta = nresn/oresn oresn = nresn else * * The stopping criteria is verified. We * stop the conjugate gradient iterations. * go to 100 endif * it = it + 1 if ( it .le. 3 * ns ) go to 20 * * The maximum number of conjugate gradient iterations * allowed is reached. * 99 continue * * When the first direction has a negative curvature, * the direction is set to minus the reduced gradient. * if( it.eq.1 ) call rg( nsu, su, dir, rgra ) * * We print out the conjugate gradient results. * 100 continue itcgmx = itcgmx + ns itcg = itcg + it * if( tprt .and. what.ge.3 ) then write(ipdevc,1002) it,ns,itcg,itcgmx,sqrt(nresn),cond 1002 format(/,' PCG direction : cg iterations ',4x,i5,' /',i5, + /,20x,'total',14x,i5,' /',i5,//,20x,'residual norm',6x, + d12.5,/,20x,'conditionning',6x,d12.5) endif * return end