leastsq9.f90

Go to the documentation of this file.

!> @file leastsq9.f90
!! @ingroup FORT1THREAD_EXAMPLES
!!
!! 
!! The model is similar to leastsq7, but we test an alternative optimal
!! solution that was found at some stage with bounded residuals.
!! The new solution is then used as a starting point for a run without
!! bounded residuals and the solution should still be optimal.
!! 
!! For efficiency, the model is smaller.
!! 
!! We solve the following nonlinear least squares model:
!! 
!! \f[
!! \min \sum_i res_{i}^2 \\
!!
!! \sum_j ( a_{ij}x_j + b_{ij}x_j^2 ) + res_i = obs_i
!! \f]
!! 
!! where \f$a\f$, \f$b\f$, and \f$obs\f$ are known data, and \f$res\f$ and \f$x\f$ are the
!! variables of the model.
!! 
!!
!! For more information about the individual callbacks, please have a look at the source code.
 
#if defined(_WIN32) && !defined(_WIN64)
#define dec_directives_win32
#endif
 
module lsq9
   Integer, parameter :: Nobs = 14
   Integer, parameter :: DimX = 10
   real*8, dimension(Nobs,DimX) :: a, b
   real*8, dimension(Nobs)      :: obs
   Integer :: run ! Run = 1 standard,
                  ! Run = 2 with wide bounds on Res
                  ! Run = 3 standard model with obs(1) changed by +1.e-3
end module lsq9
 
 
REAL FUNCTION rndx( )
!
!  Defines a pseudo random number between 0 and 1
!
      IMPLICIT NONE
 
   Integer, save :: seed = 12359
 
   seed = mod(seed*1027+25,1048576)
   rndx = float(seed)/float(1048576)
 
END FUNCTION rndx
 
subroutine defdata
!
!  Define values for A, B, and Obs
!
   Use lsq9
      IMPLICIT NONE
 
   Integer :: i, j
   real*8, Parameter ::  xtarg = -1.0
   real*8, Parameter ::  noise = 1.0
   Real, External :: Rndx
   real*8 :: o
 
   do i = 1, nobs
      o = 0.d0
      do j = 1, dimx
         a(i,j) = rndx()
         b(i,j) = rndx()
         o = o + a(i,j) * xtarg + b(i,j) * xtarg**2
      enddo
      obs(i) = o + noise * rndx()
   enddo
 
end subroutine defdata
!
!  Main program.
!
!> Main program. A simple setup and call of CONOPT
!!
Program leastsquare
 
   Use proginfo
   Use conopt
   Use lsq9
      IMPLICIT NONE
!
!  Declare the user callback routines as Integer, External:
!
   Integer, External :: lsq_readmatrix ! Mandatory Matrix definition routine defined below
   Integer, External :: lsq_fdeval     ! Function and Derivative evaluation routine
                                       ! needed a nonlinear model.
   Integer, External :: std_status     ! Standard callback for displaying solution status
   Integer, External :: std_solution   ! Standard callback for displaying solution values
   Integer, External :: std_message    ! Standard callback for managing messages
   Integer, External :: std_errmsg     ! Standard callback for managing error messages
   Integer, External :: lsq_2dlagrstr  ! Callback for second derivatives structure
   Integer, External :: lsq_2dlagrval  ! Callback for second derivatives values
   Integer, External :: lsq_option     ! Option defining callback routine
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Lsq_ReadMatrix
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Lsq_FDEval
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_Status
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_Solution
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_Message
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_ErrMsg
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Lsq_2DLagrStr
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Lsq_2DLagrVal
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Lsq_Option
#endif
!
!  Control vector
!
   INTEGER, Dimension(:), Pointer :: cntvect
   INTEGER :: coi_error
!
   Integer :: i
   Character(Len=132) :: text
   real*8 :: obj_run1, obj_run2
 
   call startup
!
!  Define data
!
   Call defdata
!
!  Create and initialize a Control Vector
!
   coi_error = coi_create( cntvect )
