      program ivmmahess

C     IVM Test driver. This master file uses IVM (Inexact Variable
C     Metric Method) to solve some the 3D-lOcation problem.
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     This version 25 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     PARAMETERS
      integer mmax
      parameter (mmax    = 3 650 000)
      integer nmax
      parameter (nmax    =     4 000)
      integer anzemax
      parameter (anzemax =  3 * mmax)

C     COMMON SCALARS
      integer anze,n

C     COMMON ARRAYS
      integer arow(anzemax),acol(anzemax)
      double precision aval(anzemax)

C     LOCAL SCALARS
      logical extrap,hess
      integer fcnt,flag,gcnt,hcnt,i,innit,j,kprint,m,maxfc,maxinnit,
     +        maxoutit,modh,outit,output,p,ntpmin,ntpmax,actconstr,
     +        inform
      double precision amax,dsupn,deucn2,epsdsn,epsden,epsmac,epsnlsn,
     +        eta,f,gamma,infty,lmin,lmax,steabs,sterel,tau,tau1,tau2,
     +        theta,tol,prob,rmin,rmax,bmin,fsol,diff
      character * 10 pname
      real time,dtime

C     LOCAL ARRAYS
      integer ngp(3)
      double precision b(mmax),l(nmax),u(nmax),x(nmax),y(mmax),step(3),
     +        sol(nmax)
      real dum(2)

C     EXTERNAL SUBROUTINES
      external ivmhess

C     INTRINSIC FUNCTIONS
      intrinsic dtime

C     COMMON BLOCKS
      common /ldata1/ arow,acol,aval,anze,n

C     SET PROBLEM DATA

      pname = 'LOCATI3D'

      ngp(1) = 11
      ngp(2) = 11
      ngp(3) = 11

      step(1) = 1.0d0
      step(2) = 1.0d0
      step(3) = 1.0d0

      prob = 0.01d0

      rmin = 1.0d0
      rmax = 1.0d0

      ntpmin = 200
      ntpmax = 200

      read(*,*) prob,ntpmin,ntpmax

c      prob   = 0.0015d0
c      ntpmin =    100
c      ntpmax =    100

      if ( prob .eq. 0.0d0 ) then
          stop
      end if

      call genpro(ngp,step,prob,rmin,rmax,ntpmin,ntpmax,m,n,arow,acol,
     +            aval,anze,b,l,u,x,sol)

      call evalf(n,sol,fsol,inform)
      write(*,*) 'There is a feasible point where f = ',fsol

C     INCLUDE THE BOUND CONSTRAINTS IN THE FORMAT A x <= b 

      if ( m + 2 * n .gt. mmax ) then
          write(*,*) 'Increase mmax and run again'
          stop
      end if

      do i = 1,n
          m = m + 1
          anze = anze + 1
          arow(anze) = m
          acol(anze) = i
          aval(anze) = - 1.0d0
          b(m) = - l(i)

          m = m + 1
          anze = anze + 1
          arow(anze) = m
          acol(anze) = i
          aval(anze) = 1.0d0
          b(m) = u(i)
      end do

      if ( anze .gt. anzemax ) then
          write(*,*) 'Increase anzemax and run again'
          stop
      end if

      write(*,*) 'PROBLEM DATA'
      write(*,*) 'Name : 3D Location Problem'
      write(*,*) 'Number of variables                    : ',n
      write(*,*) 'Number of constraints of type A x <= b : ',m
      write(*,*) 'Number of non-zero elements in matrix  : ',anze

C     SET THE INITIAL ESTIMATION OF THE LAGRANGE MULTIPLIERS

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

C     SETUP OPTIMIZER PARAMETERS 

      call iniivm(epsdsn,epsden,epsnlsn,maxoutit,maxinnit,maxfc,p,
     +            extrap,hess,output,kprint,eta,theta,tau,tol,
     +            lmin,lmax,gamma,sterel,steabs,epsmac,infty,tau1,
     +            tau2)

C     CALL THE OPTIMIZER

      time= dtime(dum)

      call 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)

      time= dtime(dum)

      time= dum(1)

C     WRITE STATISTICS

      write(*,*)
      write(*,*) 'Solution (f = ',f,' Known feasible point (f =', 
     +           fsol,')'
      diff = 0.0d0
      do i = 1,n/3
          write(*,fmt=9010) (x(3*(i-1)+j),j=1,3),(sol(3*(i-1)+j),j=1,3)
          do j = 1,3
              diff = max( diff, abs(x(3*(i-1)+j)-sol(3*(i-1)+j)) )
          end do
      end do
      write(*,*) 'Difference = ',diff

