mp_leastsq14.f90

Go to the documentation of this file.

!> @file mp_leastsq14.f90
!! @ingroup FORTOPENMP_EXAMPLES
!!
!! 
!! This model is an extension of leastsq.f90 in which we test that the
!! solution with 1, 2, and 4 threads give the same solution when
!! threadc is defined as 4.
!! 
!! We solve the following nonlinear least squares model:
!! 
!! min sum(i, res(i)**2 )
!! 
!! sum(j, a(i,j)*x(j) + b(i,j)*x(j)**2 ) + res(i) = obs(i)
!! 
!! where a, b, and obs are known data, and res and x are the
!! variables of the model.
!! 
!!
!! For more information about the individual callbacks, please have a look at the source code.
 
module lsq14
   Integer, parameter :: nobs = 700
   Integer, parameter :: dimx = 500
   real*8, dimension(Nobs,DimX) :: a, b
   real*8, dimension(Nobs)      :: obs
   real*8, dimension(Nobs+DimX) :: xprim, xprim1, xprim2, xprim4
end module lsq14
 
 
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 lsq14
   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 proginfop
   Use coidef
   Use lsq14
   Use omp_lib
   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 :: lsq_solution   ! Special callback for displaying and testing solution values
   Integer, External :: std_message    ! Standard callback for managing messages
   Integer, External :: std_errmsg     ! Standard callback for managing error messages
#if defined(itl)
!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 :: Lsq_Solution
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_Message
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_ErrMsg
#endif
!
!  Control vector
!
   INTEGER, Dimension(:), Pointer :: cntvect
   INTEGER :: coi_error
 
   real*8 time0, time1, time2, time4
   Integer :: i, maxthread
   Logical :: same12, same14, same24
 
   call startup
!
!  Define data
!
   Call defdata
