C     =================================================================
C     File: ellipsvm.f
C     =================================================================

C     =================================================================
C     Module: Spectral Projected Gradient Method. Problem definition.
C     =================================================================

C     Last update of any of the component of this module:

C     October 5, 2012.

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

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

C
C     Compile with: 
C
C     gfortran -O4 -xf77-cpp-input spgma.f spg.f ellipsvm.f -llapack \
C              -o spgma-ellipsvm
C 

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

      subroutine inip(n,x)

      implicit none

C     SCALAR ARGUMENTS
      integer n

C     ARRAY ARGUMENTS
      double precision x(n)

#include "pdatas.par"
#include "pdatas.com"

C     LOCAL SCALARS
      integer i,j
      double precision pnorm,pseed,seed

C     EXTERNAL FUNCTIONS
      double precision drand
      external drand

C     Problem data

      write(*,*) 'Enter dimension of the space (nmat <= nmatmax): '
      write(*,*) 'nmatmax = ',nmatmax
      read(*,*) nmat

      write(*,*) 'Enter the number of point (npmax = ',npmax,'): '
      read(*,*) np

      write(*,*) 'Enter seed for the points generation: '
      read(*,*) pseed

      write(*,*) 'Enter region [c]ircle, [s]quare, [r]ectangle, ',
     +           '[t]riangle (r and t only in 2D): '
      read(*,*) regtype

      if ( regtype .eq. 'c' ) then
          write(*,*) 'Enter circle radius: '
          read(*,*) rlim
      else if ( regtype .eq. 's' ) then
          write(*,*) 'Enter square half-size: '
          read(*,*) rlim
      else if ( regtype .eq. 'r' ) then
          write(*,*) 'Enter rectangle half-width: '
          read(*,*) wlim
          write(*,*) 'Enter rectangle half-height: '
          read(*,*) hlim
      else if ( regtype .eq. 't' ) then
          write(*,*) 'Enter triangle parameter r such that the ',
     +               'triangle corners are (-r,0), (0,-r), and (r,r): ' 
          read(*,*) rlim
      end if

      seed = pseed

      do j = 1,np
          do i = 1,nmat
              p(i,j) = -100.0d0 + 200.0d0 * drand(seed)
          end do

          if ( regtype .eq. 'c' ) then
              pnorm = 0.0d0
              do i = 1,nmat
                  pnorm = pnorm + p(i,j) ** 2
              end do
              if ( pnorm .le. rlim ** 2 ) then
                  t(j) = 'i'
              else
                  t(j) = 'o'
              end if

          else if ( regtype .eq. 's' ) then
              t(j) = 'i'
              do i = 1,nmat
                  if ( p(i,j) .lt. - rlim .or. p(i,j) .gt. rlim ) then
                      t(j) = 'o'
                  end if
              end do

          else if ( regtype .eq. 'r' ) then
              if ( - wlim .le. p(1,j) .and. p(1,j) .le. wlim .and.
     +             - hlim .le. p(2,j) .and. p(2,j) .le. hlim ) then
                  t(j) = 'i'
              else
                  t(j) = 'o'
              end if

          else if ( regtype .eq. 't' ) then
              if ( 2.0d0 * p(1,j) -         p(2,j) .le. rlim .and.
     +                   - p(1,j) + 2.0d0 * p(2,j) .le. rlim .and.
     +                   - p(1,j) -         p(2,j) .le. rlim ) then
                  t(j) = 'i'
              else
                  t(j) = 'o'
              end if
          end if
      end do

C     Number of variables
      n = nmat ** 2 + nmat

