C     =================================================================
C     File: sinterface.f
C     =================================================================

C     =================================================================
C     Module: Interface to convert to just-equality-constraints type
C     =================================================================

C     Last update of any of the component of this module: 
 
C     November 27, 2006.

C     Users are encouraged to download periodically updated versions of 
C     this code at the COLLECTION home page:
C
C     www.ime.usp.br/~egbirgin/collection/
C
C     and periodically updated versions of the TANGO Project solvers at
C     the TANGO home page:
C 
C     www.ime.usp.br/~egbirgin/tango/ 

C     =================================================================

C     This interface can be used to solve problems from LA CUMPARSITA
C     collection with solvers that deals just with equality constrained
C     problems. For that purpouse, this interface adds slack variables
C     to the problem definition and makes all the necessary 
C     transformations to the output of subroutines inip, evalf, evalg, 
C     evalh, evalc, evaljac, evalhc and endp. 

C     Slack variables are initialized evaluating the constraints of the 
C     original problem and giving to the slack variables the value 
C     necessary to satisfy the reformulated inequality constraints.

C     In addition, this interface also provides a subroutine to compute
C     the sparse Hessian of the Lagrangian, using the subroutines
C     provided by the user to compute the Hessian of the objective
C     function and the individual Hessians of the constraints. The
C     Hessian of the Lagrangian, which consists of the weighted sum of
C     Hessian of the objective function and the Hessians of the 
C     constraints, is computed with linear complexity in time (with
C     respect to its number of non null elements) but uses a dense
C     n x n integer matrix as auxiliary space.

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

      subroutine inips(n,x,l,u,m,lambda,rho,equatn,linear)

      implicit none

C     SCALAR ARGUMENTS
      integer m,n

C     ARRAY ARGUMENTS
      logical equatn(*),linear(*)
      double precision l(*),lambda(*),rho(*),u(*),x(*)

C     PARAMETERS
      integer mmax,nsmax
      parameter ( mmax      =  500000 )
      parameter ( nsmax     =   10000 )

C     COMMON SCALARS
      logical comphl
      integer no

C     COMMON ARRAYS
      integer sindc(mmax),pos(nsmax,nsmax)

C     LOCAL SCALARS
      integer flag,i,j
      double precision c

C     COMMON BLOCKS
      common /sladat/ sindc,no
      save   /sladat/

      common /hldata/ pos,comphl
      save   /hldata/

      call inip(no,x,l,u,m,lambda,rho,equatn,linear)

C     Add slacks for inequality constraints. Constraints of the form
C     c(x) <= 0 become c(x) + s = 0, s >= 0.

      n = no

      do i = 1,m
          if ( .not. equatn(i) ) then

C             Convert into an equality constraint
              equatn(i) = .true.

C             Add a slack variable
              n = n + 1
              sindc(i) = n

C             Add its bounds
              l(n) = 0.0d0
              u(n) = 1.0d+20

C             Initialize it
              call evalc(no,x,i,c,flag)
              if ( flag .ne. 0 ) then
                  return
              end if

              x(n) = max( l(n), min( - c, u(n) ) )

          else
              sindc(i) = 0
          end if
      end do

C     Initialize structure to compute the Hessian of the Lagrangian.

      if ( no .gt. nsmax ) then
          comphl = .false.
          write(*,100)

      else
          comphl = .true.
          do j = 1,no
              do i = 1,no
                  pos(i,j) = 0
              end do
          end do
      end if