C     TEST FEASIBILITY ( b - Ax >= 0 )

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

      actconstr = 0
      bmin = 1.0d+99
      do i = 1,m
          bmin = min( bmin, b(i) )
          if ( b(i) .le. 1.0d-6 ) then
              actconstr = actconstr + 1
          end if
      end do

      write(*,*)
      write(*,*) 'If this is positive, the final point is feasible',bmin
      write(*,*) 'Number of active constraints = ',actconstr

      open(30, file='ivmline.txt')
      write(30,fmt=9000) pname,n,m,outit,fcnt,gcnt,hcnt,modh,innit,time,
     +                   f,dsupn,sqrt(deucn2),amax,flag,actconstr
      close(30)

      stop

C     NON EXECUTABLE STATEMENTS

 9000 format(A8,' &',I7,' &',I7,' &',I6,' &',I6,' &',I6,' &',I6,' &',
     +       I6,' &',I6,' &',F9.2,' &',1PD14.6,' &',1PD9.1,' &',1PD9.1,
     +       ' &',1PD9.1,' &',I1,' &',I7)
 9010 format(6(1X,1P,D12.4))

      end

C     ******************************************************
C     ******************************************************

      subroutine evalf(n,x,f,inform)

C     This subroutine computes the objective function.
C
C     On Entry:
C
C     n     integer,
C           number of variables,
C
C     x     double precision x(n),
C           point at which the function will be evaluated.
C
C     On Return
C
C     f     double precision,
C           function value at x,
C
C     inform integer,
C           termination parameter:
C           0= the function was successfully evaluated,
C           1= some error occurs in the function evaluation.

C     SCALAR ARGUMENTS
      double precision f
      integer n,inform

C     ARRAY ARGUMENTS
      double precision x(n)

C     LOCAL SCALARS
      double precision dist
      integer i,j,npol

C     INTRINSIC FUNCTIONS
      intrinsic sqrt

      inform = 0

      npol = n / 3 - 1

      f = 0.0d0

      do i = 1,npol

          dist = 0.0d0
          do j = 1,3
              dist = dist + (x(3 * npol + j) - x(3 * (i - 1) + j)) ** 2 
          end do
          dist = sqrt( dist )

          if ( dist .le. 1.0d-4 ) then
              write (*,fmt=*)
     +          'ERROR IN PROBLEM DEFINITION (DIST TOO SMALL)'
              inform = 1
          end if

          f = f + dist

      end do

      f = f / npol

      end


C     ******************************************************
C     ******************************************************

      subroutine evalg(n,x,g,inform)

C     This subroutine computes the gradient of the
C     objective function.
C
C     On Entry:
C
C     n     integer,
C           number of variables,
C
C     x     double precision x(n),
C           point at which the gradient will be evaluated.
C
C     On Return
C
C     g     double precision g(n),
C           gradient vector at x,
C
C     inform integer,
C           termination parameter:
C           0= the gradient was successfully evaluated,
C           1= some error occurs in the gradient evaluation.

C     SCALAR ARGUMENTS
      integer n,inform

C     ARRAY ARGUMENTS
      double precision g(n),x(n)

C     LOCAL SCALARS
      double precision dist
      integer i,j,npol

C     LOCAL ARRAYS
      double precision diff(3)

C     INTRINSIC FUNCTIONS
      intrinsic sqrt

      inform = 0

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

      npol = n / 3 - 1

      do i = 1,npol

          dist = 0.0d0
          do j = 1,3
              diff(j) = x(3 * npol + j) - x(3 * (i - 1) + j)
              dist = dist + diff(j) ** 2
          end do
          dist = sqrt( dist )

          do j = 1,3
              g(3 * (i - 1) + j) = - ( diff(j) / dist ) / npol
              g(3 * npol + j) = g(3 * npol + j) - g(3 * (i - 1) + j)
          end do

      end do

      end


C     ******************************************************
C     ******************************************************

      subroutine evalH(n,x,ldh,H,inform)

C     This subroutine computes the hessian of the 
C     objective function.
C
C     On Entry:
C
C     n     integer,
C           number of variables,
C
C     x     double precision x(n),
C           point at which the Hessian will be evaluated,
C
C     lda   integer,
C           leading dimension of matrix H.
C
C     On Return
C
C     H     double precision H(nmax,n),
C           Hessian matrix at x,
C
C     inform integer,
C           termination parameter:
C           0= the Hessian was successfully evaluated,
C           1= some error occurs in the Hessian evaluation.