C     Initial point
      seed = 123456.0d0

      do i = 1,n
          x(i) = -1.0d0 + 2.0d0 * drand(seed)
      end do

      end

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

      subroutine evalf(n,x,f,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer n,flag
      double precision f

C     ARRAY ARGUMENTS
      double precision x(n)

#include "pdatas.par"
#include "pdatas.com"

C     COMMON ARRAYS
      double precision savedg(nmax)

C     COMMON BLOCKS
      common /savedgdat/ savedg
      save   /savedgdat/

      call evalfg(n,x,f,savedg,flag)

      end

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

      subroutine evalg(n,x,g,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer n,flag

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

#include "pdatas.par"
#include "pdatas.com"

C     COMMON ARRAYS
      double precision savedg(nmax)

C     LOCAL SCALARS
      integer i

C     COMMON BLOCKS
      common /savedgdat/ savedg
      save   /savedgdat/

      flag = 0

      do i = 1,n
          g(i) = savedg(i)
      end do

      end

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

      subroutine evalfg(n,x,f,g,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer n,flag
      double precision f

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

#include "pdatas.par"
#include "pdatas.com"

C     LOCAL SCALARS
      integer i,ind,j,k
      double precision fparc

C     LOCAL ARRAYS
      double precision tmp(nmatmax)

C     EXTERNAL FUNCTIONS
      integer Aind,bind
      external Aind,bind

      flag = 0

      f = 0.0d0

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

      do k = 1,np
          do i = 1,nmat
              tmp(i) = 0.0d0
          end do

          do j = 1,nmat
              do i = 1,nmat
                  tmp(i) = tmp(i) + x(Aind(i,j,nmat)) * p(j,k)
              end do
          end do

          fparc = - 1.0d0

          do i = 1,nmat
              fparc = fparc + tmp(i) * p(i,k)
          end do

          do i = 1,nmat
              fparc = fparc + x(bind(i,nmat)) * p(i,k)
          end do

          if ( ( t(k) .eq. 'i' .and. fparc .gt. 0.0d0 ) .or.
     +         ( t(k) .eq. 'o' .and. fparc .lt. 0.0d0 ) ) then

              f = f + fparc ** 2

              do i = 1,nmat
                  ind = bind(i,nmat)
                  g(ind) = g(ind) + 2.0d0 * fparc * p(i,k)
              end do

              do j = 1,nmat
                  do i = 1,nmat
                      ind = Aind(i,j,nmat)
                      g(ind) = g(ind) + 2.0d0 * fparc * p(i,k) * p(j,k)
                  end do
              end do

          end if
      end do

      f = f / np

      do i = 1,n
          g(i) = g(i) / np
      end do

      end

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

      integer function Aind(i,j,nmat)

      implicit none

C     SCALAR ARGUMENTS
      integer i,j,nmat

      Aind = ( j - 1 ) * nmat + i

      return

      end

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

      integer function bind(i,nmat)

      implicit none

C     SCALAR ARGUMENTS
      integer i,nmat

      bind = nmat ** 2 + i

      return

      end

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

      subroutine proj(n,x,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer n,flag

C     ARRAY ARGUMENTS
      double precision x(n)

#include "pdatas.par"
#include "pdatas.com"

C     LOCAL SCALARS
      character uplo,jobz
      integer i,ind,info,j,k,lwork

C     LOCAL ARRAYS
      double precision A(nmatmax,nmatmax),autval(nmatmax),
     +        work(10*nmatmax)

C     EXTERNAL SUBROUTINES
      external dsyev

C     EXTERNAL FUNCTIONS
      integer Aind,bind
      external Aind,bind

      flag = 0

      do j = 1,nmat
          do i = 1,nmat
              A(i,j) = 0.5d0 * ( x(Aind(i,j,nmat)) + x(Aind(j,i,nmat)) )
          end do
      end do

      uplo  = 'U'
      jobz  = 'V'
      lwork = 10 * nmat

      call dsyev(jobz,uplo,nmat,A,nmatmax,autval,work,lwork,info) 

      do j = 1,nmat
          do i = 1,nmat
              x(Aind(i,j,nmat)) = 0.0d0
          end do
      end do

      do k = 1,nmat
          autval(k) = max( 1.0d-04, min( autval(k), 1.0d+04 ) )
          do j = 1,nmat
              do i = 1,nmat
                  ind = Aind(i,j,nmat)
                  x(ind) = x(ind) + autval(k) * A(i,k) * A(j,k)
              end do
          end do
      end do

      end

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

      double precision function drand(ix)

      implicit none

C     This is the random number generator of Schrage:
C
C     L. Schrage, A more portable Fortran random number generator, ACM
C     Transactions on Mathematical Software 5 (1979), 132-138.

      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 drawsol(n,x)

      implicit none

C     SCALAR ARGUMENTS
      integer n

C     ARRAY ARGUMENTS
      double precision x(n)

C     This subroutine generate a metapost file with the graphical 
C     representation of the problem (just for nmat=2).

#include "pdatas.par"
#include "pdatas.com"

C     PARAMETERS
      double precision pi
      parameter ( pi = 3.14159265359d0 )

C     LOCAL SCALARS
      integer k
      double precision tmp,xmin,xmax,ymin,ymax,scale,aa,bb,cc,dd,alpha,
     +        seed,r1,r2,s1,s2

C     LOCAL ARRAYS
      double precision A(2,2),b(2),c(2),lambda(2),v(2),w(2)

C     EXTERNAL FUNCTIONS
      integer Aind,bind
      external Aind,bind

C     JUST FOR NDIM=2!

      if ( nmat .ne. 2 ) then
          return
      end if

C     SCALING

      xmin =   1.0d+99
      xmax = - 1.0d+99
      ymin =   1.0d+99
      ymax = - 1.0d+99

      do k = 1,np
          xmin = min( xmin, p(1,k) )
          xmax = max( xmax, p(1,k) )
          ymin = min( ymin, p(2,k) )
          ymax = max( ymax, p(2,k) )
      end do

      scale = min( 10.0d0 / ( xmax - xmin ), 10.0d0 / ( ymax - ymin ) )

C     DRAW THE PROBLEM

      open(unit=10,file='sol.mp')

      write(10,10)

C     ELLIPSE

C     SET A AND b

      A(1,1) = x(Aind(1,1,nmat))
      A(1,2) = x(Aind(1,2,nmat))
      A(2,1) = x(Aind(2,1,nmat))
      A(2,2) = x(Aind(2,2,nmat))

      b(1) = x(bind(1,nmat))
      b(2) = x(bind(2,nmat))

C     COMPUTE ELLIPSE'S CENTER (SOLUTION c OF 2 A c + b = 0)

      tmp = A(1,1) * A(2,2) - A(2,1) * A(1,2)

      if ( abs(A(1,1)) .gt. abs(A(2,2)) ) then
          c(2) = 0.5d0 * ( A(2,1) * b(1) - A(1,1) * b(2) ) / tmp
          c(1) = - ( A(1,2) * c(2) + 0.5d0 * b(1) ) / A(1,1)
      else
          c(1) = 0.5d0 * ( A(1,2) * b(2) - A(2,2) * b(1) ) / tmp 
          c(2) = - ( 0.5d0 * b(2) + A(2,1) * c(1) ) / A(2,2)
      end if

C     COMPUTE EIGENVALUES OF A

      aa = 1.0d0
      bb = - ( A(1,1) + A(2,2) )
      cc = tmp

      dd = sqrt( bb ** 2 - 4.0d0 * aa * cc )

      lambda(1) = ( - bb + dd ) / ( 2.0d0 * aa )
      lambda(2) = ( - bb - dd ) / ( 2.0d0 * aa )

C     COMPUTE EIGENVECTORS OF A

      if ( abs( A(2,1) ) .le. 1.0d-16 ) then
          if ( abs( lambda(1) - A(1,1) ) .le. 
     +         abs( lambda(1) - A(2,2) ) ) then
              v(1) = 1.0d0
              v(2) = 0.0d0
          else
              v(1) = 0.0d0
              v(2) = 1.0d0
          end if
      else
          if ( abs( A(2,1) ) .gt. abs( A(1,1) - lambda(1) ) ) then
              v(1) = 1.0d0
              v(2) = - ( A(1,1) - lambda(1) ) / A(2,1)
          else
              v(2) = 1.0d0
              v(1) = - A(2,1) / ( A(1,1) - lambda(1) )
          end if
      end if

      if ( abs( v(1) ) .gt. abs( v(2) ) ) then
          w(2) = 1.0d0
          w(1) = - v(2) * w(2) / v(1)
      else
          w(1) = 1.0d0
          w(2) = - v(1) * w(1) / v(2)
      end if

      tmp = sqrt( v(1) ** 2 + v(2) ** 2 )
      v(1) = v(1) / tmp
      v(2) = v(2) / tmp

      tmp = sqrt( w(1) ** 2 + w(2) ** 2 )
      w(1) = w(1) / tmp
      w(2) = w(2) / tmp

c     COMPUTE SEMIAXES

      aa = ( v(1)**2 + v(2)**2 ) * lambda(1)
      bb = v(1) * ( 2.0d0 * lambda(1) * c(1) + b(1) ) + 
     +     v(2) * ( 2.0d0 * lambda(1) * c(2) + b(2) )
      cc = c(1) * ( A(1,1) * c(1) + A(1,2) * c(2) + b(1) ) +
     +     c(2) * ( A(2,1) * c(1) + A(2,2) * c(2) + b(2) ) - 1.0d0

      dd = sqrt( bb ** 2 - 4.0d0 * aa * cc )

      r1 = ( - bb + dd ) / ( 2.0d0 * aa )
      r2 = ( - bb - dd ) / ( 2.0d0 * aa )

      s1 = abs(r1) * sqrt( v(1) ** 2 + v(2) ** 2 ) 

      aa = ( w(1)**2 + w(2)**2 ) * lambda(2)
      bb = w(1) * ( 2.0d0 * lambda(2) * c(1) + b(1) ) + 
     +     w(2) * ( 2.0d0 * lambda(2) * c(2) + b(2) )
      cc = c(1) * ( A(1,1) * c(1) + A(1,2) * c(2) + b(1) ) +
     +     c(2) * ( A(2,1) * c(1) + A(2,2) * c(2) + b(2) ) - 1.0d0

      dd = sqrt( bb ** 2 - 4.0d0 * aa * cc )

      r1 = ( - bb + dd ) / ( 2.0d0 * aa )
      r2 = ( - bb - dd ) / ( 2.0d0 * aa )

      s2 = abs(r1) * sqrt( w(1) ** 2 + w(2) ** 2 ) 

      alpha = acos(v(1)) / pi * 180.0d0
      alpha = sign(alpha,v(2))

C     DRAW ELLIPSE

      write(10,30) 2
      write(10,20) 2.0d0*s1*scale,2.0d0*s2*scale,alpha,
     +             c(1)*scale,c(2)*scale

C     REGION USED TO LABEL THE DOTS

      if ( regtype .eq. 'c' ) then
          write(10,30) 2
          write(10,20) 2.0d0*rlim*scale,
     +                 2.0d0*rlim*scale,0.0d0,0.0d0,0.0d0
      else if ( regtype .eq. 's' ) then
          write(10,30) 2
          write(10,50) -rlim*scale,-rlim*scale, rlim*scale,-rlim*scale,
     +                  rlim*scale, rlim*scale,-rlim*scale, rlim*scale 
      else if ( regtype .eq. 'r' ) then
          write(10,30) 2
          write(10,50) -wlim*scale,-hlim*scale, wlim*scale,-hlim*scale,
     +                  wlim*scale, hlim*scale,-wlim*scale, hlim*scale 
      else if ( regtype .eq. 't' ) then
          write(10,30) 2
          write(10,60) -rlim*scale,0.0d0,0.0d0,-rlim*scale,
     +                   rlim*scale,rlim*scale
      end if

C     DOTS WITHOUT LABELS

      write(10,30) 2

      do k = 1,min(np,10000)
          if ( t(k) .eq. 'i' ) then
              write(10,40) scale*p(1,k),scale*p(2,k),'blue '
           else
              write(10,40) scale*p(1,k),scale*p(2,k),'red  '
          end if
      end do

C     END

      write(10,100)

      close(10)

C     NON-EXECUTABLE STATEMENTS

  10  format('beginfig(1);'/,'u = 1.0 cm;') 
  20  format('draw fullcircle',
     +       /,5X,'xscaled  ',f20.10,'u',
     +       /,5X,'yscaled  ',f20.10,'u',
     +       /,5X,'rotated  ',f20.10,
     +       /,5X,'shifted (',f20.10,'u,',f20.10,'u);')
  30  format('pickup pencircle scaled ',I1,'pt;')
  40  format('draw (',f20.10,'u,',f20.10,'u) withcolor ',A5,';')
  50  format('draw (',f20.10,'u,',f20.10,'u)--',
     +            '(',f20.10,'u,',f20.10,'u)--',
     +            '(',f20.10,'u,',f20.10,'u)--',
     +            '(',f20.10,'u,',f20.10,'u)--cycle;')
  60  format('draw (',f20.10,'u,',f20.10,'u)--',
     +            '(',f20.10,'u,',f20.10,'u)--',
     +            '(',f20.10,'u,',f20.10,'u)--cycle;')
 100  format('endfig;'/,'end;') 

      end