!
!  Tell CONOPT about the size of the model by populating the Control Vector:
!
   coi_error = max( coi_error, coidef_numvar( cntvect, nobs + dimx ) )
   coi_error = max( coi_error, coidef_numcon( cntvect, nobs + 1 ) )
   coi_error = max( coi_error, coidef_numnz( cntvect, nobs * (dimx+2) ) )
   coi_error = max( coi_error, coidef_numnlnz( cntvect, nobs * (dimx+1) ) )
   coi_error = max( coi_error, coidef_numhess( cntvect, nobs + dimx ) )
   coi_error = max( coi_error, coidef_optdir( cntvect, -1 ) ) !  Minimize
   coi_error = max( coi_error, coidef_objcon( cntvect, nobs + 1 ) )  ! Objective is last constraint
   coi_error = max( coi_error, coidef_optfile( cntvect, 'leastsq9.opt' ) )
!
!  Tell CONOPT about the callback routines:
!
   coi_error = max( coi_error, coidef_readmatrix( cntvect, lsq_readmatrix ) )
   coi_error = max( coi_error, coidef_fdeval( cntvect, lsq_fdeval     ) )
   coi_error = max( coi_error, coidef_status( cntvect, std_status     ) )
   coi_error = max( coi_error, coidef_solution( cntvect, std_solution   ) )
   coi_error = max( coi_error, coidef_message( cntvect, std_message    ) )
   coi_error = max( coi_error, coidef_errmsg( cntvect, std_errmsg     ) )
   coi_error = max( coi_error, coidef_2dlagrstr( cntvect, lsq_2dlagrstr  ) )
   coi_error = max( coi_error, coidef_2dlagrval( cntvect, lsq_2dlagrval  ) )
   coi_error = max( coi_error, coidef_option( cntvect, lsq_option     ) )
 
#if defined(CONOPT_LICENSE_INT_1) && defined(CONOPT_LICENSE_INT_2) && defined(CONOPT_LICENSE_INT_3) && defined(CONOPT_LICENSE_TEXT)
   coi_error = max( coi_error, coidef_license( cntvect, conopt_license_int_1, conopt_license_int_2, conopt_license_int_3, conopt_license_text) )
#endif
 
   If ( coi_error .ne. 0 ) THEN
      write(*,*)
      write(*,*) '**** Fatal Error while loading CONOPT Callback routines.'
      write(*,*)
      call flog( "Skipping Solve due to setup errors", 1 )
   ENDIF
!
!  Save the solution so we can check the duals:
!
   do_allocate = .true.
!
!  Start CONOPT:
!
   run = 1
   coi_error = coi_solve( cntvect )
 
   If ( coi_error /= 0 ) then
      call flog( "Run 1: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Run 1: Status or Solution routine was not called", 1 )
   elseif ( .not. ( sstat == 1 .and. mstat == 2 ) ) then
      call flog( "Run 1: Solver or Model status was not as expected (1,2)", 1 )
!   elseif ( abs( OBJ - 19.44434311d0 ) > 1.d-7 ) then
!      call flog( "Run 1: Incorrect objective returned", 1 )
   Else
      do i  = 1, nobs
         if ( abs( udual(i) - 2*xprim(dimx+i) ) > 2.d-7 ) then
            write(text,*) 'Run 1: Dual and Primal did not match in constraint',i,' P=',xprim(dimx+i),' D=',udual(i)
            call flog( text, 1 )
         endif
      enddo
      Call checkdual( 'Run 1: Leastsq9', minimize )
   endif
   obj_run1 = obj
!
!  Repeat but with bounds on the residuals. They cannot be moved into
!  the post-triangle any more. Also use initial values close to a known
!  alternate solution.
!
   run = 2
   coi_error = coi_solve( cntvect )
 
   If ( coi_error /= 0 ) then
      call flog( "Run 2: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Run 2: Status or Solution routine was not called", 1 )
   elseif ( .not. ( sstat == 1 .and. mstat == 2 ) ) then
      call flog( "Run 2: Solver or Model status was not as expected (1,2)", 1 )