C     SCALAR ARGUMENTS
      integer n,ldh,inform

C     ARRAY ARGUMENTS
      double precision H(ldh,n),x(n)

C     LOCAL SCALARS
      integer i,j,k,npol,ind1,ind2,ind3,ind4
      double precision dist,dist2

C     LOCAL ARRAYS
      double precision diff(3)

      inform = 0

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

      npol = n / 3 - 1

      do i = 1,npol

          dist2 = 0.0d0
          do j = 1,3
              diff(j) = x(3 * npol + j) - x(3 * (i - 1) + j)
              dist2 = dist2 + diff(j) ** 2
          end do
          dist = sqrt( dist2 ) * npol

          do j = 1,3
              ind1 = 3 * (i - 1) + j
              do k = 1,3
                  ind2 = 3 * (i - 1) + k
                  H(ind1,ind2) = - diff(j) * diff(k) / dist / dist2
                  if ( j .eq. k ) then
                      H(ind1,ind2) = H(ind1,ind2) + 1.0d0 / dist
                  end if
                  ind3 = 3 * npol + j
                  ind4 = 3 * npol + k
                  H(ind1,ind4) = - H(ind1,ind2)
                  H(ind3,ind2) = - H(ind1,ind2)
                  H(ind3,ind4) =   H(ind3,ind4) + H(ind1,ind2)
              end do
          end do

      end do

      end

C     ******************************************************
C     ******************************************************

      subroutine genpro(ngp,step,prob,rmin,rmax,ntpmin,ntpmax,m,n,
     +                  arow,acol,aval,anze,b,l,u,x,sol)

C     This subroutine generates a 3D-location problem. First, a regular 
C     grid with nx, ny and nz points at edges x, y and z, respectively,
C     in the positive orthant is considered. The points of the grid start 
C     at the origin with an horizontal distances of xstep, ystep and zstep,
C     respectively. This grid will be the space at which the cities 
C     (represented by polygons) will be distributed. Before building the 
C     cities, an area (rectangular polygon) of preservation where 
C     almost nothing can be done, is defined. This area of preservation 
C     will receive, after the construction of the cities, an hydraulic 
C     plant of energy generation (to supply the energy to the cities). Then, 
C     in the rest of the space, the cities are built. At each point of the 
C     grid (out of the central region) a city (represented by a polygon)
C     will be built with probability prob. The definition of the polygon
C     uses variables ntpmin, ntpmax, rmin and rmax in a way described in
C     the genpol (generate polygon) subroutine. To transmit the energy
C     from the plant to the cities, a tower inside each city and a tower
C     inside the central region must be built. The objective of the
C     problem is to determine the location of this towers in order 
C     to minimize the sum of the distances from each city tower to the
C     central one. 
C
C     On Entry:
C
C     ngp   integer ngp(3),
C           number grid points at each one of the 3 dimensions,
C
C     step  double precision step(3),
C           distance between points of the grid at each dimension,
C
C     prob  double precision, 
C           probability of defining a city at a point of the grid
C           (0 <= prob <= 1),
C
C     ntpmin integer,
C           minimum number of half spaces for each polygon,
C
C     ntpmax integer,
C           maximum number of half spaces for each polygon,
C
C     rmin  double precision,
C           minimum ratio of the sphere contained inside each polygon,
C
C     rmax  double precision,
C           maximum ratio of the sphere contained inside each polygon.
C
C     On output:
C
C     m     integer,
C           number of constraints,
C
C     n     integer,
C           number of variables,
C
C     arow  integer arow(mmax*nmax),
C     acol  integer acol(mmax*nmax),
C     aval  double precision aval(mmax*nmax),
C           coefficients matrix of the constraints in sparse format,
C
C     anze  integer,
C           number of nonzero elements of the coefficients matrix,
C
C     b     double precision b(mmax),
C           right hand side of the linear inequality constraints.
C
C     x     double precision x(nmax),
C           an interior feasible point.

C     PARAMETERS
      integer mmax
      parameter (mmax    = 3 650 000)
      integer nmax
      parameter (nmax    =     4 000)
      integer anzemax
      parameter (anzemax =  3 * mmax)

