program fortrantest
use, intrinsic :: iso_c_binding
use highs_fortran_api
implicit none
! This illustrates the use of Highs_lpCall, the simple F90 interface to
! HiGHS. It's designed to solve the general LP problem
!
! Min c^Tx subject to L <= Ax <= U; l <= x <= u
!
! where A is a matrix with m rows and n columns
!
! The scalar n is numcol
! The scalar m is numrow
!
! The vector c is colcost
! The vector l is collower
! The vector u is colupper
! The vector L is rowlower
! The vector U is rowupper
!
! The matrix A is represented in packed column-wise form: only its
! nonzeros are stored
!
! * The number of nonzeros in A is numnz
!
! * The row indices of the nonnzeros in A are stored column-by-column
! in aindex
!
! * The values of the nonnzeros in A are stored column-by-column in
! avalue
!
! * The position in aindex/avalue of the index/value of the first
! nonzero in each column is stored in astart
!
! Note that astart(1) must be zero
!
! After a successful call to Highs_lpCall, the primal and dual
! solution, and the simplex basis are returned as follows
!
! The vector x is colvalue
! The vector Ax is rowvalue
! The vector of dual values for the variables x is coldual
! The vector of dual values for the variables Ax is rowdual
! The basic/nonbasic status of the variables x is colbasisstatus
! The basic/nonbasic status of the variables Ax is rowbasisstatus
!
! The status of the solution obtained is modelstatus
!
! To solve maximization problems, the values in c must be negated
!
! The use of Highs_lpCall is illustrated for the LP/MIP example
!
! Min f = 2x_0 + 3x_1
! s.t. x_1 <= 6
! 10 <= x_0 + 2x_1 <= 14
! 8 <= 2x_0 + x_1
! 0 <= x_0 <= 3; 1 <= x_1
integer ( c_int ), parameter :: numcol = 2
integer ( c_int ), parameter :: numrow = 3
integer ( c_int ), parameter :: numnz = 5
integer ( c_int ), parameter :: aformat_colwise = 1
integer ( c_int ), parameter :: sense = 1
real ( c_double ), parameter :: offset = 0
real ( c_double ) colcost(numcol)
real ( c_double ) collower(numcol)
real ( c_double ) colupper(numcol)
real ( c_double ) rowlower(numrow)
real ( c_double ) rowupper(numrow)
integer ( c_int ) astart(numcol)
integer ( c_int ) aindex(numnz)
real ( c_double ) avalue(numnz)
real ( c_double ) colvalue(numcol)
real ( c_double ) coldual(numcol)
real ( c_double ) rowvalue(numrow)
real ( c_double ) rowdual(numrow)
integer ( c_int ) colbasisstatus(numcol)
integer ( c_int ) rowbasisstatus(numrow)
integer modelstatus
integer runstatus
integer ( c_int ), parameter :: modelstatus_optimal = 7
integer ( c_int ), parameter :: runstatus_error = -1
integer ( c_int ), parameter :: runstatus_ok = 0
integer ( c_int ), parameter :: runstatus_warning = 1
! For the full API test
type ( c_ptr ) :: highs
! For the row-wise matrix
integer ( c_int ) arstart(numrow)
integer ( c_int ) arindex(numnz)
real ( c_double ) arvalue(numnz)
real, parameter :: inf = 1e30
integer col, row, el
integer from_el, to_el
integer iteration_count, primal_solution_status, dual_solution_status
double precision objective_function_value
double precision objective_error
integer option_type
integer dummy_info
integer alt_sense
integer scale_strategy
integer, parameter :: default_scale_strategy = 1
integer, parameter :: new_scale_strategy = 3
double precision, parameter :: dual_tolerance = 1d-6
logical ( c_bool ) write_solution_to_file
integer ( c_int ) write_solution_style
double precision, pointer :: double_null(:)
integer, pointer :: integer_null(:)
character( c_char ) :: file_name(7)
! Illustrate the solution of a QP
!
! minimize -x_2 - 3x_3 + (1/2)(2x_1^2 - 2x_1x_3 + 0.2x_2^2 + 2x_3^2)
!
! subject to x_1 + x_3 <= 2; x>=0
!
! Solution x_1 = 0.5; x_2 = 5.0; x_3 = 1.5
integer ( c_int ), parameter :: qp_numcol = 3
integer ( c_int ), parameter :: qp_numrow = 1
integer ( c_int ), parameter :: qp_numnz = 2
integer ( c_int ), parameter :: qp_hessian_numnz = 4
integer ( c_int ), parameter :: qformat_triangular = 1
real ( c_double ) qp_colcost(qp_numcol)
real ( c_double ) qp_collower(qp_numcol)
real ( c_double ) qp_colupper(qp_numcol)
real ( c_double ) qp_rowlower(qp_numrow)
real ( c_double ) qp_rowupper(qp_numrow)
integer ( c_int ) qp_astart(qp_numcol)
integer ( c_int ) qp_aindex(qp_numnz)
real ( c_double ) qp_avalue(qp_numnz)
integer ( c_int ) qp_qstart(qp_numcol)
integer ( c_int ) qp_qindex(qp_hessian_numnz)
real ( c_double ) qp_qvalue(qp_hessian_numnz)
real ( c_double ) qp_sol(qp_numcol)
real ( c_double ) qp_colvalue(qp_numcol)
real ( c_double ) qp_coldual(qp_numcol)
real ( c_double ) qp_rowvalue(qp_numrow)
real ( c_double ) qp_rowdual(qp_numrow)
integer ( c_int ) qp_colbasisstatus(qp_numcol)
integer ( c_int ) qp_rowbasisstatus(qp_numrow)
integer qp_modelstatus
integer qp_runstatus
logical, parameter :: no_highs_logging = .TRUE.
logical ( c_bool ), parameter :: logical_false = .false.
logical ( c_bool ), parameter :: logical_true = .true.
! Set up the LP/MIP example
colcost(1) = 2
colcost(2) = 3
collower(1) = 0
collower(2) = 1
colupper(1) = 3
colupper(2) = inf
rowlower(1) = -inf
rowlower(2) = 10.0
rowlower(3) = 8.0
rowupper(1) = 6.0
rowupper(2) = 14.0
rowupper(3) = inf
! The definition of sparse matrices to be passed into the FORTRAN
! interface is non-trivial, since the FORTRAN interface is a direct
! bind to the C API: no execution (esp. conversion of data) at any
! point. Hence FORTRAN users have to supply vectors that, when
! indexed from 0, are standard for the C API.
! Although the FORTRAN arrays in the example are indexed from 1 (of
! course) note that the row indices in aindex are in the interval
! [0, numrow-1]
! For this example, aindex is indexed from 1 to 5, and contains [1 2
! 0 1 2], referring to column 0 having entries in rows 1 and 2;
! column 1 having entries in rows 0, 1 and 2
! FORTRAN-wise, astart would naturally be [1 3], but it must
! indicate the starts when aindex is indexed from 0 to 4 in C, so
! the starts must have 1 subtracted from them from the values that
! would be used in FORTRAN. Hence astart is [0 2]
astart(1) = 0
astart(2) = 2
aindex(1) = 1
aindex(2) = 2
aindex(3) = 0
aindex(4) = 1
aindex(5) = 2
avalue(1) = 1
avalue(2) = 2
avalue(3) = 1
avalue(4) = 2
avalue(5) = 1
!================================================================================
! Illustrate use of Highs_lpCall to solve a given LP
print*, "*********"
print*, "Section 1"
print*, "*********"
runstatus = Highs_lpCall( numcol, numrow, numnz,&
aformat_colwise, sense, offset, &
colcost, collower, colupper, rowlower, rowupper,&
astart, aindex, avalue,&
colvalue, coldual, rowvalue, rowdual,&
colbasisstatus, rowbasisstatus, modelstatus)
if (runstatus .ne. runstatus_ok) then
write(*, '(a, i1, a, i2)')'Highs_lpCall run status is ', runstatus, ' not ', runstatus_ok
stop
endif
write(*, '(a, i1, a, i2)')'Run status = ', runstatus, '; Model status = ', modelstatus
if (modelstatus .eq. modelstatus_optimal) then
objective_function_value = 0
! Report the column primal and dual values, and basis status
do col = 1, numcol
write(*, '(a, i1, a, f10.4, a, f10.4, a, i2)') &
'Col', col, ' = ', colvalue(col), &
'; dual = ', coldual(col), &
'; status = ', colbasisstatus(col)
objective_function_value = objective_function_value + colvalue(col)*colcost(col)
enddo
! Report the row primal and dual values, and basis status
do row = 1, numrow
write(*, '(a, i1, a, f10.4, a, f10.4, a, i2)') &
'Row', row, ' = ', rowvalue(row), &
'; dual = ', rowdual(row), &
'; status = ', rowbasisstatus(row)
enddo
write(*, '(a, f10.4)')'Optimal objective value = ', objective_function_value
endif
! Illustrate use of Highs_create() to create a pointer to an
! instance of the Highs class, then Highs_passLp to pass LP to
! HiGHS, and Highs_run to solve it
highs = Highs_create()
if (no_highs_logging) then
runstatus = Highs_setBoolOptionValue(highs, "output_flag"//C_NULL_CHAR, logical_false)
endif
runstatus = Highs_passLp(highs, numcol, numrow, numnz, aformat_colwise, &
sense, offset, colcost, collower, colupper, rowlower, rowupper, &
astart, aindex, avalue)
runstatus = Highs_run(highs)
modelstatus = Highs_getModelStatus(highs);
print*, "modelstatus = ", modelstatus
call assert(runstatus .eq. runstatus_ok, "Highs_run runstatus")
call assert(modelstatus .eq. modelstatus_optimal, "Highs_run modelstatus optimal")
write(*, '(a, i1, a, i2)')'Run status = ', runstatus, '; Model status = ', modelstatus
runstatus = Highs_getDoubleInfoValue(highs, "objective_function_value"//C_NULL_CHAR, objective_function_value);
runstatus = Highs_getIntInfoValue(highs, "simplex_iteration_count"//C_NULL_CHAR, iteration_count);
write(*, '(a, f10.4, a, i6)')"Objective value = ", objective_function_value, "; Iteration count = ", iteration_count
call Highs_destroy(highs)
!================================================================================
! Illustrate use of Highs_addCols and Highs_addRows to build model,
! and then Highs_changeObjectiveSense to switch to maximization
print*, "*********"
print*, "Section 2"
print*, "*********"
highs = Highs_create()
if (no_highs_logging) then
runstatus = Highs_setBoolOptionValue(highs, "output_flag"//C_NULL_CHAR, logical_false)
endif
! Create double and integer values equal to NULL pointer
call C_F_POINTER(C_NULL_PTR, double_null, [0])
call C_F_POINTER(C_NULL_PTR, integer_null, [0])
! Add two columns to the empty LP, but no matrix. After numnz=0, can
! just pass arrays rather than NULL
runstatus = Highs_addCols(highs, numcol, colcost, collower, colupper, 0, integer_null, integer_null, double_null);
! Define the constraint matrix by adding it as three rows to the
! 2-column LP - requiring the matrix row-wise
arstart(1) = 0
arstart(2) = 1
arstart(3) = 3
arindex(1) = 1
arindex(2) = 0
arindex(3) = 1
arindex(4) = 0
arindex(5) = 1
arvalue(1) = 1
arvalue(2) = 1
arvalue(3) = 2
arvalue(4) = 2
arvalue(5) = 1
runstatus = Highs_addRows(highs, numrow, rowlower, rowupper, numnz, arstart, arindex, arvalue)
runstatus = Highs_getObjectiveSense(highs, alt_sense);
write(*, '(a, i2)')"LP problem has objective sense = ", alt_sense
call assert(alt_sense .eq. sense, "Objective sense")
alt_sense = -1 * alt_sense
runstatus = Highs_changeObjectiveSense(highs, alt_sense);
runstatus = Highs_getObjectiveSense(highs, alt_sense);
call assert(alt_sense .eq. -1, "Changed Objective sense")
! Get and set option values
runstatus = Highs_getIntOptionValue(highs, "simplex_scale_strategy"//C_NULL_CHAR, scale_strategy);
call assert(scale_strategy .eq. default_scale_strategy,&
"scale_strategy .eq. default_scale_strategy")
runstatus = Highs_setIntOptionValue(highs, "simplex_scale_strategy"//C_NULL_CHAR, new_scale_strategy)
runstatus = Highs_getIntOptionValue(highs, "simplex_scale_strategy"//C_NULL_CHAR, scale_strategy);
call assert(scale_strategy .eq. new_scale_strategy,&
"scale_strategy .eq. new_scale_strategy")
runstatus = Highs_setDoubleOptionValue(highs, "primal_feasibility_tolerance"//C_NULL_CHAR, 1d-6);
call assert(runstatus .eq. runstatus_ok, "setDoubleOptionValue runstatus")
runstatus = Highs_setDoubleOptionValue(highs, "dual_feasibility_tolerance"//C_NULL_CHAR, dual_tolerance);
call assert(runstatus .eq. runstatus_ok, "setDoubleOptionValue runstatus")
! There are some functions to check what type of option value you should provide.
runstatus = Highs_getOptionType(highs, "simplex_scale_strategy"//C_NULL_CHAR, option_type);
call assert(runstatus .eq. runstatus_ok, "getOptionType runstatus = 0")
call assert(option_type .eq. 1, "getOptionType option_type")
! This is what happens if an invalid name is passed
runstatus = Highs_getOptionType(highs, "bad_option"//C_NULL_CHAR, option_type)
call assert(runstatus .eq. runstatus_error, "getOptionType runstatus")
! Suppress HiGHS output
print*, "Suppressing all HiGHS output"
runstatus = Highs_setBoolOptionValue(highs, "output_flag"//C_NULL_CHAR, logical_false)
! Solve the LP
runstatus = Highs_run(highs);
! Get the model status
modelstatus = Highs_getModelStatus(highs);
write(*, '(a, i1, a, i2)')'Run status = ', runstatus, '; Model status = ', modelstatus
! Get solution data
runstatus = Highs_getDoubleInfoValue(highs, "objective_function_value"//C_NULL_CHAR, objective_function_value);
runstatus = Highs_getIntInfoValue(highs, "simplex_iteration_count"//C_NULL_CHAR, iteration_count);
runstatus = Highs_getIntInfoValue(highs, "primal_solution_status"//C_NULL_CHAR, primal_solution_status);
runstatus = Highs_getIntInfoValue(highs, "dual_solution_status"//C_NULL_CHAR, dual_solution_status);
! This is what happens if an invalid name is passed
runstatus = Highs_getIntInfoValue(highs, "bad_info"//C_NULL_CHAR, dummy_info)
call assert(runstatus .eq. runstatus_error, "getOptionType runstatus")
write(*, '(a, f10.4, a, i6)')"Objective value = ", objective_function_value, "; Iteration count = ", iteration_count
print*, "modelstatus = ", modelstatus
call assert(modelstatus .eq. modelstatus_optimal, "Optimal => modelstatus = modelstatus_optimal")
if (modelstatus .eq. modelstatus_optimal) then
call assert(primal_solution_status .eq. 2, "Optimal => primal_solution_status = 2")
call assert(dual_solution_status .eq. 2, "Optimal => dual_solution_status = 2")
! Get the primal and dual solution
runstatus = Highs_getSolution(highs, colvalue, coldual, rowvalue, rowdual);
! Get the basis
runstatus = Highs_getBasis(highs, colbasisstatus, rowbasisstatus);
! Report the column primal and dual values, and basis status
do col = 1, numcol
write(*, '(a, i1, a, f10.4, a, f10.4, a, i2)') &
'Col', col, ' = ', colvalue(col), &
'; dual = ', coldual(col), &
'; status = ', colbasisstatus(col)
enddo
! Report the row primal and dual values, and basis status
do row = 1, numrow
write(*, '(a, i1, a, f10.4, a, f10.4, a, i2)') &
'Row', row, ' = ', rowvalue(row), &
'; dual = ', rowdual(row), &
'; status = ', rowbasisstatus(row)
enddo
endif
! Write out model as MPS for use later
runstatus = Highs_writeModel(highs, "F90.mps"//C_NULL_CHAR)
! runstatus is runstatus_warning since there are no names
call assert(runstatus .eq. runstatus_warning, "Highs_writeModel runstatus")
call Highs_destroy(highs)
!================================================================================
! Illustrate use of Highs_readModel to read model, and
! Highs_writeSolution to write the solution
print*, "*********"
print*, "Section 3"
print*, "*********"
highs = Highs_create()
if (no_highs_logging) then
runstatus = Highs_setBoolOptionValue(highs, "output_flag"//C_NULL_CHAR, logical_false)
endif
! Read the LP
runstatus = Highs_readModel(highs, "F90.mps"//C_NULL_CHAR)
call assert(runstatus .eq. runstatus_ok, "Highs_readModel runstatus")
! Solve the LP
runstatus = Highs_run(highs);
! Get the model status
modelstatus = Highs_getModelStatus(highs);
write(*, '(a, i1, a, i2)')'Run status = ', runstatus, '; Model status = ', modelstatus
runstatus = Highs_getDoubleInfoValue(highs, "objective_function_value"//C_NULL_CHAR, objective_function_value);
runstatus = Highs_getIntInfoValue(highs, "simplex_iteration_count"//C_NULL_CHAR, iteration_count);
write(*, '(a, f10.4, a, i6)')"Objective value = ", objective_function_value, "; Iteration count = ", iteration_count
! Write solution to the screen
runstatus = Highs_writeSolutionPretty(highs, ""//C_NULL_CHAR)
call Highs_destroy(highs)
!================================================================================
! Illustrate use of setting bool options and string options
! (solution_file) so only run(highs) is required
print*, "*********"
print*, "Section 4"
print*, "*********"
highs = Highs_create()
if (no_highs_logging) then
runstatus = Highs_setBoolOptionValue(highs, "output_flag"//C_NULL_CHAR, logical_false)
endif
! Get and set string options
runstatus = Highs_getStringOptionValue(highs, "solution_file"//C_NULL_CHAR, file_name)
print*, "Default solution_file is |", file_name, "|"
runstatus = Highs_setStringOptionValue(highs, "solution_file"//C_NULL_CHAR, "F90.sol"//C_NULL_CHAR)
runstatus = Highs_getStringOptionValue(highs, "solution_file"//C_NULL_CHAR, file_name)
print*, "New solution_file is |", file_name, "|"
! Get and set bool options. NB Cannot pass .true. as it's 4-byte
runstatus = Highs_getBoolOptionValue(highs, "write_solution_to_file"//C_NULL_CHAR, write_solution_to_file)
print*, "Default write_solution_to_file = ", write_solution_to_file
write_solution_to_file = .true.
runstatus = Highs_setBoolOptionValue(highs, "write_solution_to_file"//C_NULL_CHAR, write_solution_to_file)
runstatus = Highs_getIntOptionValue(highs, "write_solution_style"//C_NULL_CHAR, write_solution_style)
print*, "Default write_solution_style = ", write_solution_style
write_solution_style = 1;
runstatus = Highs_setIntOptionValue(highs, "write_solution_style"//C_NULL_CHAR, write_solution_style)
! Report all the deviations from default options
runstatus = Highs_writeOptionsDeviations(highs, "OptionsDeviations.set"//C_NULL_CHAR)
! Reset all the options
runstatus = Highs_resetOptions(highs)
! Report all the options
runstatus = Highs_writeOptions(highs, "Options.set"//C_NULL_CHAR)
if (no_highs_logging) then
runstatus = Highs_setBoolOptionValue(highs, "output_flag"//C_NULL_CHAR, logical_false)
endif
! Read the LP
runstatus = Highs_readModel(highs, "F90.mps"//C_NULL_CHAR)
! Solve the LP
runstatus = Highs_run(highs);
! Get the model status
modelstatus = Highs_getModelStatus(highs);
call assert(modelstatus .eq. modelstatus_optimal, "Highs_run modelstatus optimal")
write(*, '(a, i1, a, i2)')'Run status = ', runstatus, '; Model status = ', modelstatus
runstatus = Highs_getDoubleInfoValue(highs, "objective_function_value"//C_NULL_CHAR, objective_function_value);
runstatus = Highs_getIntInfoValue(highs, "simplex_iteration_count"//C_NULL_CHAR, iteration_count);
write(*, '(a, f10.4, a, i6)')"Objective value = ", objective_function_value, "; Iteration count = ", iteration_count
call Highs_destroy(highs)
!================================================================================
! Illustrate use of Highs_qpCall to solve a given QP
print*, "**********"
print*, "QP Example"
print*, "**********"
qp_colcost(1) = 0
qp_colcost(2) = -1
qp_colcost(3) = -3
qp_collower(1) = 0
qp_collower(2) = 0
qp_collower(3) = 0
qp_colupper(1) = inf
qp_colupper(2) = inf
qp_colupper(3) = inf
qp_rowlower(1) = -inf
qp_rowupper(1) = 2
qp_astart(1) = 0
qp_astart(2) = 1
qp_astart(3) = 1
qp_aindex(1) = 0
qp_aindex(2) = 0
qp_avalue(1) = 1
qp_avalue(2) = 1
qp_qstart(1) = 0
qp_qstart(2) = 2
qp_qstart(3) = 3
qp_qindex(1) = 0
qp_qindex(2) = 2
qp_qindex(3) = 1
qp_qindex(4) = 2
qp_qvalue(1) = 2.0
qp_qvalue(2) = -1.0
qp_qvalue(3) = 0.2
qp_qvalue(4) = 2.0
qp_sol(1) = 0.5
qp_sol(2) = 5.0
qp_sol(3) = 1.5
runstatus = Highs_qpCall( qp_numcol, qp_numrow, qp_numnz, qp_hessian_numnz,&
aformat_colwise, qformat_triangular, sense, offset,&
qp_colcost, qp_collower, qp_colupper, qp_rowlower, qp_rowupper,&
qp_astart, qp_aindex, qp_avalue,&
qp_qstart, qp_qindex, qp_qvalue,&
qp_colvalue, qp_coldual, qp_rowvalue, qp_rowdual,&
qp_colbasisstatus, qp_rowbasisstatus, modelstatus)
if (runstatus .ne. runstatus_ok) then
write(*, '(a, i1, a, i2)')'Highs_lpCall run status is ', runstatus, ' not ', runstatus_ok
stop
endif
write(*, '(a, i1, a, i2)')'Run status = ', runstatus, '; Model status = ', modelstatus
if (modelstatus .eq. modelstatus_optimal) then
objective_function_value = 0
! Report the column primal and dual values, and basis status
do col = 1, qp_numcol
write(*, '(a, i1, a, f10.4, a, f10.4, a, i2)') &
'Col', col, ' = ', qp_colvalue(col), &
'; dual = ', qp_coldual(col)
call assert(abs(qp_colvalue(col)-qp_sol(col)) .le. 1e-4, "Solution error")
enddo
! Report the row primal and dual values, and basis status
do row = 1, qp_numrow
write(*, '(a, i1, a, f10.4, a, f10.4, a, i2)') &
'Row', row, ' = ', qp_rowvalue(row), &
'; dual = ', qp_rowdual(row)
enddo
do col = 1, qp_numcol
objective_function_value = objective_function_value + qp_colvalue(col)*qp_colcost(col)
enddo
do col = 1, qp_numcol
from_el = qp_qstart(col)
if (col < qp_numcol) then
to_el = qp_qstart(col+1)-1
else
to_el = qp_hessian_numnz-1
endif
objective_function_value = &
objective_function_value + 0.5 * qp_colvalue(col) * qp_qvalue(from_el+1) * qp_colvalue(col)
do el = from_el+1, to_el
row = qp_qindex(el+1)+1
objective_function_value = &
objective_function_value + qp_colvalue(col) * qp_qvalue(el+1) * qp_colvalue(row)
enddo
enddo
write(*, '(a, f10.4)')'Optimal objective value = ', objective_function_value
objective_error = abs(objective_function_value+5.25)
call assert(objective_error .le. 1e-4, "Objective error")
endif
end program fortrantest
subroutine assert ( logic, message)
logical logic
character*(*) message
if (.not.logic) then
write(*, '(a, a)')'assert fail for ', message
stop
endif
end subroutine assert