* ************************************************************************* subroutine cg( itcg, itcgmx, rgra, dir, ds, gds, res, dw1, dw2, + x, fuval, gptr, lg, negval, hptr, fcalc, lfc, nqa, + nsu, su, nba, ba, elst, elvar, elptr, bep, beptr, + bgx, 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 truncated conjugate gradient method. 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. * 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. * dw2 ( dble ) * input : array used as workspace. * output : meaningless. * 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. * bgx ( dble ) * input : array with the gradient at x. * output : unmodified (only used for Dembo's strategy). * 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, max, min. * dsetvl, dnrm2, 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(*), rgra(*), dir(*), ds(*), gds(*), res(*), + dw1(*), dw2(*), fuval(*), bgx(*) * Internal variables integer ns, it, ik, k double precision etatn, cond, rayl, oresn, nresn, sqrem, + crmax, crmin, rgnorm, dsnor2, alpha, beta, + curv, dnrm2 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) ds(k) = -rgra(k) 10 continue * * We define the stopping criteria ETATN for * the truncated conjugate gradient. * rgnorm = dnrm2( ns, res, 1 ) oresn = rgnorm etatn = max( sqrem , min( tenm2 , rgnorm ) ) * * 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) = -res(ik) + beta * ds(k) 30 continue dsnor2 = oresn**2 + 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, dw2, elvar, + elptr, elst, ba, nba, su, ns, bep, beptr, + ftstem, kdist, bgx ) 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**2 / 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. 1.0d0 ) 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 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) 40 continue * * We update the residual norm. * nresn = dnrm2( ns, res, 1 ) * * We test if the stopping criteria is verified. * if( nresn/rgnorm .gt. etatn ) then * * The stopping criteria is not verified. * A new conjugate gradient iteration is * executed. * beta = (nresn/oresn)**2 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,nresn,cond 1002 format(/,' CG 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