C     SCALAR ARGUMENTS
      double precision prob,rmax,rmin
      integer m,n,ntpmin,ntpmax,anze

C     ARRAY ARGUMENTS
      integer ngp(3),arow(anzemax),acol(anzemax)
      double precision step(3),aval(anzemax),b(mmax),l(nmax),u(nmax),
     +        x(nmax),sol(nmax)

C     LOCAL SCALARS
      double precision r,seed
      integer i,j,k,nac,q

C     LOCAL ARRAYS
      double precision c(3),cbox(3),cpl(3),cpu(3)

C     EXTERNAL FUNCTIONS
      double precision drand
      external drand

C     EXTERNAL SUBROUTINES
      external genpol

C     INITIALIZE NUMBER OF VARIABLES AND CONSTRAINTS

      m = 0
      n = 0

C     DEFINE CENTRAL POLYGON. THE SIX CONSTRAINTS WHICH DEFINE THE
C     CENTRAL POLYGON ARE, IN FACT, BOUND CONSTRAINTS.

      do i = 1,3

C         The lower bound for component i
          cpl(i) = 0.4d0 * ( ngp(i) - 1 ) * step(i)

C         The upper bound for component i
          cpu(i) = 0.6d0 * ( ngp(i) - 1 ) * step(i)

C         The center of the central polygon
          cbox(i)  = ( cpl(i) + cpu(i) ) / 2.0d0
      end do

C     SEED FOR THE RANDOM GENERATION

      seed = 760013

C     DEFINE CITY-POLYGONS CENTERED AT THE GRID POINTS AND 'SAFELY'
C     OUTSIDE THE CENTRAL REGION. SAFELY MEANS rmax FAR FROM THE 
C     CENTRAL REGION.

      anze = 0

      do i = 0,ngp(1) - 1
          c(1) = i * step(1)

          do j = 0,ngp(2) - 1
              c(2) = j * step(2)

              do k = 0,ngp(3) - 1
                  c(3) = k * step(3)

                  if ( c(1) .le. cpl(1) - 2.0d0 * rmax .or. 
     +                 c(1) .ge. cpu(1) + 2.0d0 * rmax .or. 
     +                 c(2) .le. cpl(2) - 2.0d0 * rmax .or. 
     +                 c(2) .ge. cpu(2) + 2.0d0 * rmax .or. 
     +                 c(3) .le. cpl(3) - 2.0d0 * rmax .or. 
     +                 c(3) .ge. cpu(3) + 2.0d0 * rmax ) then

C                      GENERATE A NEW POLYGON

                       if ( drand( seed ) .le. prob ) then

c                      if ( 
c    +
c    +     ( c(1).eq.2.0d0.and.c(2).eq.5.0d0.and.c(3).eq.5.0d0 ) .or.
c    +     ( c(1).eq.8.0d0.and.c(2).eq.5.0d0.and.c(3).eq.5.0d0 ) .or.
c    +     ( c(1).eq.5.0d0.and.c(2).eq.2.0d0.and.c(3).eq.5.0d0 ) .or.
c    +     ( c(1).eq.5.0d0.and.c(2).eq.8.0d0.and.c(3).eq.5.0d0 ) .or.
c    +     ( c(1).eq.5.0d0.and.c(2).eq.5.0d0.and.c(3).eq.2.0d0 ) .or.
c    +     ( c(1).eq.5.0d0.and.c(2).eq.5.0d0.and.c(3).eq.8.0d0 ) 
c    +
c    +                    ) then

C                         SET THE COMPONENTS OF THE FEASIBLE INTERIOR POINT

                          if ( n + 3 .gt. nmax ) then
                              write(*,*) 'Increase nmax and run again'
                              stop
                          end if

                          do q = 1,3
                              x(n+q) = c(q)
                          end do

C                         GENERATE THE RATIO OF THE SPHERE

                          r = rmin + ( rmax - rmin ) * drand( seed )

C                         GENERATE THE SIX CONSTRAINTS WHICH DETERMINE 
C                         THE SMALLEST BOX WHICH CONTAINS THE SPHERE OF 
C                         RATIO r CENTERED AT c. THIS SIX CONSTRAINTS 
C                         ARE, IN FACT, BOUND CONSTRAINTS

                          do q = 1,3

C                             The lower bound for component i
                              l(n+q) = c(q) - r

C                             The upper bound for component i
                              u(n+q) = c(q) + r

                          end do