!   elseif ( abs( OBJ - OBJ_Run1 ) > 1.d-7 ) then
!      call flog( "Run 2: Incorrect objective returned", 1 )
   Else
      do i  = 1, nobs
         if ( abs( udual(i) - 2*xprim(dimx+i) ) > 2.d-7 ) then
            write(text,*) 'Run 2: Dual and Primal did not match in constraint',i,' P=',xprim(dimx+i),' D=',udual(i)
            call flog( text, 1 )
         endif
      enddo
      Call checkdual( 'Run 2: Leastsq9', minimize )
   endif
   obj_run2 = obj
!
!  Repeat but with the solution from run2 used as an initial solution.
!
   run = 3
   coi_error = coi_solve( cntvect )
 
   If ( coi_error /= 0 ) then
      call flog( "Run 3: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Run 3: Status or Solution routine was not called", 1 )
   elseif ( .not. ( sstat == 1 .and. mstat == 2 ) ) then
      call flog( "Run 3: Solver or Model status was not as expected (1,2)", 1 )
   elseif ( abs( obj - obj_run2 ) > 1.d-6 ) then !
      call flog( "Run 3: Incorrect objective returned in run3 ", 1 )
   else
      Call checkdual( 'Run 3: Leastsq9', minimize )
   endif
 
   write(*,*)
   write(*,*) 'End of Least Square example 7. Return code=',coi_error
 
   if ( coi_free(cntvect) /= 0 ) call flog( "Error while freeing control vector",1)
 
   call flog( "Successful Solve", 0 )
 
End Program leastsquare
!
! ============================================================================
!  Define information about the model:
!
 