!
!  Create and initialize a Control Vector
!
   coi_error = coi_createfort( 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_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, 'mp_leastsq14.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, lsq_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_debugfv( cntvect, 0 ) )  ! Debug Fdeval on or off
 
#if defined(LICENSE_INT_1) && defined(LICENSE_INT_2) && defined(LICENSE_INT_3) && defined(LICENSE_TEXT)
   coi_error = max( coi_error, coidef_license( cntvect, license_int_1, license_int_2, license_int_3, 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
   maxthread = omp_get_max_threads() ! Define ThreadC as the largest number of threads we can get
   coi_error = max( coi_error, coidef_threadc( cntvect, maxthread ) )
!
!  Start CONOPT -- First time using a single thread:
!
   coi_error = max( coi_error, coidef_threads( cntvect, 1 ) )  ! 1 means use 1 thread
   time0 = omp_get_wtime()
   coi_error = coi_solve( cntvect )
   time1 = omp_get_wtime() - time0
 
   write(*,*)
   write(*,*) 'End of Least Square example 1 with 1 thread. Return code=',coi_error
 
   If ( coi_error /= 0 ) then
      call flog( "One thread: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "One thread: Status or Solution routine was not called", 1 )
   elseif ( .not. ( sstat == 1 .and. mstat == 2 ) ) then
      call flog( "One thread: Solver or Model status was not as expected (1,2)", 1 )
   elseif ( abs( obj - 19.44434311d0 ) > 1.d-7 ) then
      call flog( "One thread: Incorrect objective returned", 1 )
   endif
   xprim1 = xprim
!
!  Start CONOPT again -- this time using 2 threads:
!
   coi_error = max( coi_error, coidef_threads( cntvect, 2 ) )  ! 2 means use 2 threads
   time0 = omp_get_wtime()
   coi_error = coi_solve( cntvect )
   time2 = omp_get_wtime() - time0
 
   write(*,*)
   write(*,*) 'Two threads: End of Least Square example 1. Return code=',coi_error
 
   If ( coi_error /= 0 ) then
      call flog( "Two threads: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Two threads: Status or Solution routine was not called", 1 )
   elseif ( .not. ( sstat == 1 .and. mstat == 2 ) ) then
      call flog( "Two threads: Solver or Model status was not as expected (1,2)", 1 )
   elseif ( abs( obj - 19.44434311d0 ) > 1.d-7 ) then
      call flog( "Two threads: Incorrect objective returned", 1 )
   endif
   xprim2 = xprim
   if ( maxthread >= 4 ) then
!
!  Start CONOPT again -- this time using 4 threads (if available on this machine):
!
   coi_error = max( coi_error, coidef_threads( cntvect, 4 ) )  ! 4 means use 4 threads
   time0 = omp_get_wtime()
   coi_error = coi_solve( cntvect )
   time4 = omp_get_wtime() - time0
 
   write(*,*)
   write(*,*) 'Four threads: End of Least Square example 1. Return code=',coi_error
 
   If ( coi_error /= 0 ) then
      call flog( "Four threads: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Four threads: Status or Solution routine was not called", 1 )
   elseif ( .not. ( sstat == 1 .and. mstat == 2 ) ) then
      call flog( "Four threads: Solver or Model status was not as expected (1,2)", 1 )
   elseif ( abs( obj - 19.44434311d0 ) > 1.d-7 ) then
      call flog( "Four threads: Incorrect objective returned", 1 )
   endif
   xprim4 = xprim
   Endif
 
   same12 = .true.
   same14 = .true.
   same24 = .true.
   Do i = 1, nobs+dimx
      If ( xprim2(i) /= xprim1(i) ) then
         if ( same12 ) then
            write(*,*) 'Solution 2 is different from Solution 1. First diff for index',i
         endif
         same12 = .false.
      endif
      if ( maxthread >= 4 ) then
         If ( xprim4(i) /= xprim1(i) ) then
            if ( same14 ) then
               write(*,*) 'Solution 4 is different from Solution 1. First diff for index',i
            endif
            same14 = .false.
         endif
         If ( xprim4(i) /= xprim2(i) ) then
            if ( same24 ) then
               write(*,*) 'Solution 4 is different from Solution 2. First diff for index',i
            endif
            same24 = .false.
         endif
      endif
   Enddo
   write(*,*)
   write(*,"('Time for single thread ',f10.3)") time1
   write(*,"('Time for two    threads',f10.3)") time2
   if ( maxthread >= 4 ) &
   write(*,"('Time for four   threads',f10.3)") time4
   write(*,"('Speedup    two  threads',f10.3)") time1/time2
   write(*,"('Efficiency two  threads',f10.3)") time1/time2/2
   if ( maxthread >= 4 ) &
   write(*,"('Speedup    four threads',f10.3)") time1/time4
   if ( maxthread >= 4 ) &
   write(*,"('Efficiency four threads',f10.3)") time1/time4/4
 
#if defined (nocomp)
#else
   if ( same12 .and. same14 .and. same24 ) then
      write(*,*) 'All solutions are the same'
   else
      write(*,*) 'The solutions are NOT the same'
      call flog( "Solutions are not the same", 1 )
   Endif
#endif
 
   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 )
#if defined(itl)
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Lsq_ReadMatrix
#endif
   Use lsq14
   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
!
!  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. We use the standard method with
!  Rowno, Value, Nlflag and Colsta and we do not use Colno.
!
!  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 )
#if defined(itl)
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Lsq_FDEval
#endif
   use lsq14
   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
 
Integer Function lsq_solution( XVAL, XMAR, XBAS, XSTA, YVAL, YMAR, YBAS, YSTA, N, M, USRMEM )
#if defined(itl)
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Lsq_Solution
#endif
!
!  Simple implementation in which we write the solution values to
!  the 'Documentation file' on unit 10.
!
   Use proginfop
   Use lsq14
   IMPLICIT NONE
   INTEGER, Intent(IN)                  :: n, m
   INTEGER, Intent(IN),    Dimension(N) :: xbas, xsta
   INTEGER, Intent(IN),    Dimension(M) :: ybas, ysta
   real*8,  Intent(IN),    Dimension(N) :: xval, xmar
   real*8,  Intent(IN),    Dimension(M) :: yval, ymar
   real*8,  Intent(IN OUT)              :: usrmem(*)
 
   INTEGER i
   CHARACTER*5, Parameter, Dimension(4) :: stat = (/ 'Lower','Upper','Basic','Super' /)
 
   WRITE(10,"(/' Variable    Solution value    Reduced cost    B-stat')")
   DO i = 1, n
      WRITE(10,"(1P,I7,E20.6,E16.6,4X,A5 )") i, xval(i), xmar(i), stat(1+xbas(i))
   ENDDO
 
   WRITE(10,"(/' Constrnt    Activity level    Marginal cost   B-stat')")
   DO i = 1, m
      WRITE(10,"(1P,I7,E20.6,E16.6,4X,A5 )") i, yval(i), ymar(i), stat(1+ybas(i))
   ENDDO
   xprim(1:n) = xval(1:n)
#if defined (wei)
   Call flush(10)
#endif
 
   solcalls = solcalls + 1
   lsq_solution = 0
 
END Function lsq_solution