C                         GENERATE THE OTHER CONSTRAINTS WHICH DEFINE
C                         THE POLYGON

                          call genpol(cbox,c,r,ntpmin,ntpmax,seed,m,n,
     +                                arow,acol,aval,anze,b,nac,sol)

C                         INCREASE THE NUMBER OF VARIABLES AND CONSTRAINTS

                          m = m + nac
                          n = n + 3

                      end if
                  end if
              end do
          end do
      end do

C     DEFINE CENTRAL POLYGON. THE SIX CONSTRAINTS WHICH DEFINE THE
C     CENTRAL POLYGON ARE, IN FACT, BOUND CONSTRAINTS.

      do i = 1,3

C         The lower bound for component i
          l(n+i) = cpl(i)

C         The upper bound for component i
          u(n+i) = cpu(i)

          x(n+i) = ( l(n+i) + u(n+i) ) / 2.0d0

          sol(n+i) = x(n+i)
      end do

C     INCREASE THE NUMBER OF VARIABLES

      n = n + 3

      end


C     ******************************************************
C     ******************************************************

      subroutine genpol(cbox,c,r,ntpmin,ntpmax,seed,m,n,arow,acol,aval,
     +                  anze,b,nac,sol)

C     This subroutine generates a polygon in R^3 as the intersection
C     of half spaces determined by planes tangent to the sphere with
C     center c and ratio r. The number of tangent planes is randomly 
C     generated satisfying ntpmin <= number of tangent planes <= ntpmax. 
C 
C     On Entry:
C
C     c     double precision c(3),
C           center of the sphere,
C
C     r     double precision,
C           ratio of the sphere,
C
C     ntpmin integer,
C           minimum number of tangent planes,
C
C     ntpmax integer,
C           maximum number of tangent planes,
C
C     seed  double precision,
C           seed for the random generation,
C
C     m     integer,
C           current number of constraints,
C
C     n     integer,
C           current number of variables.
C
C     arow  integer arow(mmax*nmax),
C     acol  integer acol(mmax*nmax),
C     aval  double precision aval(mmax*nmax),
C           coefficients matrix of the current constraints in 
C           sparse format,
C
C     anze  integer,
C           current number of non zero elements of matrix A,
C
C     b     double precision b(mmax),
C           right hand side of the current linear inequality 
C           constraints A x <= b.
C
C     On Return:
C
C     arow  integer arow(mmax*nmax),
C     acol  integer acol(mmax*nmax),
C     aval  double precision aval(mmax*nmax),
C           coefficients matrix of the current constraints in 
C           sparse format,
C
C     anze  integer,
C           current number of non zero elements of matrix A,
C
C     b     double precision b(mmax),
C           right hand side of the current linear inequality 
C           constraints A x <= b,
C
C     nac  integer,
C          number of added constraints.

C     PARAMETERS
      integer mmax
      parameter (mmax    = 3 650 000)
      integer nmax
      parameter (nmax    =     4 000)
      integer anzemax
      parameter (anzemax =  3 * mmax)

C     SCALAR ARGUMENTS
      integer ntpmin,ntpmax,m,n,anze,nac
      double precision r,seed

C     ARRAY ARGUMENTS
      integer arow(anzemax),acol(anzemax)
      double precision cbox(3),c(3),aval(anzemax),b(mmax),sol(nmax)

C     LOCAL SCALARS
      double precision penorm,rhs
      integer i,j

C     LOCAL ARRAYS
      double precision p(3)

C     EXTERNAL FUNCTIONS
      double precision drand
      external drand

C     INTRINSIC FUNCTIONS
      intrinsic int,sqrt

C     GENERATE THE NUMBER OF TANGENT PLANES 

      nac = ntpmin + int( ( ntpmax - ntpmin + 1 ) * drand( seed ) )

      if ( m + nac .gt. mmax ) then
          write(*,*) 'Increase mmax and run again'
          stop
      end if

C     GENERATE THE TANGENT PLANES

C     GENERATE A SPECIFIC PLANE

      if ( nac  .ge. 1 ) then

          penorm = 0.0d0
          do j = 1,3
              p(j) = cbox(j) - c(j)
              penorm = penorm + p(j) ** 2
          end do
          penorm = sqrt( penorm )

          rhs = r ** 2
          do j = 1,3
              p(j) = p(j) / penorm * r
              rhs = rhs + p(j) * c(j) 
          end do

          do j = 1,3
              sol(n+j) = c(j) + p(j)
          end do

          do j = 1,3
              anze = anze + 1

              if ( anze .gt. anzemax ) then
                  write(*,*) 'Increase anzemax and run again'
                  stop
              end if

              arow(anze) = m + 1
              acol(anze) = n + j
              aval(anze) = p(j)
          end do

      end if

