      subroutine ivmhess(n,m,x,y,arow,acol,aval,anze,b,epsdsn,
     +                   epsden,epsnlsn,maxoutit,maxinnit,maxfc,p,
     +                   extrap,hess,output,kprint,eta,theta,tau,tol,
     +                   lmin,lmax,gamma,sterel,steabs,epsmac,infty,
     +                   tau1,tau2,f,dsupn,deucn2,amax,outit,innit,fcnt,
     +                   gcnt,hcnt,modh,flag)

C     Subroutine IVM implements the Inexact Variable Metric method
C     to find the local minimizers of a given function with convex 
C     constraints (linear inequality constraints in this implementation), 
C     described in:
C
C     R. Andreani, E. G. Birgin, J. M. Martínez and J. Yuan, "Spectral 
C     projected gradient and variable metric methods for optimization 
C     with linear inequalities", IMA Journal of Numerical Analysis 25, 
C     pp. 221-252, 2005.
C
C     The user must supply the external subroutines evalf, evalg, 
C     and evalH to evaluate the objective function, its gradient and
C     its hessian, respectively. Depending on parameter hess, evalH
C     subroutine may be an empty subroutine.
C
C     This version 22 SET 2003 by E.G.Birgin, J.M.Martinez and J.Yuan.
C     Final revision ... by E.G.Birgin, J.M.Martinez and J. Yuan.
C
C     ON ENTRY:
C
C     n     integer,
C           number of variables,
C
C     m     integer,
C           number of linear inequality constraints,
C
C     x     double precision x(n),
C           initial estimation of the solution,
C
C     y     double precision y(m),
C           initial estimation of the Lagrange multipliers,
C
C     lda   integer
C           leading dimension of matrix A (see below),
C
C     A     double precision A(m,n),
C           linear inequalities coefficients matrix,
C
C     b     double precision b(m),
C           linear inequalities right hand side,
C
C     epsdsn double precision,
C           tolerance for the stopping criterion || amax d ||_sup < epsdsn,
C           where d is the IVM search direction and amax is
C           the maximum step which preserves feasibility,
C
C     epsden double precision,
C           tolerance for the stopping criterion || amax d ||_2 < epsden,
C           where d is the IVM search direction and amax is
C           the maximum step which preserves feasibility,
C
C     epsnlsn double precision,
C           tolerance for the stopping criterion of the inner
C           solver based on the sup-norm of the projected gradient
C           of the Lagrangian function,
C     
C     maxoutit integer,
C           maximum number of outer iterations,
C
C     maxinnit integer,
C           maximum number of inner iterations,
C
C     maxfc integer,
C           maximum number of function evaluations,
C
C     p     integer,
C           number of previous function values to be considered 
C           in the nonmonotone line search,
C
C           RECOMMENDED = 10
C
C     extrap logical,
C           if extrap = TRUE then, when solving the subproblem (using
C           an active set strategy), extrapolations are tried to increase
C           the dimension of the active set (it may be recomended for high
C           dimensional subproblems which comes from original problems with
C           many constraints); else extrapolations are not used at all.
C
C           RECOMMENDED = .true.
C
C     hess  logical,
C           if hess = TRUE then the true Hessian is used in the
C           quadratic subproblems (the Identity is used if the function
C           is linear); else the spectral approximation is used. 
C
C           RECOMMENDED = try both!
C
C     output integer,
C           < 0 means no print at all;
C           = 0 just initial and final information;
C           = 1 information at each outer iteration;
C           = 2 initial and final information for the inner
C               jobs: inner subproblems and line searches;
C           = 3 information of the inner iterations and 
C               the line search iterations;
C           = 4 detailed information of the inner iterations.
C
C     kprint integer,
C           maximum number of components to be printed when an
C           array must be showed,
C
C           RECOMMENDED = 5
C
C     eta   double precision,
C           eta in (0,1] controls the inexactness for solving the subproblem. 
C           eta = 1 means "solve the subproblem exactly".
C
C           RECOMMENDED: 0.9
C
C     theta double precision,
C           theta in (eta,1) is used to warranty that the iterates are
C           interior enough.
C
C           RECOMMENDED: 0.95
C
C     tau   double precision,
C           tau > 0 is used to compute the search direction of a 
C           subproblem of type q(y) = 0.5 y^t B y + b^t y as 
C           (B + sigma I) dk = - \nabla q(yk). See formulae (44) and
C           the end of Section 7 for details.
C
C           RECOMMENDED = 10^-6 
C
C     tol   double precision,
C           tol in (0,1) is used, i the algorithm to solve the subproblems
C           (which uses an active set strategy) to decide if the current 
C           face must be abandoned or not. Near 0 means to abandon the
C           face as fast as possible. Near 1 means try to find the solution
C           inside the current face.
C
C           RECOMMENDED: 0.9  
C
C     lmin  double precision,
C           lmin > 0, lower bound for the safeguarded spectral 
C           steplenght,
C
C           RECOMMENDED: 10^-3
C
C     lmax  double precision,
C           lmax >= lmin, upper bound for the safeguarded spectral 
C           steplenght,
C
C           RECOMMENDED: 10^3
C
C     gamma double precision,
C           gamma > 0, constant for the Armijo-like stooping criterion 
C           of the line search,
C
C           RECOMMENDED: 10^-4
C
C     sterel double precision,
C           epsrel > 0, a relative small number used for finite difference
C           approximation of derivatives,
C
C           RECOMMENDED: 10^-7
C
C     steabs double precision,
C           epsabs > 0, an absolute small number used for finite difference
C           approximation of derivatives,
C
C           RECOMMENDED: 10^-10
C
C     epsmac double precision,
C           epsmac > 0, machine epsilon, i.e., the smallest positive value 
C           such that 1 + epsmac > 1.
C
C           RECOMMENDED: In IEEE double precision floating point
C           arithmetic, the machine epsilon is something near 10^-16.
C           But the recomendation is to fix its value running the   
C           following subroutine:         
C
C           subroutine mcheps(epsmac)
C           double precision epsmac
C           epsmac = 1.0d0
C      20   if ( ( 1.0d0 + epsmac ) .ne. 1.0d0 ) then
C               epsmac = epsmac / 2.0d0
C               go to 20
C           end if
C           epsmac = 2.0d0 * epsmac
C           end
C
C     infty double precision,
C           a big number,
C
C           RECOMMENDED = 10^20
C
C     The next two parameters are parameters for subroutine modchl
C     which computes the "Eskow-Schnabel" modified Cholesky factorization
C     of a given symetric matrix, taken from:
C
C     ALGORITHM 695, COLLECTED ALGORITHMS FROM ACM.
C     PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C     VOL. 17, NO. 3, SEPTEMBER, 1991, PP. 306-312.
C
C     tau1  double precision,
C           tau1 > 0, tolerance used for determining when to switch phase 2
C
C           RECOMMENDED: 10^-6 (approx. epsmac ^ 1/3, where epsmac
C           is the machine epsilon (approx. 10^-16).
C
C     tau2  double precision,
C           tolerance used for determining the maximum condition 
C           number of the final 2X2 submatrix.
C
C           RECOMMENDED: 10^-6 (approx. epsmac ^ 1/3, where epsmac
C           tau2 > 0, is the machine epsilon (approx. 10^-16).
C
C     ON RETURN:
C
C     x     double precision x(n),
C           final estimation of to the local minimizer,
C
C     y     double precision y(m),
C           final estimation of the Lagrange multipliers,
C
C     f     double precision,
C           function value at the approximation to the local
C           minimizer,
C
C     dsupn double precision,
C           || d ||_sup at the final iteration,
C           where d is the IVM search direction,
C
C     deucn2 double precision,
C           || d ||_2^2 at the final iteration,
C           where d is the IVM search direction,
C
C     amax double precision,
C           maximum step to preserve feasibility in the last
C           computed direction,
C
C     outit  integer,
C           number of outer iterations,
C
C     innit  integer,
C           number of inner iterations,
C
C     fcnt  integer,
C           number of function evaluations,
C
C     gcnt  integer,
C           number of gradient evaluations,
C
C     hcnt  integer,
C           number of Hessian evaluations,
C
C     flag  integer,
C           termination parameter:
C           0= convergence with sup-norm of the IVM search direction,
C           1= convergence with euclidian-norm of the IVM search direction,

C           6= too many outer iterations,
C           7= too many inner iterations,
C           8= too many functional evaluations,

C          -1= error in evalf subroutine,
C          -2= error in evalg subroutine,
C          -3= error in evalH subroutine,
C          -4= error in cholcol subroutine,
C          -5= error in sscol subroutine.

C     PARAMETERS
      integer mmax
      parameter (mmax = 3 650 000)
      integer nmax
      parameter (nmax =     4 000)
      integer pmax
      parameter (pmax =       100)

C     SCALAR ARGUMENTS
      logical hess,extrap
      integer anze,fcnt,flag,gcnt,hcnt,innit,kprint,m,maxfc,
     +        maxinnit,maxoutit,modh,n,outit,output,p
      double precision amax,dsupn,deucn2,epsdsn,epsden,epsnlsn,f,eta,
     +        theta,tau,tol,lmin,lmax,gamma,tau1,tau2,sterel,steabs,
     +        epsmac,infty

C     ARRAY ARGUMENTS
      integer arow(anze),acol(anze)
      double precision aval(anze),b(m),x(n),y(m)

C     LOCAL SCALARS
      integer i,j,k,inform,nind,nprint,mprint,adyind
      double precision fbest,fnew,gtd,redc,L,Q,nLcsupn,nLpsupn,num,
     +        den,alpha,atemp,fmax,adymax,adymaxi,maxelem,lambda,sts,
     +        sty,gtg,tsmall,Ltmp,Qtmp,Lnew,Qnew,adytmp,adynew,amaxtmp,
     +        amaxnew,redctmp,redcnew,nLpi,w1max,w2max,nLmax,dymax,
     +        atamax,Hsafemax,sigma

C     LOCAL ARRAYS
      integer ind(mmax),arini(mmax+1)
      double precision C(nmax,nmax),H(nmax,nmax),Hsafe(nmax,nmax),
     +        E(nmax),xnew(nmax),xtmp(nmax),xbest(nmax),g(nmax),
     +        gnew(nmax),gtmp(nmax),d(nmax),dnew(nmax),dtmp(nmax),
     +        ynew(mmax),ytmp(mmax),ybest(mmax),nL(mmax),nLnew(mmax),
     +        nLtmp(mmax),nLc(mmax),dy(mmax),lastfv(0:pmax-1),
     +        bbar(mmax),w1(nmax),w2(nmax),w3(nmax),w4(nmax),w5(mmax)

C     EXTERNAL SUBROUTINES
      external cholcol,sscol,evallnlhess

C     INTRINSIC FUNCTIONS
      intrinsic abs,max,min,mod,sqrt

C     ===============================================================
C     INEXACT VARIABLE METRIC METHOD
C     ===============================================================

C     ===============================================================
C     INITALIZATION
C     ===============================================================

C     PRINT PRESENTATION INFORMATION

      if ( output .ge. 0 ) then
          write(*, fmt=9010)
          write(10,fmt=9010)
      end if

C     CONSTANTS

      nprint = min(kprint,n)
      mprint = min(kprint,m)

C     COUNTERS

      outit = 0
      innit = 0
      fcnt  = 0
      gcnt  = 0
      hcnt  = 0
      modh  = 0

C     SET WHERE EACH LINE OF SPARSE MATRIX A STARTS

      do i = 1,m
          arini(i) = 0
      end do
      arini(m+1) = anze + 1

      do i = 1,anze
          if ( arini(arow(i)) .eq. 0 ) then
              arini(arow(i)) = i
          end if
      end do

C     INITIALIZE LAST FUNCTIONAL VALUES VECTOR

      do i = 0,p - 1
          lastfv(i) = - infty
      end do

C     EVALUATE THE FUNCTION

      call evalf(n,x,f,inform)
      fcnt = fcnt + 1

      if ( inform .ne. 0 ) then
        ! ERROR IN EVALF SUBROUTINE, STOP
          flag = -1
          go to 900
      end if

C     EVALUATE THE GRADIENT

      call evalg(n,x,g,inform)
      gcnt = gcnt + 1

      if ( inform .ne. 0 ) then
        ! ERROR IN EVALG SUBROUTINE, STOP
          flag = -2
          go to 900
      end if

C     STORE FUNCTION VALUE FOR THE NONMONOTONE LINE SEARCH

      lastfv(0) = f

C     INITIALIZE BEST SOLUTION, BEST FUNCTIONAL VALUE AND ITS
C     ASSOCIATED LAGRANGE MULTIPLIERS

      fbest = f

      do i = 1,n
          xbest(i) = x(i)
      end do

      do i = 1,m
          ybest(i) = y(i)
      end do

C     DEFINE INITIAL SPECTRAL STEPLENGHT

      tsmall = 0.0d0
      do i = 1,n
          tsmall = max( tsmall, abs( x(i) ) )
      end do

      tsmall = max( sterel * tsmall, steabs )          

      do i = 1,n
          xtmp(i) = x(i) - tsmall * g(i)
      end do

      call evalg(n,xtmp,gtmp,inform)
      gcnt = gcnt + 1

      gtg = 0.0d0
      do i = 1,n
          gtg = gtg + g(i) ** 2
      end do

      sts = tsmall ** 2 * gtg

      sty = 0.0d0
      do i = 1,n
          sty = sty + g(i) * ( gtmp(i) - g(i) )
      end do
      sty = - tsmall * sty

      if ( sty .le. 0.0d0 ) then
          lambda = lmax

      else
          lambda = min( lmax, max( lmin, sts / sty ) )
      end if

C     ===============================================================
C     MAIN OUTER LOOP
C     ===============================================================

 100  continue

      outit = outit + 1

C     ===============================================================
C     STEP 1: COMPUTE THE SEARCH DIRECTION
C     ===============================================================

C     ---------------------------------------------------------------
C     Inner subproblem
C     ---------------------------------------------------------------

C     To compute the search direction we want to 
C
C     Minimize (approx) Qk(d) = 0.5 d^t Hk d + gk^t d 
C     subject to A(xk + d) <= b.

C     Equivalently (via duality), we want to 
C
C     Minimize (approx) -Lk(y) = 0.5 (gk + A^t y) Hk^-1 (gk + A^t y) + 
C                                y^t (b - A xk)
C     subject to y >= 0.

C     The primal solution is then obtained by: 
C
C     d = - Hk^-1 (gk + A^t y)

C     ---------------------------------------------------------------
C     Initialization
C     ---------------------------------------------------------------

      if ( output .ge. 2 ) then
          write(*, fmt=9200)
          write(10,fmt=9200)
      end if

C     Compute bbar = b - A xk

      do i = 1,m
          bbar(i) = b(i)
      end do

      do i = 1,anze
          bbar(arow(i)) = bbar(arow(i)) - aval(i) * x(acol(i))
      end do

c     Just to avoid bbar(i) = -10^-16 due to rounding errors
      do i = 1,m
          bbar(i) = max( 0.0d0, bbar(i) )
      end do

C     Evaluate the Hessian,

      if ( hess ) then

          call evalH(n,x,nmax,H,inform)
          hcnt = hcnt + 1

          if ( inform .ne. 0 ) then
            ! Error in evalh subroutine, stop
              flag = -3
              go to 900
          end if

C         Verify if the function is a linear function

          maxelem = 0.0d0
          do j = 1,n
              do i = 1,n
                  maxelem = max( maxelem, abs( H(i,j) ) )
              end do
          end do

          if ( maxelem .le. epsmac ) then
              do j = 1,n 
                  do i = 1,n
                      H(i,j) = 0.0d0
                  end do
                  H(j,j) = 1.0d0
              end do

              if ( output .ge. 2 ) then
                  write(*, fmt=9201)
                  write(10,fmt=9201)
              end if
          else
              if ( output .ge. 2 ) then
                  write(*, fmt=9202)
                  write(10,fmt=9202)
              end if
          end if

C     or use the spectral approximation

      else

          do j = 1,n
              do i = 1,n
                  H(i,j) = 0.0d0
              end do
              H(j,j) = 1.0d0 / lambda
          end do

          if ( output .ge. 2 ) then
              write(*, fmt=9203) lambda
              write(10,fmt=9203) lambda
          end if

      end if

C     Save the Hessian

      do j = 1,n
          do i = 1,n
              Hsafe(i,j) = H(i,j)
          end do
      end do

C     Modified Cholesky factorization of Hk

      call cholcolms(nmax,n,H,epsmac,E,inform)
      if ( inform .ne. 0 ) then
          flag = -4
          go to 900
      end if

C     COUNT THE NUMBER OF TIMES THE HESSIAN IS ACTUALLY MODIFIED

      if ( E(n) .gt. 0.0d0 ) then
          modh = modh + 1
          if ( hess ) then
              if (output .ge. 2 ) then
                  write(*, 9204) E(n)
                  write(10,9204) E(n)
              end if
          end if
      end if

C     Add E to the saved Hessian

      do i = 1,n
         Hsafe(i,i) = Hsafe(i,i) + E(i)
      end do

C     ---------------------------------------------------------------
C     Given y, this subroutine computes:
C
C     (a) The corresponding primal point: d = - Hk^-1 (gk + A^t y);
C
C     (b) The maximum step to mantain feasibility of xk + d in the
C         primal feasible set A( xk + alpha d ) <= b;
C
C     (c) The reduction constant to warranty that xk + redc * d is
C         an (enough) interior point;
C
C     (d) Lk(d,y) and Qk(redc d);
C
C     (e) The gradient of -L as a function of y.
C     ---------------------------------------------------------------

      call evallnlhess(m,n,nmax,y,g,arow,acol,aval,anze,bbar,H,d,
     +                 amax,redc,L,Q,nL,w1,w5,w3,w4,theta,infty,epsmac,
     +                 inform)

      if ( inform .ne. 0 ) then
          flag = -4
          go to 900
      end if

C     ---------------------------------------------------------------
C     Inner main loop
C     ---------------------------------------------------------------

 200  continue

C     Compute chopped antigradient and projected antigradient norms; 
C     and store in nind the number of free variables and in array nind 
C     its identifiers.

      nind = 0
      nLcsupn = 0.0d0
      nLpsupn = 0.0d0
      do i = 1,m
          if ( y(i) .gt. 0.0d0 ) then
              nind = nind + 1
              ind(i) =   nind
              nLc(i) =   0.0d0
              nLpi   = - nL(i)
          else
              ind(i) =   0
              nLc(i) =   max( 0.0d0, - nL(i) )
              nLpi   =   max( 0.0d0, - nL(i) )
          end if
          nLcsupn = max( nLcsupn, abs( nLc(i) ) )
          nLpsupn = max( nLpsupn, abs( nLpi   ) )
      end do

C     Print inner iteration information

      if ( output .ge. 3 ) then
          if ( output .ge. 4 ) then
              write(*, fmt=9210) innit,mprint,m,(y(i),i=1,mprint)
              write(*, fmt=9220) nind,m
              write(*, fmt=9240) -L,nLcsupn,nLpsupn
              write(*, fmt=9250) nprint,n,(d(i),i=1,nprint)
              write(*, fmt=9260) amax,redc,Q

              write(10,fmt=9210) innit,mprint,m,(y(i),i=1,mprint)
              write(10,fmt=9220) nind,m
              write(10,fmt=9240) -L,nLcsupn,nLpsupn
              write(10,fmt=9250) nprint,n,(d(i),i=1,nprint)
              write(10,fmt=9260) amax,redc,Q
          else
              write(*, fmt=9205) innit,-L,nLpsupn,nind
              write(10,fmt=9205) innit,-L,nLpsupn,nind
          end if
      end if

C     ---------------------------------------------------------------
C     Inner stooping criteria
C     ---------------------------------------------------------------

C     Testing precognized inner stooping criterion: 
C     Qk(redc d) <= eta L(d,y)

      if ( Q .le. eta * L ) then
        ! Precongnized stooping criterion satisfied, 
        ! stop the inner algorithm
          if ( output .ge. 2 ) then
              write(*, fmt=9270)
              write(10,fmt=9270)
          end if

          go to 300

      end if

C     Testing projected-gradient norm inner stooping criterion

      if ( nLpsupn .le. epsnlsn ) then
        ! Stationary point of the subproblem, 
        ! stop the inner algorithm
          if ( output .ge. 2 ) then
              write(*, fmt=9280)
              write(10,fmt=9280)
          end if
 
          go to 300
 
      end if

C     Testing number of inner iterations

      if ( innit .ge. maxinnit ) then
        ! Maximum number of inner iterations exceeded, 
        ! stop the algorithm
          flag = 7
          go to 900
      end if

C     ---------------------------------------------------------------
C     Doing an inner iteration
C     ---------------------------------------------------------------

      innit = innit + 1

C     ---------------------------------------------------------------
C     Testing whether we leave the current face or not and compute
C     the corresponding direction
C     ---------------------------------------------------------------

      if ( nLpsupn * tol .le. nLcsupn ) then

C         -----------------------------------------------------------
C         Leaving the face in the chopped gradient direcion
C         -----------------------------------------------------------

          if ( output .ge. 4 ) then
              write(*, fmt=9290)
              write(10,fmt=9290)
          end if

          do i = 1,m
              dy(i) = nLc(i)
          end do

      else

C         -----------------------------------------------------------
C         Interior to the face iteration (all the operations of this 
C         job are done into the reduced space)
C         -----------------------------------------------------------

          if ( output .ge. 4 ) then
              write(*, fmt=9300)
              write(10,fmt=9300)
          end if

C         Compute direction dy = - (A Hk^-1 A^t + sigma I)^-1 nL
C         in two steps:
C
C         First, solving (sigma Hk + A^t A) r = - A^t nL; and
C         second dy = - ( 1 / sigma ) ( nL + A r ).

C         Compute C = sigma Hk + A^t A

C         Compute A^t A

          do j = 1,n
              do i = j,n
                  C(i,j) = 0.0d0
              end do
          end do

          do k = 1,m
              if ( ind(k) .ne. 0 ) then
                  do j = arini(k),arini(k+1)-1
                      do i = j,arini(k+1)-1
                          C(acol(i),acol(j)) = 
     +                    C(acol(i),acol(j)) + aval(i) * aval(j) 
                      end do 
                  end do
              end if
          end do

          atamax = 0.0d0
          do j = 1,n
              do i = j,n
                  atamax = max( atamax, abs( C(i,j) ) )
              end do
          end do

C         Compute sigma

          Hsafemax = 0.0d0
          do j = 1,n
              do i = j,n
                  Hsafemax = max( Hsafemax, abs( Hsafe(i,j) ) )
              end do
          end do
          
          sigma = tau * atamax / Hsafemax

C         Add sigma Hk

          do j = 1,n
              do i = j,n
                  C(i,j) = C(i,j) + sigma * Hsafe(i,j)
                  C(j,i) = C(i,j)
              end do
          end do

C         Cholesky factorization of C

          call cholcols(nmax,n,C,epsmac,inform)
          if ( inform .ne. 0 ) then
              flag = -4
              go to 900
          end if

C         Solve system C w2 = - A^t nL

C         Compute w1 = - A^t nL

          do j = 1,n
              w1(j) = 0.0d0
          end do

          do i = 1,anze
              if ( ind(arow(i)) .ne. 0 ) then
                  w1(acol(i)) = w1(acol(i)) - aval(i) * nL(arow(i))
              end if
          end do
  
C         Solve system C w2 = w1

          call sscols(nmax,n,C,w1,w2,epsmac,inform)
          if ( inform .ne. 0 ) then
              flag = -5
              go to 900
          end if

C         Compute dy = - ( 1 / sigma ) ( nL + A w2 )

          do i = 1,m
              if (ind(i) .ne. 0 ) then
                  dy(i) = nL(i)
              else
                  dy(i) = 0.0d0
              end if
          end do

          do i = 1,anze
              if ( ind(arow(i)) .ne. 0 ) then
                  dy(arow(i)) = dy(arow(i)) + aval(i) * w2(acol(i))
              end if
          end do

          do i = 1,m
              dy(i) = - dy(i) / sigma
          end do

      end if

      if ( output .ge. 4 ) then
          write(*, fmt=9310) mprint,(dy(i),i=1,mprint) 
          write(10,fmt=9310) mprint,(dy(i),i=1,mprint) 
      end if

C     ---------------------------------------------------------------
C     Compute maximum step
C     ---------------------------------------------------------------

      adymax = infty
      adyind = 0
      do i = 1,m
          if ( dy(i) .lt. 0.0d0 ) then
              adymaxi = - y(i) / dy(i)
              if ( adymaxi .le. adymax ) then
                  adymax = adymaxi
                  adyind = i
              end if
          end if
      end do 

      if ( output .ge. 4 ) then
          if ( adyind .ne. 0 ) then
              write(*, fmt=9320) adymax,adyind,y(adyind)
              write(10,fmt=9320) adymax,adyind,y(adyind)
          else
              write(*, fmt=9321) adymax
              write(10,fmt=9321) adymax
          end if
      end if

C     ---------------------------------------------------------------
C     Compute the one-dimensional minimizer of the quadratic -Lk(d,y) 
C     along the dy direction (as the primal is feasible the dual is 
C     bounded and, so, this minimizer exists).
C     ---------------------------------------------------------------

C     alphady = - nL^t dy / ( dy^t A Hk^-1 A^t dy )

C     Compute w2 = Hk^-1 A^t dy

C     w1 = A^t dy

      do j = 1,n
          w1(j) = 0.0d0
      end do

      do i = 1,anze
          w1(acol(i)) = w1(acol(i)) + aval(i) * dy(arow(i))
      end do

C     Solve system Hk w2 = w1

      call sscols(nmax,n,H,w1,w2,epsmac,flag)
      if ( flag .ne. 0 ) then
          flag = -4
          go to 900
      end if

C     Compute the one-dimensional minimizer

C     num = nL^t dy
 
      nLmax = 0.0d0
      dymax = 0.0d0
      do i = 1,m
          nLmax = max( nLmax, abs( nL(i) ) )
          dymax = max( dymax, abs( dy(i) ) )
      end do

      if ( nLmax .eq. 0.0d0 .or. dymax .eq. 0.0d0 ) then
          num = 0.0d0
      else
          num = 0.0d0
          do i = 1,m
              num = num + ( nL(i) / nLmax ) * ( dy(i) / dymax )
          end do
          num = num * nLmax * dymax
      end if

C     den = v^t r

      w1max = 0.0d0
      w2max = 0.0d0
      do i = 1,n
          w1max = max( w1max, abs( w1(i) ) )
          w2max = max( w2max, abs( w2(i) ) )
      end do

      if ( w1max .eq. 0.0d0 .or. w2max .eq. 0.0d0 ) then
          den = 0.0d0
      else
          den = 0.0d0
          do i = 1,n
              den = den + ( w1(i) / w1max ) * ( w2(i) / w2max )
          end do
          den = den * w1max * w2max
      end if

C     One-dimensional minimizer

      adynew = - num / den

      if ( output .ge. 4 ) then
          write(*, fmt=9340) adynew        
          write(10,fmt=9340) adynew        
c         write(*, fmt=9341) adynew,num,den
c         write(10,fmt=9341) adynew,num,den
      end if

      adynew = min( adynew, adymax )

C     ---------------------------------------------------------------
C     Compute the first candidate to be the new iterate (without 
C     considering the extrapolation)
C     ---------------------------------------------------------------

      do i = 1,m
          ynew(i) = y(i) + adynew * dy(i)
      end do

      if ( adynew .eq. adymax .and. adyind .ne. 0 ) then
          ynew(adyind) = 0.0d0
      end if

      do i = 1,m
          if ( ynew(i) .le. epsmac ) then
              ynew(i) = 0.0d0
          end if
      end do

C     ---------------------------------------------------------------
C     Given y, this subroutine computes:
C
C     (a) The corresponding primal point: d = - Hk^-1 (gk + A^t y);
C
C     (b) The maximum step to mantain feasibility of xk + d in the
C         primal feasible set A( xk + alpha d ) <= b;
C
C     (c) The reduction constant to warranty that xk + redc * d is
C         an (enough) interior point;
C
C     (d) Lk(d,y) and Qk(redc d);
C
C     (e) The gradient of -L as a function of y.
C     ---------------------------------------------------------------

      call evallnlhess(m,n,nmax,ynew,g,arow,acol,aval,anze,bbar,H,
     +                 dnew,amaxnew,redcnew,Lnew,Qnew,nLnew,w1,w5,w3,
     +                 w4,theta,infty,epsmac,inform)

      if ( inform .ne. 0 ) then
          flag = -4
          go to 900
      end if

C     ---------------------------------------------------------------
C     Extrapolation
C     ---------------------------------------------------------------

      if ( extrap ) then

 250      continue

          adytmp = 2.0d0 * adynew

          do i = 1,m
              ytmp(i) = max( 0.0d0, y(i) + adytmp * dy(i) )
          end do

          call evallnlhess(m,n,nmax,ytmp,g,arow,acol,aval,anze,bbar,
     +                     H,dtmp,amaxtmp,redctmp,Ltmp,Qtmp,nLtmp,w1,
     +                     w5,w3,w4,theta,infty,epsmac,inform)

          if ( inform .ne. 0 ) then
              flag = -4
              go to 900
          end if

          if ( -Ltmp .le. -Lnew ) then

              adynew = adytmp

              do i = 1,m
                  ynew(i) = ytmp(i)
                  nLnew(i) = nLtmp(i)
              end do

              do i = 1,n
                  dnew(i) = dtmp(i)
              end do

              amaxnew = amaxtmp
              redcnew = redctmp

              Lnew = Ltmp
              Qnew = Qtmp

              go to 250

          end if

      end if

C     ---------------------------------------------------------------
C     Set next iterate
C     ---------------------------------------------------------------

      do i = 1,m
          y(i) = ynew(i)
          nL(i) = nLnew(i)
      end do

      do i = 1,n
          d(i) = dnew(i)
      end do

      amax = amaxnew
      redc = redcnew

      L = Lnew
      Q = Qnew

C     ---------------------------------------------------------------
C     End of the inner iteration
C     ---------------------------------------------------------------

      go to 200

C     ---------------------------------------------------------------
C     End of the inner main loop
C     ---------------------------------------------------------------

C     ===============================================================
C     PRINT OUTER ITERATION INFORMATION AND TEST STOPPING CRITERIA
C     ===============================================================

 300  continue

C     ---------------------------------------------------------------
C     Compute d norms
C     ---------------------------------------------------------------

      dsupn  = 0.0d0
      deucn2 = 0.0d0
      do i = 1,n
          dsupn  = max( dsupn, abs( d(i) ) )
          deucn2 = deucn2 + d(i) ** 2
      end do

C     ---------------------------------------------------------------
C     Print outer iteration information
C     ---------------------------------------------------------------

      if ( output .ge. 1 ) then
          write(*, fmt=9020) outit,f,dsupn,amax
          write(10,fmt=9020) outit,f,dsupn,amax
      end if

C     ---------------------------------------------------------------
C     Stooping criterion
C     ---------------------------------------------------------------

      if ( abs( redc ) * dsupn .le. epsdsn ) then
        ! IVM SEARCH DIRECTION SUP-NORM STOPPING CRITERION SATISFIED, STOP
          flag = 0
          go to 900
      end if

      if ( redc ** 2 * deucn2 .le. epsden ** 2 ) then
        ! IVM SEARCH DIRECTION 2-NORM STOPPING CRITERION SATISFIED, STOP
          flag = 1
          go to 900
      end if

      if ( outit .gt. maxoutit ) then
        ! MAXIMUM NUMBER OF OUTER ITERATIONS EXCEEDED, STOP
          flag = 6
          go to 900
      end if

C     ===============================================================
C     STEP 2: NONMONOTONE LINE SEARCH
C     ===============================================================

C     ---------------------------------------------------------------
C     Line search initialization
C     ---------------------------------------------------------------

      if ( output .ge. 2 ) then
          write(*, fmt=9500)
          write(10,fmt=9500)
      end if

C     Compute g^t d

      gtd = 0.0d0
      do i = 1,n
          gtd    = gtd + g(i) * d(i)
      end do

C     Compute the maximum functional value of the last p iterations

      fmax = lastfv(0)
      do i = 1,p - 1
          fmax = max(fmax,lastfv(i))
      end do

C     Compute first trial

      alpha = redc

      do i = 1,n
          xnew(i) = x(i) + alpha * d(i)
      end do

      call evalf(n,xnew,fnew,inform)
      fcnt = fcnt + 1

      if (inform .ne. 0) then
        ! Error in evalf subroutine, stop
          flag = -1
          go to 900
      end if

C     ---------------------------------------------------------------
C     Line search main loop
C     ---------------------------------------------------------------

 400  continue

C     Print iteration information

      if ( output .ge. 3 ) then
          write(*, fmt=9510) alpha,fnew
          write(10,fmt=9510) alpha,fnew
      end if

C     ---------------------------------------------------------------
C     Line search stooping criteria
C     ---------------------------------------------------------------

      if ( fnew .le. fmax + gamma * alpha * gtd ) then
        ! Nonmonotone Armijo-like stopping criterion satisfied, 
        ! stop the line search

          if ( output .ge. 2 ) then
              write(*, fmt=9520)
              write(10,fmt=9520)
          end if

          go to 500
      end if

      if ( fcnt .ge. maxfc ) then
        ! Maximum number of functional evaluations exceeded,
        ! stop the algorithm
          flag = 8
          go to 900
      end if 

C     ---------------------------------------------------------------
C     Line search iteration
C     ---------------------------------------------------------------

C     Safeguarded quadratic interpolation

      if ( alpha .le. 0.1d0 ) then
          alpha = 0.5d0 * alpha

      else
          atemp = ( - gtd * alpha ** 2 ) / 
     +            ( 2.0d0 * ( fnew - f - alpha * gtd ) )
          if ( atemp .lt. 0.1d0 .or. atemp .gt. 0.9d0 * alpha ) then
              atemp = 0.5d0 * alpha
          end if
          alpha = atemp
      end if

C     Compute trial point

      do i = 1,n
          xnew(i) = x(i) + alpha * d(i)
      end do

C     Evaluate the function

      call evalf(n,xnew,fnew,inform)
      fcnt = fcnt + 1

      if (inform .ne. 0) then
        ! Error in evalf subroutine, stop
          flag = -1
          go to 900
      end if

C     ---------------------------------------------------------------
C     End of the line search iteration
C     ---------------------------------------------------------------

      go to 400

C     ---------------------------------------------------------------
C     End of the line search main loop
C     ---------------------------------------------------------------

C     ===============================================================
C     PREPARING FOR THE NEXT OUTER ITERATION
C     ===============================================================

 500  continue

C     EVALUATE THE GRADIENT AT THE NEW ITERATE

      call evalg(n,xnew,gnew,inform)
      gcnt = gcnt + 1

      if ( inform .ne. 0 ) then
        ! ERROR IN EVALG SUBROUTINE, STOP
          flag = -2
          go to 900
      end if

C     COMPUTE THE SPECTRAL STEPLENGHT

      sts = 0.0d0
      sty = 0.0d0
      do i = 1,n
          sts = sts + ( xnew(i) - x(i) ) ** 2
          sty = sty + ( xnew(i) - x(i) ) * ( gnew(i) - g(i) )
      end do

      if ( sty .le. 0.0d0 ) then
          lambda = lmax

      else
          lambda = min( lmax, max( lmin, sts / sty ) )
      end if

C     SET THE NEW CURRENT POINT AND ITS GRADIENT

      do i = 1,n
          x(i) = xnew(i)
          g(i) = gnew(i)
      end do

C     SET NEW FUNCTION VALUE AND SAVE IT FOR THE NONMONOTONE
C     LINE SEARCH

      f = fnew
      lastfv(mod(outit,p)) = f

C     COMPARE THE NEW FUNCTION VALUE AGAINST THE BEST FUNCTION 
C     VALUE  AND, IF SMALLER, UPDATE THE BEST FUNCTION 
C     VALUE, THE CORRESPONDING BEST POINT AND ITS ASSOCIATED
C     LAGRANGE MULTIPLIERS

      if ( f .lt. fbest ) then
          fbest = f
          do i = 1,n
              xbest(i) = x(i)
          end do
          do i = 1,m
              ybest(i) = y(i)
          end do
      end if

C     ===============================================================
C     END OF THE OUTER ITERATION
C     ===============================================================

      go to 100

C     ===============================================================
C     END OF THE MAIN LOOP
C     ===============================================================

C     ===============================================================
C     RETURN
C     ===============================================================

 900  continue

C     SET X, Y AND F WITH THE BEST SOLUTION, ITS ASSOCIATED LAGRANGE
C     MULTIPLIERS AND ITS FUNCTIONAL VALUE

      f = fbest

      do i = 1,n
          x(i) = xbest(i)
      end do

      do i = 1,m
          y(i) = ybest(i)
      end do

C     WRITE OUTPUT STATUS

      if ( output .ge. 0 ) then

          if ( flag .eq. 0 ) then
              write(*, fmt=9600) flag,epsdsn
              write(10,fmt=9600) flag,epsdsn

          else if ( flag .eq. 1 ) then
              write(*, fmt=9610) flag,epsden
              write(10,fmt=9610) flag,epsden

          else if ( flag .eq. 6 ) then
              write(*, fmt=9620) flag,maxoutit
              write(10,fmt=9620) flag,maxoutit

          else if ( flag .eq. 7 ) then
              write(*, fmt=9630) flag,maxinnit
              write(10,fmt=9630) flag,maxinnit

          else if ( flag .eq. 8 ) then
              write(*, fmt=9640) flag,maxfc
              write(10,fmt=9640) flag,maxfc

          else if ( flag .eq. -1 ) then
              write(*, fmt=9650) flag
              write(10,fmt=9650) flag

          else if ( flag .eq. -2 ) then
              write(*, fmt=9660) flag
              write(10,fmt=9660) flag

          else if ( flag .eq. -3 ) then
              write(*, fmt=9670) flag
              write(10,fmt=9670) flag

          else if ( flag .eq. -4 ) then
              write(*, fmt=9680) flag
              write(10,fmt=9680) flag

          else if ( flag .eq. -5 ) then
              write(*, fmt=9690) flag
              write(10,fmt=9690) flag

          end if
      end if

C     WRITE STATISTICS

      if ( output .ge. 0 ) then
          write(*, fmt=9700) f,dsupn,sqrt(deucn2),amax,flag
          write(*, fmt=9710) outit,fcnt,gcnt,hcnt,modh,innit

          write(10,fmt=9700) f,dsupn,sqrt(deucn2),amax,flag
          write(10,fmt=9710) outit,fcnt,gcnt,hcnt,modh,innit
      end if

      return

C     NON-EXECUTABLE STATEMENTS

 9010 format(/,1X,'INEXACT VARIABLE METRIC METHOD')
 9020 format(/,1X,'OUTIT= ',I10,' F= ',1P,D14.7,' DSUPN= ',1P,D14.7,
     +            ' AMAX= ',1P,D14.7)

 9200 format(/,6X,'INNER SOLVER (COMPUTING THE SEARCH DIRECTION FOR ',
     +             'THE NEXT ITERATION)')
 9201 format(  6X,'(The function is linear, using Gk = I)')
 9202 format(  6X,'(Using Gk = the true Hessian)')
 9203 format(  6X,'(Using Gk = (1 / lambda) I, where lambda = ',
     +         1P,D14.7,' is the spectral',
     +       /,6X,'steplenght)')
 9204 format(  6X,'(',1P,D14.7,' was added to the diagonal)')
 9205 format(  6X,'INNIT = ',I6,' -L =',1P,D14.7,' NLPSUPN = ',
     +         1P,D14.7,' FACDIM = ',I6)
 9210 format(/,6X,'INNER ITERATION: ',I6,
     +       /,6X,'Current y (first ',I6,' components over a total of ',
     +         I6,'):',
     +       /,1(6X,6(1X,1PD11.4)))
 9220 format(  6X,'Number of free variables at this point: ',I6,
     +            ' (over a total of ',I6, ')')
 9230 format(  6X,'Indices of the free variables (first ',I6,'): ',
     +       /,1(6X,10(1X,I6)))
 9240 format(  6X,'Functional value                = ',1PD14.7,
     +       /,6X,'Chopped gradient sup-norm       = ',1PD14.7,
     +       /,6X,'Projected gradient sup-norm     = ',1PD14.7)
 9250 format(  6X,'Associated primal point (first ',I6,' components ',
     +            'over a total of ',I6,'):',
     +       /,1(6X,6(1X,1PD11.4)))
 9260 format(  6X,'Maximum primal step             = ',1PD14.7,
     +       /,6X,'Including constant of reduction = ',1PD14.7,
     +       /,6X,'Q(redc d)                       = ',1PD14.7)
 9270 format(  6X,'Precognized stooping criterion satisfied')
 9280 format(  6X,'Projected gradient norm stooping criterion ',
     +            'satisfied')
 9290 format(  6X,'Leaving the face')
 9300 format(  6X,'Remaining into the face')
 9310 format(  6X,'dy direction (first ',I6, ' components): ',
     +       /,1(6X,6(1X,1PD11.4)))
 9320 format(  6X,'Maximum step: ',1PD14.7,' Variable index: ',I6,
     +            ' value: ',1PD14.7)
 9321 format(  6X,'There is no maximum step. Maximum step: ',1PD14.7)
 9330 format(  6X,'The one-dimensional minimizer does not exist')
 9340 format(  6X,'One-dimensional minimizer: ',1PD14.7)
 9341 format(  6X,'One-dimensional minimizer: ',1PD14.7,1PD14.7,1PD14.7)

 9500 format(/,6X,'LINE SEARCH')
 9510 format(  6X,'Step = ',1P,D14.7,' F = ',1P,D14.7)
 9520 format(  6X,'Armijo-like stooping criterion satisfied')
 9600 format(/,1X,'IVM terminated with flag ',I2,': ',
     +            'convergence with the sup-norm of the search ',
     +       /,1X,'direction smaller than ',1P,D14.7,
     +            ' (parameter epsdsn)')
 9610 format(/,1X,'IVM terminated with flag ',I2,': ',
     +            'convergence with the euclidian-norm of the ',
     +            'search',
     +       /,1X,'direction smaller than ',1P,D14.7,
     +            ' (parameter epsden)')
 9620 format(/,1X,'IVM terminated with flag ',I2,': ',
     +            'too many outer iterations ',
     +       /,1X,'(more than parameter maxoutit = ',I6,')')
 9630 format(/,1X,'IVM terminated with flag ',I2,': ',
     +            'too many inner iterations ',
     +       /,1X,'(more than parameter maxinnit = ',I6,')')
 9640 format(/,1X,'IVM terminated with flag ',I2,': ',
     +            'too many functional evaluations ',
     +       /,1X,'(more than parameter maxfc = ',I6,')')
 9650 format(/,1X,'IVM terminated with flag ',I2,': ',
     +            'Error in evalf subroutine')
 9660 format(/,1X,'IVM terminated with flag ',I2,': ',
     +            'Error in evalg subroutine')
 9670 format(/,1X,'IVM terminated with flag ',I2,': ',
     +            'Error in evalH subroutine')
 9680 format(/,1X,'IVM terminated with flag ',I2,': ',
     +            'Error in cholcol subroutine')
 9690 format(/,1X,'IVM terminated with flag ',I2,': ',
     +            'Error in sscol subroutine')
 9700 format(/,1X,'F        = ',   1P,D17.10,
     +       /,1X,'DSUPNORM = ',1X,1P,D16.10,
     +       /,1X,'DEUCNORM = ',1X,1P,D16.10,
     +       /,1X,'AMAX     = ',1X,1P,D16.10,
     +       /,1X,'FLAG     = ',15X,I2)
 9710 format(/,1X,'OUTIT    = ',1X,I16,
     +       /,1X,'FCNT     = ',1X,I16,
     +       /,1X,'GCNT     = ',1X,I16,
     +       /,1X,'HCNT     = ',1X,I16,
     +       /,1X,'MODH     = ',1X,I16,
     +       /,1X,'INNIT    = ',1X,I16)

 9999 format(9(1X,1P,D12.4))

      end

c     ******************************************************
c     ******************************************************

      subroutine evallnlhess(m,n,ldh,y,g,arow,acol,aval,anze,bbar,H,d,
     +           amax,redc,L,Q,nL,w1,w2,w3,w4,theta,infty,epsmac,inform)

C     SCALAR ARGUMENTS
      integer anze,m,n,ldh,inform
      double precision amax,redc,L,Q,theta,infty,epsmac

C     ARRAY ARGUMENTS
      integer arow(anze),acol(anze)
      double precision y(m),g(n),aval(anze),H(ldh,n),d(n),
     +        bbar(m),nL(m),w1(n),w2(m),w3(n),w4(n)

C     LOCAL SCALARS
      integer i,j
      double precision tmp,dmax,gmax,w1max,ymax,w2mbbarmax

C     ---------------------------------------------------------------
C     Computing corresponding primal point d = - Hk^-1 ( gk + A^t y )
C     ---------------------------------------------------------------

C     w1 = - ( gk + A^t y )

      do j = 1,n
          w1(j) = - g(j)
      end do

      do i = 1,anze
          w1(acol(i)) = w1(acol(i)) - aval(i) * y(arow(i))
      end do

C     Solve system Hk d = w1

      call sscols(ldh,n,H,w1,d,epsmac,inform)
      if ( inform .ne. 0 ) then
          return
      end if

C     ---------------------------------------------------------------
C     Computing maximum step to mantain feasibility of xk + d
C     in the primal feasible set A( xk + alpha d ) <= b
C     ---------------------------------------------------------------

C     w2  = A d

      do i = 1,m
          w2(i) = 0.0d0
      end do

      do i = 1,anze
          w2(arow(i)) = w2(arow(i)) + aval(i) * d(acol(i))
      end do

C     amax = min_{i | [A d]_i > 0} { [bbar]_i / [A d]_i }

      amax = infty
      do i = 1,m
          if ( w2(i) .gt. 0.0d0 ) then
              amax = min( amax, bbar(i) / w2(i) )
          end if
      end do

C     ---------------------------------------------------------------
C     Computing the reduction constant to warranty that xk + redc * d 
C     is an (enough) interior point
C     ---------------------------------------------------------------

      redc = min( 1.0d0, theta * amax )

C     ---------------------------------------------------------------
C     Computing Lk(d,y) and Qk(redc d) 
C     ---------------------------------------------------------------

      dmax  = 0.0d0
      gmax  = 0.0d0
      w1max = 0.0d0
      do i = 1,n
          dmax  = max( dmax, abs( d(i) ) )
          gmax  = max( gmax, abs( d(i) ) )
          w1max = max( w1max, abs( w1(i) ) )
      end do

      ymax  = 0.0d0
      w2mbbarmax = 0.0d0
      do i = 1,m
          ymax  = max( ymax, abs( y(i) ) )
          w2mbbarmax = max( w2mbbarmax, abs( w2(i) - bbar(i) ) )
      end do

      dmax = 1.0d0
      gmax = 1.0d0
      w1max = 1.0d0
      ymax = 1.0d0
      w2mbbarmax = 1.0d0

C     Computing 0.5 d^t Hk d = - d^t (gk + A^t y) = d^t w1

      if ( dmax .eq. 0.0d0 .or. w1max .eq. 0.0d0 ) then
          tmp = 0.0d0
      else
          tmp = 0.0d0
          do i = 1,n
              tmp = tmp + ( d(i) / dmax ) * ( w1(i) / w1max )
          end do
          tmp = tmp * dmax * w1max
      end if

      L = 0.5d0 * tmp
      Q = 0.5d0 * redc ** 2 * tmp

C     Computing gk^t d

      if ( gmax .eq. 0.0d0 .or. dmax .eq. 0.0d0 ) then
          tmp = 0.0d0
      else
          tmp = 0.0d0
          do i = 1,n
              tmp = tmp + ( g(i) /gmax ) * ( d(i) / dmax )
          end do
          tmp = tmp * gmax * dmax
      end if

      L = L + tmp
      Q = Q + redc * tmp

C     Computing y^t ( A xk - b + A d  ) = y^t ( w2 - bbar )

      if ( ymax .eq. 0.0d0 .or. w2mbbarmax .eq. 0.0d0 ) then
          tmp = 0.0d0
      else
          tmp = 0.0d0
          do i = 1,m
          tmp = tmp + ( y(i) / ymax ) * ( w2(i) - bbar(i) ) / w2mbbarmax
          end do
          tmp = tmp * ymax * w2mbbarmax
      end if

      L = L + tmp

C     ---------------------------------------------------------------
C     Computing - \nabla L(d,y) = A Hk^-1 (gk + A^t y) + bbar = bbar - A d
C     ---------------------------------------------------------------

C     nL is the antigradient of L. As we are minimizing -L, nL is the
C     positive gradient of the function being minimized.

      do i = 1,m
          nL(i) = bbar(i) - w2(i)
      end do

      end


c     ******************************************************
c     ******************************************************

      subroutine cholcolms(lda,n,A,epsmac,E,inform)

C     SCALAR ARGUMENTS
      integer lda,n,inform
      double precision epsmac

C     ARRAY ARGUMENTS
      double precision A(lda,n),E(n)

C     LOCAL SCALARS
      integer i,k,base,lastb,npol
      double precision gamma

      inform = 0

C     Number of Hi submatrices

      npol = n / 3 - 1

      lastb = 3 * npol

C     Add gamma to the submatrices diagonals to make them positive definite

      gamma = 1.0d-4

      do k = 1,npol

          base = 3 * ( k - 1 )

          do i = 1,3
              E(base+i) = gamma
              A(base+i,base+i) = A(base+i,base+i) + E(base+i)
          end do

      end do

      do i = 1,3
          E(lastb+i) = ( npol + 1 ) * gamma
          A(lastb+i,lastb+i) = A(lastb+i,lastb+i) + E(lastb+i)
      end do

      call cholcols(lda,n,A,epsmac,inform)

      return

      end

c     ******************************************************
c     ******************************************************

      subroutine cholcols(lda,n,A,epsmac,inform)

C     SCALAR ARGUMENTS
      integer lda,n,inform
      double precision epsmac

C     ARRAY ARGUMENTS
      double precision A(lda,n)

C     LOCAL SCALARS
      integer i,j,k,l,base,lastb,npol

C     LOCAL ARRAYS
      double precision M(3,3),B(3,3),x(3)

      inform = 0

C     Number of Hi+Ai submatrices

      npol = n / 3 - 1

      lastb = 3 * npol

C     Compute the Cholesky decomposition of each submatrix.

      do k = 1,npol + 1

          base = 3 * ( k - 1 )

C         Copy Hi+Ai into M

          do j = 1,3
              do i = 1,3
                  M(i,j) = A(base+i,base+j)
              end do
          end do

C         Cholesky decomposition of M
      
          call cholcol(3,3,M,epsmac,inform)
          if ( inform .ne. 0 ) then
              return
          end if

C         Delete the elements above the diagonal

          do j = 1,3
              do i = 1,j-1
                  M(i,j) = 0.0d0
              end do
          end do

C         Mount the Cholesky decomposition of matrix A

          do j = 1,3
              do i = 1,3
                  A(base+i,base+j) = M(i,j)
              end do
          end do

          if ( k .ne. npol + 1 ) then

              do j = 1,3

                  do i = 1,3
                      x(i) = A(lastb+i,base+j)
                  end do

                  call sscol1(3,3,M,x,epsmac,inform)
                  if ( inform .ne. 0 ) then
                      return
                  end if

                  do i = 1,3
                      B(i,j) = x(i)
                  end do

              end do

              do j = 1,3
                  do i = 1,3
                      A(lastb+i,base+j) = B(j,i)
                  end do
              end do

              do i = 1,3
                  do j = 1,3
                      do l = 1,3
                          A(lastb+i,lastb+j) =
     +                    A(lastb+i,lastb+j) - B(l,i) * B(l,j)
                      end do
                  end do
              end do
     
          end if

      end do

      return

      end

c     ******************************************************
c     ******************************************************

      subroutine cholcol(lda,n,A,epsmac,inform)

      integer lda,n,inform
      double precision A(lda,n),epsmac

      integer i,j,k

      inform = 0

      do k = 1,n

          if ( A(k,k) .le. epsmac ) then
              inform = -1
              return
          end if

          A(k,k) = sqrt( A(k,k) )

          do i = k+1,n
              A(i,k) = A(i,k) / A(k,k)
          end do
 
          do j = k+1,n
              do i = j,n
                  A(i,j) = A(i,j) - A(i,k) * A(j,k)
              end do
          end do
      end do

      return

      end

c     ******************************************************
c     ******************************************************

      subroutine sscol(lda,n,A,b,x,epsmac,inform)

      integer lda,n,inform
      double precision A(lda,n),b(n),x(n),epsmac

      integer i

      inform = 0

      do i = 1,n
          x(i) = b(i)
      end do

C     Solving G y = b

      call sscol1(lda,n,A,x,epsmac,inform)
      if ( inform .ne. 0 ) then
          return
      end if

C     Solving G^t x = y

      call sscol2(lda,n,A,x,epsmac,inform)
      if ( inform .ne. 0 ) then
          return
      end if

      return

      end

c     ******************************************************
c     ******************************************************

      subroutine sscols(lda,n,A,b,x,epsmac,inform)

C     SCALAR ARGUMENTS
      integer lda,n,inform
      double precision epsmac

C     ARRAY ARGUMENTS
      double precision A(lda,n),b(n),x(n)

C     LOCAL SCALARS
      integer i,j,k,base,lastb,npol

C     LOCAL ARRAYS
      double precision M(3,3),y(3)

      inform = 0

      do i = 1,n
          x(i) = b(i)
      end do

      npol = n / 3 - 1

      lastb = 3 * npol

C     Solving G y = b

      do k = 1,npol + 1

          base = 3 * ( k - 1 )

          do j = 1,3
              do i = j,3
                  M(i,j) = A(base+i,base+j)
              end do
          end do

          do i = 1,3
              y(i) = x(base+i)
          end do

          call sscol1(3,3,M,y,epsmac,inform)
          if ( inform .ne. 0 ) then
              return
          end if

          do i = 1,3
              x(base+i) = y(i)
          end do

          if ( k .ne. npol + 1 ) then

              do j = 1,3
                  do i = 1,3
                      x(lastb+i) = x(lastb+i) - A(lastb+i,base+j) * y(j)
                  end do
              end do

          end if

      end do

C     Solving G^t x = y

      do k = npol + 1,1,-1

          base = 3 * ( k - 1 )

          do j = 1,3
              do i = j,3
                  M(i,j) = A(base+i,base+j)
              end do
          end do

          do i = 1,3
              y(i) = x(base+i)
          end do

          if ( k .ne. npol + 1 ) then

              do j = 1,3
                  do i = 1,3
                      y(i) = y(i) - A(lastb+j,base+i) * x(lastb+j)
                  end do
              end do

          end if

          call sscol2(3,3,M,y,epsmac,inform)
          if ( inform .ne. 0 ) then
              return
          end if

          do i = 1,3
              x(base+i) = y(i)
          end do

      end do

      return

      end

c     ******************************************************
c     ******************************************************

      subroutine sscol1(lda,n,A,x,epsmac,inform)

      integer lda,n,inform
      double precision A(lda,n),x(n),epsmac

      integer i,j

C     Solving G y = b

      do j = 1,n

          if ( A(j,j) .le. epsmac ) then
              inform = -1
              return
          end if

          x(j) = x(j) / A(j,j)

          do i = j+1,n
              x(i) = x(i) - x(j) * A(i,j)
          end do

      end do

      return

      end

c     ******************************************************
c     ******************************************************

      subroutine sscol2(lda,n,A,x,epsmac,inform)

      integer lda,n,inform
      double precision A(lda,n),x(n),epsmac

      integer i,j

C     Solving G^t x = y

      do i = n,1,-1

          do j = i+1,n
              x(i) = x(i) - A(j,i) * x(j)
          end do

          x(i) = x(i) / A(i,i)

          if ( A(i,i) .le. epsmac ) then
              inform = -1
              return
          end if

      end do

      return

      end