!>  Define information about the model
!!
!! @include{doc} readMatrix_params.dox
Integer Function lsq_readmatrix( lower, curr, upper, vsta, type, rhs, esta,      &
                                 colsta, rowno, value, nlflag, n, m, nz,  &
                                 usrmem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Lsq_ReadMatrix
#endif
   Use proginfo
   Use lsq9
      IMPLICIT NONE
   integer, intent (in)                      :: n      ! number of variables
   integer, intent (in)                      :: m      ! number of constraints
   integer, intent (in)                      :: nz     ! number of nonzeros
   real*8,  intent (in out), dimension(n)    :: lower  ! vector of lower bounds
   real*8,  intent (in out), dimension(n)    :: curr   ! vector of initial values
   real*8,  intent (in out), dimension(n)    :: upper  ! vector of upper bounds
   integer, intent (in out), dimension(n)    :: vsta   ! vector of initial variable status
                                                       ! (not defined here)
   integer, intent (out),    dimension(m)    :: type   ! vector of equation types
   integer, intent (in out), dimension(m)    :: esta   ! vector of initial equation status
                                                       ! (not defined here)
   real*8,  intent (in out), dimension(m)    :: rhs    ! vector of right hand sides
   integer, intent (in out), dimension(n+1)  :: colsta ! vector with start of column indices
   integer, intent (out),    dimension(nz)   :: rowno  ! vector of row numbers
   integer, intent (in out), dimension(nz)   :: nlflag ! vector of nonlinearity flags
   real*8,  intent (in out), dimension(nz)   :: value  ! vector of matrix values
   real*8   usrmem(*)                                  ! optional user memory
 
   Integer :: i, j, k
!
!  Information about Variables:
!  Default: Lower = -Inf, Curr = 0, and Upper = +inf.
!  Default: the status information in Vsta is not used.
!
   do i = 1, dimx
      curr(i) = -0.8d0
   enddo
   if ( run == 2 ) Then  ! Add large bound on Res that will prevent the
      do i = 1, nobs     ! post triangle from being identified.
         lower(dimx+i) = -2.0d0
         upper(dimx+i) = +2.0d0
      enddo
      curr = &
      (/ -1.280997e+00, -3.966357e-01, -1.418473e+00, -7.497624e-01, -1.034839e+00, &
         -1.006460e+00, -1.090813e+00, -1.181425e+00,  2.505639e-01, -1.139540e+00, &
          1.275053e-02, -3.982541e-02,  1.602648e-01, -8.255030e-02, -9.448025e-02, &
         -9.631539e-03,  1.857009e-01, -2.428132e-01,  2.244003e-01, -6.741576e-02, &
         -1.749643e-01, -7.359635e-02,  1.151421e-03,  6.899964e-02 /)
   endif
   if ( run == 3 ) Then    ! In run 3 use the previous solution as the
      curr = xprim
   endif
!
!  Information about Constraints:
!  Default: Rhs = 0
!  Default: the status information in Esta and the function
!           value in FV are not used.
!  Default: Type: There is no default.
!           0 = Equality,
!           1 = Greater than or equal,
!           2 = Less than or equal,
!           3 = Non binding.
!
!  Constraints 1 to Nobs:
!  Rhs = Obs(i) and type Equality
!
   do i = 1, nobs
      rhs(i)  = obs(i)
      type(i) = 0
   enddo
!
!  Constraint Nobs + 1 (Objective)
!  Rhs = 0 and type Non binding
!
   type(nobs+1) = 3
!
!  Information about the Jacobian. CONOPT expects a columnwise
!  representation in Rowno, Value, Nlflag and Colsta.
!
!  Colsta = Start of column indices (No Defaults):
!  Rowno  = Row indices
!  Value  = Value of derivative (by default only linear
!                                derivatives are used)
!  Nlflag = 0 for linear and 1 for nonlinear derivative
!           (not needed for completely linear models)
!
!
!  Indices
!       x(i)  res(j)
!  j:    NL   L=1
!  obj:       NL
!
   k = 1
   do j = 1, dimx
      colsta(j) =  k
      do i = 1, nobs
         rowno(k)  = i
         nlflag(k) = 1
         k         = k + 1
      enddo
   enddo
   do i = 1, nobs
      colsta(dimx+i) = k
      rowno(k)  = i
      nlflag(k) = 0
      value(k)  = 1.d0
      k         = k + 1
      rowno(k)  = nobs+1
      nlflag(k) = 1
      k         = k + 1
   enddo
   colsta(dimx+nobs+1) = k
 
   lsq_readmatrix = 0 ! Return value means OK
 
end Function lsq_readmatrix
!
!==========================================================================
!  Compute nonlinear terms and non-constant Jacobian elements
!
 
!>  Compute nonlinear terms and non-constant Jacobian elements
!!
!! @include{doc} fdeval_params.dox
Integer Function lsq_fdeval( x, g, jac, rowno, jcnm, mode, ignerr, errcnt, &
                             n, nz, thread, usrmem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Lsq_FDEval
#endif
   use lsq9
      IMPLICIT NONE
   integer, intent (in)                      :: n      ! number of variables
   integer, intent (in)                      :: rowno  ! number of the row to be evaluated
   integer, intent (in)                      :: nz     ! number of nonzeros in this row
   real*8,  intent (in),     dimension(n)    :: x      ! vector of current solution values
   real*8,  intent (in out)                  :: g      ! constraint value
   real*8,  intent (in out), dimension(n)    :: jac    ! vector of derivatives for current constraint
   integer, intent (in),     dimension(nz)   :: jcnm   ! list of variables that appear nonlinearly
                                                       ! in this row. Ffor information only.
   integer, intent (in)                      :: mode   ! evaluation mode: 1 = function value
                                                       ! 2 = derivatives, 3 = both
   integer, intent (in)                      :: ignerr ! if 1 then errors can be ignored as long
                                                       ! as errcnt is incremented
   integer, intent (in out)                  :: errcnt ! error counter to be incremented in case
                                                       ! of function evaluation errors.
   integer, intent (in)                      :: thread
   real*8   usrmem(*)                                  ! optional user memory
 
   integer :: i,j
   real*8  :: s
!
!  Row Nobs+1: the objective function is nonlinear
!
   if ( rowno .eq. nobs+1 ) then
!
!  Mode = 1 or 3. Function value: G = sum(i, res(i)**2)
!
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         s    = 0.d0
         do i = 1, nobs
            s = s + x(dimx+i)**2
         enddo
         g    = s
      endif
!
!  Mode = 2 or 3: Derivative values:
!
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         do i = 1, nobs
            jac(dimx+i) = 2*x(dimx+i)
         enddo
      endif
!
!  Row 1 to Nobs: The observation definitions:
!
   else
!
!  Mode = 1 or 3: Function value - x-part only
!
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         s = 0.d0
         do j = 1, dimx
            s = s + a(rowno,j)*x(j) + b(rowno,j)*x(j)**2
         enddo
         g = s
      endif
!
!  Mode = 2 or 3: Derivatives
!
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         do j = 1, dimx
            jac(j) = a(rowno,j) + 2.d0*b(rowno,j)*x(j)
         enddo
      endif
   endif
   lsq_fdeval = 0
 
end Function lsq_fdeval
 
 
!>  Specify the structure of the Lagrangian of the Hessian
!!
!! @include{doc} 2DLagrStr_params.dox
Integer Function lsq_2dlagrstr( HSRW, HSCL, NODRV, N, M, NHESS, UsrMem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Lsq_2DLagrStr
#endif
   use lsq9
   Implicit None
 
   Integer, Intent (IN)                    :: n, m, nhess
   Integer, Intent (IN OUT)                :: nodrv
   Integer, Dimension(Nhess), Intent (Out) :: hsrw, hscl
   real*8,  Intent(IN OUT)                 :: usrmem(*)
 
   integer           :: i
!
!  We assume that this is a Large Residual model, where the most
!  important 2nd derivatives are those that appear directly in the
!  objective.
!  We will therefore only define the 2nd derivatives corresponding to
!  the objective constraint where the 2nd derivatives is 2*I.
!  CONOPT will by default complain that some nonlinear variables do
!  not appear in the Hessian. We will therefore define one dummy element
!  on the diagonal of the Hessian for each of the nonlinear X-variables
!  and give it the numerical value 0.
!
!
!  Define structure of Hessian
!
   do i = 1, dimx+nobs
      hsrw(i) = i
      hscl(i) = i
   enddo
 
   lsq_2dlagrstr = 0
 
End Function lsq_2dlagrstr
 
 
!>  Compute the Lagrangian of the Hessian
!!
!! @include{doc} 2DLagrVal_params.dox
Integer Function lsq_2dlagrval( X, U, HSRW, HSCL, HSVL, NODRV, N, M, NHESS, UsrMem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Lsq_2DLagrVal
#endif
   use lsq9
   Implicit None
 
   Integer, Intent (IN)                    :: n, m, nhess
   Integer, Intent (IN OUT)                :: nodrv
   real*8,  Dimension(N), Intent (IN)      :: x
   real*8,  Dimension(M), Intent (IN)      :: u
   Integer, Dimension(Nhess), Intent (In)  :: hsrw, hscl
   real*8,  Dimension(NHess), Intent (Out) :: hsvl
   real*8,  Intent(IN OUT)                 :: usrmem(*)
 
   integer           :: i
!
!  We assume that this is a Large Residual model, where the most
!  important 2nd derivatives are those that appear directly in the
!  objective.
!  We will therefore only define the 2nd derivatives corresponding to
!  the objective constraint where the 2nd derivatives is 2*I.
!  CONOPT will by default complain that some nonlinear variables do
!  not appear in the Hessian. We will therefore define one dummy element
!  on the diagonal of the Hessian for each of the nonlinear X-variables
!  and give it the numerical value 0.
!
!  Normal Evaluation mode
!
   do i = 1, dimx
      hsvl(i) = 0.d0
   enddo
   do i = 1, nobs
      hsvl(dimx+i) = 2.d0*u(nobs+1)
   enddo
 
   lsq_2dlagrval = 0
 
End Function lsq_2dlagrval
 
 
!>  Sets runtime options
!!
!! @include{doc} option_params.dox
Integer Function lsq_option( ncall, rval, ival, lval, usrmem, name )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Lsq_Option
#endif
   integer  ncall, ival, lval
   character(Len=*) :: name
   real*8   rval
   real*8   usrmem(*)                                  ! optional user memory
 
   Select case (ncall)
   case (1)
      name = 'lkdeb2'  ! Turn the debugger off
      ival = 0
   case default
      name = ' '
   end Select
   lsq_option = 0
 
end Function lsq_option