C     SET THE RIGHT HAND SIDE VECTOR b

      b(m+1) = rhs

      do i = 2,nac

C         GENERATE A POINT p IN THE BORDER OF THE SPHERE OF RATIO r
C         CENTERED AT THE ORIGIN. THE PLANE TANGENT TO THAT POINT,
C         TRANSLATED TO THE SPHERE OF CENTER c IS DEFINED BY
C         p^t x = p^t ( p + c ) = r^2 + p^t c. THE HALF SPACE WHICH
C         CONTAINS c IS p^t x <= r^2 + p^t c.

          penorm = 0.0d0
          do j = 1,3
              p(j) = drand(seed)
              penorm = penorm + p(j) ** 2
          end do
          penorm = sqrt( penorm )

          rhs = r ** 2
          do j = 1,3
              p(j) = p(j) / penorm * r
              rhs = rhs + p(j) * c(j) 
          end do

C         ADD THE COEFFICIENT TO THE SPARSE MATRIX A

          do j = 1,3
              anze = anze + 1

              if ( anze .gt. anzemax ) then
                  write(*,*) 'Increase anzemax and run again'
                  stop
              end if

              arow(anze) = m + i
              acol(anze) = n + j
              aval(anze) = p(j)
          end do

C         SET THE RIGHT HAND SIDE VECTOR b

          b(m+i) = rhs

      end do

      end

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

      double precision function drand(ix)
      double precision ix

      double precision a,p,b15,b16,xhi,xalo,leftlo,fhi,k
      data a/16807.d0/,b15/32768.d0/,b16/65536.d0/,p/2147483647.d0/

      xhi= ix/b16
      xhi= xhi - dmod(xhi,1.d0)
      xalo= (ix-xhi*b16)*a
      leftlo= xalo/b16
      leftlo= leftlo - dmod(leftlo,1.d0)
      fhi= xhi*a + leftlo
      k= fhi/b15
      k= k - dmod(k,1.d0)
      ix= (((xalo-leftlo*b16)-p)+(fhi-k*b15)*b16)+k
      if (ix.lt.0) ix= ix + p
      drand= ix*4.656612875d-10

      return
      end


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

      subroutine iniivm(epsdsn,epsden,epsnlsn,maxoutit,maxinnit,maxfc,p,
     +                  extrap,hess,output,kprint,eta,theta,tau,tol,
     +                  lmin,lmax,gamma,sterel,steabs,epsmac,infty,tau1,
     +                  tau2)

C     This subroutine set the IVM parameters with its default values

C     SCALAR ARGUMENTS
      logical hess,extrap
      integer kprint,maxfc,maxinnit,maxoutit,output,p
      double precision epsdsn,epsden,epsnlsn,eta,theta,tau,tol,lmin,
     +        lmax,gamma,tau1,tau2,sterel,steabs,epsmac,infty

      epsdsn   =  1.0d-06
      epsden   =  1.0d-06
      epsnlsn  =  1.0d-08
      maxoutit =      100 
      maxinnit =   100000   
      maxfc    =     2000
      p        =       10
      extrap   =    .true.
      hess     =    .true.
      output   =        3
      kprint   =        5
c     eta      =   0.68d0
c     theta    =   0.85d0
      eta      =   0.90d0
      theta    =   0.95d0
      tau      =  1.0d-01
      tol      =   0.99d0
      lmin     =  1.0d-03
      lmax     =  1.0d+03
      gamma    =  1.0d-04
      sterel   =  1.0d-07
      steabs   =  1.0d-10
      infty    =  1.0d+20

      call mcheps(epsmac)

      tau1 = epsmac ** ( 1.0d0 / 3.0d0 )
      tau2 = epsmac ** ( 1.0d0 / 3.0d0 )

      end

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

      subroutine mcheps(epsmac)

C     This subroutine computes an approximation of machine epsilon.

C     SCALAR ARGUMENTS
      double precision epsmac

      epsmac = 1.0d0

 20   if ( ( 1.0d0 + epsmac ) .ne. 1.0d0 ) then
          epsmac = epsmac / 2.0d0
          go to 20
      end if

      epsmac = 2.0d0 * epsmac

      end