C     Non-executable statements

 100  format(/,1X,'WARNING: There will be not enough space to compute ',
     +         'the sparse Hessian of the',/,1X,'Lagrangian. If your ',
     +         'solver needs it (ALGENCAN does not) then increase the ',
     +         'value of',/,1X,'nsmax in the file sinterface.f, ',
     +         'compile and run again.')

      end

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

      subroutine evalfs(n,x,f,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,n
      double precision f

C     ARRAY ARGUMENTS
      double precision x(n)

C     PARAMETERS
      integer mmax
      parameter ( mmax      =  500000 )

C     COMMON SCALARS
      integer no

C     COMMON ARRAYS
      integer sindc(mmax)

C     COMMON BLOCKS
      common /sladat/ sindc,no
      save   /sladat/

      call evalf(no,x,f,flag)

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalgs(n,x,g,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,n

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

C     PARAMETERS
      integer mmax
      parameter ( mmax      =  500000 )

C     COMMON SCALARS
      integer no

C     COMMON ARRAYS
      integer sindc(mmax)

C     LOCAL SCALARS
      integer i

C     COMMON BLOCKS
      common /sladat/ sindc,no
      save   /sladat/

      call evalg(no,x,g,flag)

      do i = no + 1,n
          g(i) = 0.0d0
      end do

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalhs(n,x,hlin,hcol,hval,hnnz,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,n,hnnz

C     ARRAY ARGUMENTS
      integer hcol(*),hlin(*)
      double precision hval(*),x(n)

C     PARAMETERS
      integer mmax
      parameter ( mmax      =  500000 )

C     COMMON SCALARS
      integer no

C     COMMON ARRAYS
      integer sindc(mmax)

C     COMMON BLOCKS
      common /sladat/ sindc,no
      save   /sladat/

      call evalh(no,x,hlin,hcol,hval,hnnz,flag)

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalcs(n,x,ind,c,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer ind,flag,n
      double precision c

C     ARRAY ARGUMENTS
      double precision x(n)

C     PARAMETERS
      integer mmax
      parameter ( mmax      =  500000 )

C     COMMON SCALARS
      integer no

C     COMMON ARRAYS
      integer sindc(mmax)

C     COMMON BLOCKS
      common /sladat/ sindc,no
      save   /sladat/

      call evalc(no,x,ind,c,flag)

      if ( sindc(ind) .ne. 0 ) then
          c = c + x(sindc(ind))
      end if

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evaljacs(n,x,ind,jcvar,jcval,jcnnz,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,ind,n,jcnnz

C     ARRAY ARGUMENTS
      integer jcvar(n)
      double precision x(n),jcval(n)

C     PARAMETERS
      integer mmax
      parameter ( mmax      =  500000 )

C     COMMON SCALARS
      integer no

C     COMMON ARRAYS
      integer sindc(mmax)

C     COMMON BLOCKS
      common /sladat/ sindc,no
      save   /sladat/

      call evaljac(no,x,ind,jcvar,jcval,jcnnz,flag)

      if ( sindc(ind) .ne. 0 ) then
          jcnnz = jcnnz + 1
          jcvar(jcnnz) = sindc(ind)
          jcval(jcnnz) = 1.0d0
      end if

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalhcs(n,x,ind,hclin,hccol,hcval,hcnnz,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,ind,n,hcnnz

C     ARRAY ARGUMENTS
      integer hccol(*),hclin(*)
      double precision hcval(*),x(n)

C     PARAMETERS
      integer mmax
      parameter ( mmax      =  500000 )

C     COMMON SCALARS
      integer no

C     COMMON ARRAYS
      integer sindc(mmax)

C     COMMON BLOCKS
      common /sladat/ sindc,no
      save   /sladat/

      call evalhc(no,x,ind,hclin,hccol,hcval,hcnnz,flag)

      end

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

      subroutine evalhls(n,x,m,lambda,hllin,hlcol,hlval,hlnnz,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,hlnnz,m,n

C     ARRAY ARGUMENTS
      integer hllin(*),hlcol(*)
      double precision hlval(*),lambda(m),x(n)

C     This subroutine computes the Hessian of the lagrangian. Its 
C     memory complexity is n^2 while its time complexity is linear
C     in the number of non null elements of the Hessian being 
C     computed. It is assumed that pos is an n x n integer matrix 
C     with all its components initialized with zero.
 
C     PARAMETERS
      integer hnnzmax,nsmax
      parameter ( hnnzmax   = 5000000 )
      parameter ( nsmax     =   10000 )

C     COMMON SCALARS
      logical comphl

C     COMMON ARRAYS
      integer pos(nsmax,nsmax)

C     LOCAL SCALARS
      integer i,j,k,hnnz
      double precision mult

C     LOCAL ARRAYS
      integer hlin(hnnzmax),hcol(hnnzmax)
      double precision hval(hnnzmax)

C     COMMON BLOCKS
      common /hldata/ pos,comphl
      save   /hldata/

C     Verify whether there is enough space or not.

      if ( .not. comphl ) then
          write(*,100)
          flag = - 1
          return
      end if

C     Compute Hessian of the Lagrangian.

      hlnnz = 0

      do j = 0,m
          if ( j .eq. 0 ) then
              mult = 1.0d0
              call evalhs(n,x,hlin,hcol,hval,hnnz,flag)
              if ( flag .ne. 0 ) then
                  return
              end if
          else
              mult = lambda(j)
              call evalhcs(n,x,j,hlin,hcol,hval,hnnz,flag)
              if ( flag .ne. 0 ) then
                  return
              end if
          end if

          do i = 1,hnnz
              k = pos(hlin(i),hcol(i))
              if ( k .ne. 0 ) then
                  hlval(k) = hlval(k) + hval(i) * mult
              else
                  hlnnz = hlnnz + 1
                  pos(hlin(i),hcol(i)) = hlnnz

                  k = pos(hlin(i),hcol(i))
                  hllin(k) = hlin(i)
                  hlcol(k) = hcol(i)
                  hlval(k) = hval(i) * mult
              end if
          end do
      end do

C     Clean matrix pos.

      do i = 1,hlnnz
          pos(hllin(i),hlcol(i)) = 0
      end do

C     Non-executable statements

 100  format(/,1X,'ERROR: Your solver called subroutine evalhls to ',
     +         'compute the sparse Hessian of the',/,1X,'Lagrangian. ',
     +         'This subroutine has a time complexity that is linear ',
     +         'in the number of non',/,1X,'null elements of the ',
     +         'Hessian being computed but uses n^2 space. Increase ',
     +         'the value of',/,1X,'nsmax within file sinterface.f, ',
     +         'compile and run again.')

      end

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

      subroutine evalhlps(n,x,m,lambda,p,hp,goth,flag)

      implicit none

C     SCALAR ARGUMENTS
      logical goth
      integer flag,m,n

C     ARRAY ARGUMENTS
      double precision hp(n),lambda(m),p(n),x(n)

C     PARAMETERS
      integer mmax
      parameter ( mmax      =  500000 )
C     integer hnnzmax
C     parameter ( hnnzmax   = 5000000 )

C     COMMON SCALARS
      integer no

C     COMMON ARRAYS
      integer sindc(mmax)

C     LOCAL SCALARS
C     integer i,hlnnz,lin,col
C     double precision val

C     LOCAL ARRAYS
C     integer hllin(hnnzmax),hlcol(hnnzmax)
C     double precision hlval(hnnzmax)

C     COMMON BLOCKS
      common /sladat/ sindc,no
      save   /sladat/

      call evalhlp(no,x,m,lambda,p,hp,goth,flag)

C     The code below can be used to test the correctness of 
C     subroutine evalhls.

C     if ( .not. goth ) then
C         goth = .true.

C         call evalhls(n,x,m,lambda,hllin,hlcol,hlval,hlnnz,flag)
C         if ( flag .ne. 0 ) then
C             return
C         end if
C     end if

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

C     do i = 1,hlnnz
C         lin = hllin(i)
C         col = hlcol(i)
C         val = hlval(i)

C         hp(lin) = hp(lin) + val * p(col)
C         if ( lin .ne.col ) then
C             hp(col) = hp(col) + val * p(lin)
C         end if
C     end do

      end

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

      subroutine endps(n,x,l,u,m,lambda,rho,equatn,linear)

      implicit none

C     SCALAR ARGUMENTS
      integer m,n

C     ARRAY ARGUMENTS
      logical equatn(m),linear(m)
      double precision l(n),lambda(m),rho(m),u(n),x(n)

C     PARAMETERS
      integer mmax
      parameter ( mmax      =  500000 )

C     COMMON SCALARS
      integer no

C     COMMON ARRAYS
      integer sindc(mmax)

C     COMMON BLOCKS
      common /sladat/ sindc,no
      save   /sladat/

      call endp(no,x,l,u,m,lambda,rho,equatn,linear)

      end

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

      subroutine endpr(n,x,l,u,m,lambda,rho,equatn,linear)

      implicit none

C     SCALAR ARGUMENTS
      integer m,n

C     ARRAY ARGUMENTS
      logical equatn(m),linear(m)
      double precision l(n),lambda(m),rho(m),u(n),x(n)

C     This subroutine should be used by a solver that deals with
C     equality and inequality constraints (La Cumparsita original 
C     problem formulation) but needs to save the solution as it had 
C     been solved the problem using the interface subroutines of the
C     the present file, i.e., adding slacks variables. Moreover, this 
C     subroutine also computes and save the Lagrange multipliers for 
C     the bound constraints (including the bounds for the slacks).

C     PARAMETERS
      integer mmax,nmax
      parameter ( mmax      =  500000 )
      parameter ( nmax      =  500000 )

C     COMMON SCALARS
      integer no

C     COMMON ARRAYS
      integer sindc(mmax)

C     LOCAL SCALARS
      integer flag,i,j,jcnnz
      double precision c

C     LOCAL ARRAYS
      integer jcvar(nmax)
      double precision jcval(nmax),nl(nmax),llambda(nmax),ulambda(nmax)

C     COMMON BLOCKS
      common /sladat/ sindc,no
      save   /sladat/

C     Add slack variables

      no = n

      do i = 1,m
          if ( .not. equatn(i) ) then

C             Add a slack variable
              n = n + 1
              sindc(i) = n

C             Add its bounds
              l(n) = 0.0d0
              u(n) = 1.0d+20

C             Initialize it
              call evalc(no,x,i,c,flag)
              if ( flag .ne. 0 ) then
                  return
              end if

              x(n) = max( l(n), min( - c, u(n) ) )
          else
              sindc(i) = 0
          end if
      end do

C     Compute the gradient of the Lagrangian

      call evalgs(n,x,nl,flag)
      if ( flag .ne. 0 ) then
          return
      end if

      do j = 1,m
          call evaljacs(n,x,j,jcvar,jcval,jcnnz,flag)
          if ( flag .ne. 0 ) then
              return
          end if
      
          do i = 1,jcnnz
              nl(jcvar(i)) = nl(jcvar(i)) + jcval(i) * lambda(j)
          end do
      end do

C     Compute Lagrange multipliers of the bound constraints

      do i = 1,n
          llambda(i) = 0.0d0
          ulambda(i) = 0.0d0

          if ( x(i) .le. l(i) + 1.0d-04 ) then
              llambda(i) = nl(i)

          else if ( x(i) .ge. u(i) - 1.0d-04 ) then
              ulambda(i) = - nl(i)
          end if
      end do

C     Save solution

      open(10,file='inipoint.txt')

      write(10,*) n,m

      do i = 1,n
          write(10,100) i,x(i)
      end do

      do i = 1,m
          write(10,100) i,lambda(i)
      end do

      do i = 1,n
          if ( l(i) .gt. - 1.0d+20 ) then
              write(10,100) i,llambda(i)
          end if
      end do

      do i = 1,n
          if ( u(i) .lt. 1.0d+20 ) then
              write(10,100) i,ulambda(i)
          end if
      end do

      close(10)

C     Restore original number of variables

      n = no

C     Non-executable statements

 100  format(I7,1X,1P,D24.16)

      end
