conmax_module Module

Module for CONMAX.

Notes

  • We note here that the only active write statements in this package are in the sample driver program, but some of those which have been commented out (along with the statements nwrit=output_unit) elsewhere in the program could prove useful, especially the statement write(nwrit,1100)... in subroutine conmax.)

Reference

  • E. H. Kaufman Jr., D. J. Leeming & G. D. Taylor, "An ODE-based approach to nonlinearly constrained minimax problems", Numerical Algorithms, Volume 9, pages 25-37 (1995). Note that since that paper was written, tolcon was deleted from the argument list of conmax.

History

  • CONMAX, Version 1.7. This package was last revised on December 4, 1996. from: http://www.netlib.org/opt/conmax.f
  • Jacob Williams, August 2021 : Refactored into modern Fortran.

Uses

  • module~~conmax_module~~UsesGraph module~conmax_module conmax_module iso_fortran_env iso_fortran_env module~conmax_module->iso_fortran_env

Variables

Type Visibility Attributes Name Initial
integer, private, parameter :: nwrit = output_unit
real(kind=wp), private, parameter :: zero = 0.0_wp
real(kind=wp), private, parameter :: one = 1.0_wp
real(kind=wp), private, parameter :: two = 2.0_wp
real(kind=wp), private, parameter :: three = 3.0_wp
real(kind=wp), private, parameter :: four = 4.0_wp
real(kind=wp), private, parameter :: ten = 10.0_wp
real(kind=wp), private, parameter :: spcmn = real(radix(1.0_wp), wp)**(-digits(1.0_wp))

d1mach3: the smallest relative spacing

real(kind=wp), private, parameter :: big = one/spcmn

Abstract Interfaces

abstract interface

  • private subroutine func(me, Nparm, Numgr, Pttbl, Iptb, Indm, Param, Ipt, Indfn, Icntyp, Confun)

    Interface for the fnset function.

    The first eight variables in the calling sequence for fnset are for input to fnset, with the first five variables being exactly as the user set them in the driver program. if the ten thousands digit of ioptn was set to 0, fnset should be written to place the appropriate values in icntyp and confun using the parameters in param, as follows:

    • icntyp(ipt) = the type of the ipt(th) constraint (i.e. 2, 1, -1, or -2), or the user can set icntyp(ipt)=0 as a signal to ignore constraint ipt.

    • confun(ipt,1) = the appropriate value as discussed above. (this can be left undefined if icntyp(ipt)=0.)

    if indfn=1 (which is the only possibility other than indfn=0) then in addition to the above, for j=1,...,nparm, fnset should compute

    confun(ipt,j+1) = the value of the partial derivative with respect to param(j) of the function whose value was computed in confun(ipt,1) (unless icntyp(ipt)=0, in which case these values need not be computed).

    Arguments

    Type IntentOptional Attributes Name
    class(conmax_solver), intent(inout) :: me
    integer, intent(in) :: Nparm
    integer, intent(in) :: Numgr
    real(kind=wp), intent(in), dimension(Iptb, Indm) :: Pttbl
    integer, intent(in) :: Iptb
    integer, intent(in) :: Indm
    real(kind=wp), intent(in), dimension(Nparm) :: Param
    integer, intent(in) :: Ipt
    integer, intent(in) :: Indfn
    integer, intent(out), dimension(Numgr) :: Icntyp
    real(kind=wp), intent(out), dimension(Numgr,Nparm+1) :: Confun

Derived Types

type, public, abstract ::  conmax_solver

main CONMAX solver class. This class must be extended to define the users's fnset function.

Type-Bound Procedures

procedure, public :: solve => conmax ../../

main solver routine

procedure, public :: muller
procedure, public :: searsl
procedure(func), public, deferred :: fnset
procedure, private :: ercmp1
procedure, private :: rkcon
procedure, private :: slpcon
procedure, private :: corrct
procedure, private :: rkpar
procedure, private :: searcr
procedure, private :: derst
procedure, private :: setu1
procedure, private :: pmtst

Functions

private pure function iloc(Iarr, Nparm, Numgr)

This function subprogram returns the subscript of the first element of array iarr relative to iwork (if the array is integer, i.e. 13 <= iarr <= 23) or relative to work (if the array is floating point, i.e. 1 <= iarr <= 12 or 24 <= iarr <= 48).

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: Iarr
integer, intent(in) :: Nparm
integer, intent(in) :: Numgr

Return Value integer

private pure function dotprd(Lgth, Vec1, Vec2, Nparm) result(dd)

This subprogram computes the dot product of vectors vec1 and vec2 of length lgth. vec1 and vec2 do not appear in function iloc since they are used only as input names for this subprogram, and so they don't need to have space reserved for them in the array work.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: Lgth
real(kind=wp), intent(in), dimension(Nparm + 1) :: Vec1
real(kind=wp), intent(in), dimension(Nparm + 1) :: Vec2
integer, intent(in) :: Nparm

Return Value real(kind=wp)


Subroutines

private subroutine conmax(me, ioptn, nparm, numgr, itlim, fun, ifun, pttbl, iptb, Indm, iwork, liwrk, work, lwrk, iter, param, error)

CONMAX consists of two programs for solving the problem

Read more…

Arguments

Type IntentOptional Attributes Name
class(conmax_solver), intent(inout) :: me
integer, intent(in) :: ioptn

THIS IS THE OPTION SWITCH, WHICH SHOULD BE SET TO 0 UNLESS ONE OR MORE OF THE EXTRA OPTIONS DESCRIBED BELOW IS USED. THE USER HAS SEVERAL EXTRA OPTIONS WHICH ARE ACTIVATED BY SETTING IOPTN TO A VALUE OTHER THAN 0; MORE THAN ONE AT A TIME CAN BE USED. IN PARTICULAR:

Read more…
integer, intent(in) :: nparm

THIS IS THE NUMBER OF PARAMETERS IN THE PROBLEM. (THEY ARE STORED IN PARAM -- SEE BELOW.)

integer, intent(in) :: numgr

THIS IS THE TOTAL NUMBER OF CONSTRAINTS.

integer, intent(in) :: itlim

THIS IS THE LIMIT ON THE NUMBER OF ITERATIONS, I.E. THE LIMIT ON THE NUMBER OF TIMES THE PROGRAM REDUCES W. IF ITLIM IS SET TO 0 THE PROGRAM WILL COMPUTE THE ERRORS FOR THE INITIAL APPROXIMATION AND STOP WITHOUT CHECKING FEASIBILITY.

real(kind=wp), intent(in) :: fun(Ifun)

(VECTOR ARRAY OF DIMENSION IFUN) THIS IS A VECTOR ARRAY IN WHICH DATA OR FUNCTION VALUES IN TYPE 2 CONSTRAINTS (SEE ABOVE) CAN BE STORED. FUN(I) NEED NOT BE ASSIGNED A VALUE IF IT IS NOT GOING TO BE USED.

integer, intent(in) :: ifun

THIS IS THE DIMENSION OF FUN IN THE DRIVER PROGRAM. IT MUST BE >= THE LARGEST INDEX I FOR WHICH FUN(I) IS USED UNLESS NO FUN(I) IS USED, IN WHICH CASE IT MUST BE >= 1.

real(kind=wp), intent(in) :: pttbl(Iptb,Indm)

(MATRIX ARRAY OF DIMENSION (IPTB,INDM)) ROW I OF THIS ARRAY NORMALLY CONTAINS A POINT USED IN THE ITH CONSTRAINT. THE ENTRIES IN ROW I NEED NOT BE ASSIGNED VALUES IF SUCH A POINT IS NOT USED IN THE ITH CONSTRAINT. (EXAMPLE: IF THE LEFT SIDE OF CONSTRAINT I IS A POLYNOMIAL IN ONE INDEPENDENT VARIABLE, THEN THE VALUE OF THE INDEPENDENT VARIABLE SHOULD BE IN PTTBL(I,1), AND THE COEFFICIENTS SHOULD BE IN PARAM.) PTTBL CAN ALSO BE USED TO PASS OTHER INFORMATION FROM THE DRIVER PROGRAM TO SUBROUTINE FNSET.

integer, intent(in) :: iptb

THIS IS THE FIRST DIMENSION OF PTTBL IN THE DRIVER PROGRAM. IT MUST BE >= THE LARGEST SUBSCRIPT I FOR WHICH A VALUE PTTBL(I,J) IS USED, AND MUST BE >= 1 IF NO SUCH VALUES ARE USED.

integer, intent(in) :: Indm

THIS IS THE SECOND DIMENSION OF PTTBL IN THE DRIVER PROGRAM. IT MUST BE >= THE LARGEST SUBSCRIPT J FOR WHICH A VALUE PTTBL(I,J) IS USED, AND MUST BE >= 1 IF NO SUCH VALUES ARE USED.

integer, intent(inout) :: iwork(Liwrk)

(VECTOR ARRAY OF DIMENSION LIWRK) THIS IS AN INTEGER WORK ARRAY. THE USER NEED NOT PLACE ANY VALUES IN IT, EXCEPT POSSIBLY CERTAIN OPTIONAL INFORMATION AS DESCRIBED BELOW.

integer, intent(in) :: liwrk

THIS IS THE DIMENSION OF IWORK IN THE DRIVER PROGRAM. IT MUST BE AT LEAST 7NUMGR + 7NPARM + 3. IF NOT, CONMAX WILL RETURN WITH THIS MINIMUM VALUE MULTIPLIED BY -1 AS A WARNING.

real(kind=wp), intent(inout) :: work(Lwrk)

(VECTOR ARRAY OF DIMENSION LWRK) THIS IS A FLOATING POINT WORK ARRAY. THE USER NEED NOT PLACE ANY VALUES IN IT, EXCEPT POSSIBLY CERTAIN OPTIONAL INFORMATION AS DESCRIBED BELOW.

integer, intent(in) :: lwrk

THIS IS THE DIMENSION OF WORK IN THE DRIVER PROGRAM. IT MUST BE AT LEAST 2NPARM2 + 4NUMGRNPARM + 11NUMGR + 27*NPARM + 13. IF NOT, CONMAX WILL RETURN WITH THIS MINIMUM VALUE MULTIPLIED BY -1 AS A WARNING.

integer, intent(out) :: iter

THIS IS THE NUMBER OF ITERATIONS PERFORMED BY CONMAX, INCLUDING THOSE USED IN ATTEMPTING TO GAIN FEASIBILITY, UNTIL EITHER IT CAN NO LONGER IMPROVE THE SITUATION OR THE ITERATION LIMIT IS REACHED. IF ITER=ITLIM IT IS POSSIBLE THAT THE PROGRAM COULD FURTHER REDUCE W IF RESTARTED (POSSIBLY WITH THE NEW PARAMETERS).

Read more…
real(kind=wp), intent(inout) :: param(Nparm)

(VECTOR ARRAY OF DIMENSION AT LEAST NPARM IN THE DRIVER PROGRAM) THE USER SHOULD PLACE AN INITIAL GUESS FOR THE PARAMETERS IN PARAM, AND ON OUTPUT PARAM WILL CONTAIN THE BEST PARAMETERS CONMAX HAS BEEN ABLE TO FIND. IF THE INITIAL PARAM IS NOT FEASIBLE THE PROGRAM WILL FIRST TRY TO FIND A FEASIBLE PARAM.

real(kind=wp), intent(out) :: error(Numgr+3)

(VECTOR ARRAY OF DIMENSION AT LEAST NUMGR + 3 IN THE DRIVER PROGRAM) FOR I=1,...,NUMGR, CONMAX WILL PLACE IN ERROR(I) THE ERROR IN CONSTRAINT I (DEFINED TO BE THE VALUE OF THE LEFT SIDE OF CONSTRAINT I, EXCEPT WITHOUT THE ABSOLUTE VALUE IN TYPE 2 CONSTRAINTS). FURTHER,

Read more…

private subroutine derst(me, Ioptn, Nparm, Numgr, Pttbl, Iptb, Indm, Param, Ipt, Param1, v, Kcntyp, Confun)

This subroutine uses fnset to compute confun(i,1) and the partial derivatives of the function whose value is in confun(i,1) for certain value(s) of i. note that we do not want the icntyp computed by fnset to override the icntyp (or jcntyp) carried into this subroutine in icntyp, so we use kcntyp when we call fnset. (the icntyp computed by fnset was stored earlier through a call to ercmp1 from conmax.)

Arguments

Type IntentOptional Attributes Name
class(conmax_solver), intent(inout) :: me
integer, intent(in) :: Ioptn
integer :: Nparm
integer :: Numgr
real(kind=wp), dimension(Iptb, Indm) :: Pttbl
integer :: Iptb
integer :: Indm
real(kind=wp), dimension(Nparm) :: Param
integer :: Ipt
real(kind=wp), dimension(Nparm) :: Param1
real(kind=wp), dimension(Numgr + 2*Nparm + 1, Nparm + 2) :: v
integer, dimension(Numgr) :: Kcntyp
real(kind=wp), dimension(Numgr, Nparm + 1) :: Confun

private subroutine slpcon(me, Ioptn, Nparm, Numgr, Fun, Ifun, Pttbl, Iptb, Indm, Tolcon, Rchin, Irk, Itypm1, Itypm2, Icntyp, Rchdwn, Numlim, Itersl, Prjslp, Funtbl, Iyrct, x, Mact1, Iact1, Jcntyp, Iphse, Enchg, Iwork, Liwrk, Work, Lwrk, Parser, Isucc, Param, Error)

Arguments

Type IntentOptional Attributes Name
class(conmax_solver), intent(inout) :: me
integer :: Ioptn
integer, intent(in) :: Nparm
integer, intent(in) :: Numgr
real(kind=wp) :: Fun(Ifun)
integer, intent(in) :: Ifun
real(kind=wp) :: Pttbl(Iptb,Indm)
integer, intent(in) :: Iptb
integer, intent(in) :: Indm
real(kind=wp) :: Tolcon
real(kind=wp) :: Rchin
integer :: Irk
integer :: Itypm1
integer :: Itypm2
integer :: Icntyp(Numgr)
real(kind=wp) :: Rchdwn
integer :: Numlim
integer :: Itersl
real(kind=wp) :: Prjslp
real(kind=wp) :: Funtbl(Numgr,Nparm+1)
integer :: Iyrct(Numgr+2*Nparm)
real(kind=wp) :: x(Nparm+1)
integer :: Mact1
integer :: Iact1(Numgr)
integer :: Jcntyp(Numgr)
integer :: Iphse
real(kind=wp) :: Enchg
integer :: Iwork(Liwrk)
integer, intent(in) :: Liwrk
real(kind=wp) :: Work(Lwrk)
integer, intent(in) :: Lwrk
real(kind=wp) :: Parser(Nparm)
integer :: Isucc
real(kind=wp) :: Param(Nparm)
real(kind=wp) :: Error(Numgr+3)

private subroutine bndset(Nparm, x, Itersl, Numin, Prjslp, Cofbnd, Xkeep, Bndkp)

This subroutine sets the bounds on the coefficient changes in slnpro.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: Nparm
real(kind=wp), intent(in), dimension(Nparm + 1) :: x
integer, intent(in) :: Itersl
integer, intent(in) :: Numin
real(kind=wp), intent(in) :: Prjslp
real(kind=wp), intent(inout), dimension(Nparm) :: Cofbnd
real(kind=wp), intent(inout), dimension(Nparm + 1) :: Xkeep
real(kind=wp), intent(inout), dimension(Nparm) :: Bndkp

private subroutine setu1(me, Ioptn, Numgr, Nparm, Numin, Rchin, Pttbl, Iptb, Indm, Fun, Ifun, Funtbl, Cofbnd, Param, Icntyp, Rchdwn, Error, Mact1, Iact1, Bndlgt, Iyrct, Iphse, Iwork, Liwrk, Work, Lwrk, Confun, Iact, v, m)

This subroutine sets up v for slnpro to solve a modified linearized (about the old parameters in param) version of our problem.

Arguments

Type IntentOptional Attributes Name
class(conmax_solver), intent(inout) :: me
integer :: Ioptn
integer :: Numgr
integer :: Nparm
integer :: Numin
real(kind=wp) :: Rchin
real(kind=wp), dimension(iptb, indm) :: Pttbl
integer :: Iptb
integer :: Indm
real(kind=wp), dimension(ifun) :: Fun
integer :: Ifun
real(kind=wp), dimension(numgr, nparm + 1) :: Funtbl
real(kind=wp), dimension(nparm) :: Cofbnd
real(kind=wp), dimension(nparm) :: Param
integer, dimension(numgr) :: Icntyp
real(kind=wp) :: Rchdwn
real(kind=wp), dimension(numgr + 3) :: Error
integer :: Mact1
integer, dimension(numgr) :: Iact1
real(kind=wp) :: Bndlgt
integer, dimension(numgr + 2*nparm) :: Iyrct
integer :: Iphse
integer, dimension(liwrk) :: Iwork
integer :: Liwrk
real(kind=wp), dimension(lwrk) :: Work
integer :: Lwrk
real(kind=wp), dimension(numgr, nparm + 1) :: Confun
integer, dimension(numgr) :: Iact
real(kind=wp), dimension(numgr + 2*nparm + 1, nparm + 2) :: v
integer :: m

public subroutine slnpro(v, m, n, Iyrct, y, Ixrct, Iycct, Nparm, Numgr, x, Indic)

This subroutine solves the linear programming problem

Read more…

Arguments

Type IntentOptional Attributes Name
real(kind=wp) :: v(Numgr+2*Nparm+1,Nparm+2)
integer :: m
integer :: n
integer :: Iyrct(Numgr+2*Nparm)
real(kind=wp) :: y(Numgr+2*Nparm)
integer :: Ixrct(Numgr+2*Nparm)
integer :: Iycct(Nparm+1)
integer, intent(in) :: Nparm
integer, intent(in) :: Numgr
real(kind=wp) :: x(Nparm+1)
integer :: Indic

private subroutine sjelim(l, Ll, k, Ir, Is, Nparm, Numgr, v)

This subroutine performs a modified jordan elimination on the l-ll+1 by k matrix consisting of rows ll through l of v and columns 1 through k of v. The resolvent is v(ir,is).

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: l
integer, intent(in) :: Ll
integer, intent(in) :: k
integer, intent(in) :: Ir
integer, intent(in) :: Is
integer, intent(in) :: Nparm
integer, intent(in) :: Numgr
real(kind=wp), intent(inout) :: v(Numgr+2*Nparm+1,Nparm+2)

private subroutine searsl(me, Ioptn, Numgr, Nparm, Prjlim, Tol1, x, Fun, Ifun, Pttbl, Iptb, Indm, Param, Error, Rchdwn, Mact, Iact, Iphse, Unit, Tolcon, Rchin, Itypm1, Itypm2, Iwork, Liwrk, Work, Lwrk, Err1, Parprj, Projct, Emin, Emin1, Parser, Nsrch)

This subroutine uses a modified quadratic fitting process to search for the minimum of a function f. it requres an initial guess in projct, a tolerance tol1 on the search interval length, an upper bound prjlim on the minimizing point (which should be set very large if no upper bound is desired), and a way to compute f(x) for a given x. the subroutine will print a warning and return if it would need to compute f more than initlm times in the initialization or more than nadd additional times in the main part of the program. when the subroutine returns, it will have put the minimum location in projct, the minimum f value in emin, the f value for the initial projct in emin1, and the number of times it computed f in nsrch.

Arguments

Type IntentOptional Attributes Name
class(conmax_solver), intent(inout) :: me
integer :: Ioptn
integer :: Numgr
integer :: Nparm
real(kind=wp) :: Prjlim
real(kind=wp) :: Tol1
real(kind=wp), dimension(nparm + 1) :: x
real(kind=wp), dimension(ifun) :: Fun
integer :: Ifun
real(kind=wp), dimension(iptb, indm) :: Pttbl
integer :: Iptb
integer :: Indm
real(kind=wp), dimension(nparm) :: Param
real(kind=wp), dimension(numgr + 3) :: Error
real(kind=wp) :: Rchdwn
integer :: Mact
integer, dimension(numgr) :: Iact
integer :: Iphse
real(kind=wp) :: Unit
real(kind=wp) :: Tolcon
real(kind=wp) :: Rchin
integer :: Itypm1
integer :: Itypm2
integer, dimension(liwrk) :: Iwork
integer :: Liwrk
real(kind=wp), dimension(lwrk) :: Work
integer :: Lwrk
real(kind=wp), dimension(numgr + 3) :: Err1
real(kind=wp), dimension(nparm) :: Parprj
real(kind=wp) :: Projct
real(kind=wp) :: Emin
real(kind=wp) :: Emin1
real(kind=wp), dimension(nparm) :: Parser
integer :: Nsrch

private subroutine ercmp1(me, Ioptn, Nparm, Numgr, Fun, Ifun, Pttbl, Iptb, Indm, Param, Icnuse, Iphse, Iwork, Liwrk, Confun, Icntyp, Ipmax, Ismax, Error)

Arguments

Type IntentOptional Attributes Name
class(conmax_solver), intent(inout) :: me
integer, intent(in) :: Ioptn
integer, intent(in) :: Nparm
integer, intent(in) :: Numgr
real(kind=wp), intent(in) :: Fun(Ifun)
integer, intent(in) :: Ifun
real(kind=wp), intent(in) :: Pttbl(Iptb,Indm)
integer, intent(in) :: Iptb
integer, intent(in) :: Indm
real(kind=wp), intent(in) :: Param(Nparm)
integer, intent(in) :: Icnuse
integer, intent(in) :: Iphse
integer :: Iwork(Liwrk)
integer, intent(in) :: Liwrk
real(kind=wp) :: Confun(Numgr,Nparm+1)
integer :: Icntyp(Numgr)
integer, intent(out) :: Ipmax
integer, intent(out) :: Ismax
real(kind=wp) :: Error(Numgr+3)

private subroutine rkcon(me, Ioptn, Nparm, Numgr, Fun, Ifun, Pttbl, Iptb, Indm, Tolcon, Rchin, Iter, Irk, Ityp2, Ityp1, Itypm1, Itypm2, Icntyp, Projct, Rchdwn, Nstep, Iphse, Enchg, Enc1, Pmat, Funtbl, Iwork, Liwrk, Work, Lwrk, Iact, Actdif, Parprj, Parser, Xrk, Err1, Confun, Isucc, Param, Error)

Arguments

Type IntentOptional Attributes Name
class(conmax_solver), intent(inout) :: me
integer :: Ioptn
integer, intent(in) :: Nparm
integer, intent(in) :: Numgr
real(kind=wp) :: Fun(Ifun)
integer, intent(in) :: Ifun
real(kind=wp) :: Pttbl(Iptb,Indm)
integer, intent(in) :: Iptb
integer :: Indm
real(kind=wp) :: Tolcon
real(kind=wp) :: Rchin
integer :: Iter
integer :: Irk
integer :: Ityp2
integer :: Ityp1
integer :: Itypm1
integer :: Itypm2
integer :: Icntyp(Numgr)
real(kind=wp) :: Projct
real(kind=wp) :: Rchdwn
integer :: Nstep
integer :: Iphse
real(kind=wp) :: Enchg
real(kind=wp) :: Enc1
real(kind=wp) :: Pmat(Nparm+1,Numgr)
real(kind=wp) :: Funtbl(Numgr,Nparm+1)
integer :: Iwork(Liwrk)
integer, intent(in) :: Liwrk
real(kind=wp) :: Work(Lwrk)
integer, intent(in) :: Lwrk
integer :: Iact(Numgr)
real(kind=wp) :: Actdif(Numgr)
real(kind=wp) :: Parprj(Nparm)
real(kind=wp) :: Parser(Nparm)
real(kind=wp) :: Xrk(Nparm+1)
real(kind=wp) :: Err1(Numgr+3)
real(kind=wp) :: Confun(Numgr,Nparm+1)
integer :: Isucc
real(kind=wp) :: Param(Nparm)
real(kind=wp) :: Error(Numgr+3)

private subroutine rksact(Ioptn, Numgr, Icntyp, Rchdwn, Rchin, Conup, Projct, Error, Mactrk, Actdif, Iact)

This subroutine puts the (signed) indices of the mactrk active constraints in iact. it also sets the right side vector actdif for the wolfe subproblem.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: Ioptn
integer, intent(in) :: Numgr
integer, intent(in) :: Icntyp(Numgr)
real(kind=wp), intent(in) :: Rchdwn
real(kind=wp), intent(in) :: Rchin
real(kind=wp), intent(in) :: Conup
real(kind=wp), intent(in) :: Projct
real(kind=wp), intent(in) :: Error(Numgr+3)
integer, intent(out) :: Mactrk
real(kind=wp) :: Actdif(Numgr)
integer :: Iact(Numgr)

private subroutine pmtst(me, Ioptn, Numgr, Nparm, Param, Icntyp, Mactrk, Iact, Pttbl, Iptb, Indm, Actdif, Iphse, Iwork, Liwrk, Work, Lwrk, Confun, Pmat)

This subroutine sets up the (nparm+1) by mactrk matrix pmat. for 1 <= j <= mactrk, the top nparm elements of column j of pmat will contain the negative of the gradient of active constraint j (if constraint j is of type 2, i.e. of the form abs(f(x) - f(parwrk,x)) <= w, the left side will be treated as f(x) - f(parwrk,x) if this quantity is nonnegative and will be treated as f(parwrk,x) - f(x) otherwise). the (nparm+1)st row of pmat will contain actdif, the right side of the inequalities gradient.vector >= actdif.

Arguments

Type IntentOptional Attributes Name
class(conmax_solver), intent(inout) :: me
integer, intent(in) :: Ioptn
integer, intent(in) :: Numgr
integer, intent(in) :: Nparm
real(kind=wp) :: Param(Nparm)
integer :: Icntyp(Numgr)
integer, intent(in) :: Mactrk
integer :: Iact(Numgr)
real(kind=wp) :: Pttbl(Iptb,Indm)
integer, intent(in) :: Iptb
integer, intent(in) :: Indm
real(kind=wp), intent(in) :: Actdif(Numgr)
integer, intent(in) :: Iphse
integer :: Iwork(Liwrk)
integer, intent(in) :: Liwrk
real(kind=wp) :: Work(Lwrk)
integer, intent(in) :: Lwrk
real(kind=wp) :: Confun(Numgr,Nparm+1)
real(kind=wp) :: Pmat(Nparm+1,Numgr)

private subroutine rkpar(me, Ioptn, Numgr, Nparm, Icntyp, Mactrk, Iact, Actdif, Projct, Param, Fun, Ifun, Pttbl, Iptb, Indm, Vder, Pmat, Ncor, s, Itypm1, Itypm2, Unit, Tolcon, Rchin, Nstep, Error, Iphse, Iwork, Liwrk, Work, Lwrk, Confun, Vdern, Vders, Wvec, Parprj, Ifrkpr)

This subroutine computes a parameter vector parprj using fourth order runge kutta with h=-projct. h is negative since we want to approximate the parameters resulting from decreasing w by abs(h). if we do nstep steps then h=-projct/nstep.

Arguments

Type IntentOptional Attributes Name
class(conmax_solver), intent(inout) :: me
integer :: Ioptn
integer, intent(in) :: Numgr
integer, intent(in) :: Nparm
integer :: Icntyp(Numgr)
integer :: Mactrk
integer :: Iact(Numgr)
real(kind=wp) :: Actdif(Numgr)
real(kind=wp) :: Projct
real(kind=wp) :: Param(Nparm)
real(kind=wp) :: Fun(Ifun)
integer, intent(in) :: Ifun
real(kind=wp) :: Pttbl(Iptb,Indm)
integer, intent(in) :: Iptb
integer, intent(in) :: Indm
real(kind=wp) :: Vder(Nparm)
real(kind=wp) :: Pmat(Nparm+1,Numgr)
integer :: Ncor
real(kind=wp) :: s
integer :: Itypm1
integer :: Itypm2
real(kind=wp) :: Unit
real(kind=wp) :: Tolcon
real(kind=wp) :: Rchin
integer :: Nstep
real(kind=wp) :: Error(Numgr+3)
integer :: Iphse
integer :: Iwork(Liwrk)
integer, intent(in) :: Liwrk
real(kind=wp) :: Work(Lwrk)
integer, intent(in) :: Lwrk
real(kind=wp) :: Confun(Numgr,Nparm+1)
real(kind=wp) :: Vdern(Nparm)
real(kind=wp) :: Vders(Nparm)
real(kind=wp) :: Wvec(Nparm)
real(kind=wp) :: Parprj(Nparm)
integer :: Ifrkpr

private subroutine corrct(me, Ioptn, Nparm, Numgr, Fun, Ifun, Pttbl, Iptb, Indm, Icntyp, Unit, Tolcon, Rchin, Error, Mact, Iact, Projct, Iphse, Iwork, Liwrk, Work, Lwrk, Parwrk, Err1, Dvec, Pmat, Confun, Zwork, Jcntyp, Parprj, Icorct)

This subroutine determines whether parprj violates any type -2 or type -1 (i.e. standard) constraints by more than tolcon, and if so it attempts to correct back to the feasible region. if it is successful it sets icorct=0 and replaces parprj by the corrected vector. if it is not successful it sets icorct=1 and leaves parprj unchanged. if no correction was needed it sets icorct=-1 and leaves parprj unchanged.

Arguments

Type IntentOptional Attributes Name
class(conmax_solver), intent(inout) :: me
integer :: Ioptn
integer, intent(in) :: Nparm
integer, intent(in) :: Numgr
real(kind=wp) :: Fun(Ifun)
integer, intent(in) :: Ifun
real(kind=wp) :: Pttbl(Iptb,Indm)
integer, intent(in) :: Iptb
integer :: Indm
integer :: Icntyp(Numgr)
real(kind=wp) :: Unit
real(kind=wp) :: Tolcon
real(kind=wp) :: Rchin
real(kind=wp) :: Error(Numgr+3)
integer :: Mact
integer :: Iact(Numgr)
real(kind=wp) :: Projct
integer :: Iphse
integer :: Iwork(Liwrk)
integer, intent(in) :: Liwrk
real(kind=wp) :: Work(Lwrk)
integer, intent(in) :: Lwrk
real(kind=wp) :: Parwrk(Nparm)
real(kind=wp) :: Err1(Numgr+3)
real(kind=wp) :: Dvec(Nparm)
real(kind=wp) :: Pmat(Nparm+1,Numgr)
real(kind=wp) :: Confun(Numgr,Nparm+1)
real(kind=wp) :: Zwork(Nparm)
integer :: Jcntyp(Numgr)
real(kind=wp) :: Parprj(Nparm)
integer :: Icorct

private subroutine searcr(me, Ioptn, Nparm, Numgr, Dvec, Fun, Ifun, Pttbl, Iptb, Indm, Zwork, Tolcon, Iphse, Iwork, Liwrk, Work, Lwrk, Parwrk, Err1, p1, f1, Procor, Emin, Isrcr)

This subroutine uses a modified quadratic fitting process to search for a projection factor procor for which the maximum of the left sides of the type -2 and -1 constraints evaluated at zwork + procor*dvec is <= tolcon. note that when corrct calls this subroutine it will have lumped the type -1 constraints in with the type -2 constraints using jcntyp, which is carried through this subroutine into subroutine ercmp1 in iwork. if searcr is able to force this maximum <= tolcon it will return with isrcr=0, with the minimum value found for the maximum in emin, with the corresponding projection factor in procor, with the number of times the maximum was computed in nsrch, and with the closest point found to the left with the maximum > tolcon in (p1,f1). the subroutine will begin by computing the maxima for procor = 1.0, 0.5, and 2.0, and if none of these maxima is <= tolcon and it is not the case that the maximum at 1.0 is <= the other two maxima the subroutine will return with the warning isrcr=1. the subroutine will also return with isrcr=1 if it would need to compute f more than limscr times, or the search interval length drops below tol1, or the quadratic fit becomes too flat. even in the event of a return with isrcr=1, emin, procor, and nsrch will be as above.

Arguments

Type IntentOptional Attributes Name
class(conmax_solver), intent(inout) :: me
integer :: Ioptn
integer, intent(in) :: Nparm
integer :: Numgr
real(kind=wp) :: Dvec(Nparm)
real(kind=wp) :: Fun(Ifun)
integer, intent(in) :: Ifun
real(kind=wp) :: Pttbl(Iptb,Indm)
integer, intent(in) :: Iptb
integer, intent(in) :: Indm
real(kind=wp) :: Zwork(Nparm)
real(kind=wp) :: Tolcon
integer :: Iphse
integer :: Iwork(Liwrk)
integer, intent(in) :: Liwrk
real(kind=wp) :: Work(Lwrk)
integer, intent(in) :: Lwrk
real(kind=wp) :: Parwrk(Nparm)
real(kind=wp) :: Err1(Numgr+3)
real(kind=wp) :: p1
real(kind=wp) :: f1
real(kind=wp) :: Procor
real(kind=wp) :: Emin
integer :: Isrcr

private subroutine muller(me, Ioptn, Nparm, Numgr, Dvec, Fun, Ifun, Pttbl, Iptb, Indm, Zwork, Tolcon, Iphse, Iwork, Liwrk, Work, Lwrk, Parwrk, Err1, p1, f1, Procor, Emin)

In this subroutine we are given a base vector zwork, a direction vector dvec, a scalar procor with emin = f(procor) = (the maximum type -2 and -1 error with parameters zwork + procor*dvec) < -tolcon, and a scalar p1 with p1 < procor and f1 = f(p1) > tolcon. we do a revised mullers method approach (with a solution contained in a shrinking interval) to attempt to adjust procor so that -tolcon <= f(procor) <= tolcon, but if we are not successful we return with the leftmost procor found satisfying emin = f(procor) < -tolcon on the theory that overcorrection is better than no correction. note that when corrct calls this subroutine it will have lumped the type -1 constraints in with the type -2 constraints using jcntyp, which is carried through this subroutine into subroutine ercmp1 in iwork.

Arguments

Type IntentOptional Attributes Name
class(conmax_solver), intent(inout) :: me
integer :: Ioptn
integer, intent(in) :: Nparm
integer, intent(in) :: Numgr
real(kind=wp) :: Dvec(Nparm)
real(kind=wp) :: Fun(Ifun)
integer, intent(in) :: Ifun
real(kind=wp) :: Pttbl(Iptb,Indm)
integer, intent(in) :: Iptb
integer, intent(in) :: Indm
real(kind=wp) :: Zwork(Nparm)
real(kind=wp) :: Tolcon
integer :: Iphse
integer :: Iwork(Liwrk)
integer :: Liwrk
real(kind=wp) :: Work(Lwrk)
integer, intent(in) :: Lwrk
real(kind=wp) :: Parwrk(Nparm)
real(kind=wp) :: Err1(Numgr+3)
real(kind=wp) :: p1
real(kind=wp) :: f1
real(kind=wp) :: Procor
real(kind=wp) :: Emin

private subroutine rchmod(Numgr, Error, Err1, Icntyp, Mact, Iact, Ipmax, Ismax, Unit, Irch, Rchdwn, Rchin)

This subroutine increases rchdwn or rchin if it appears some constraints which should have been declared active were not so declared.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: Numgr
real(kind=wp), intent(in) :: Error(Numgr+3)
real(kind=wp), intent(in) :: Err1(Numgr+3)
integer, intent(in) :: Icntyp(Numgr)
integer, intent(in) :: Mact
integer, intent(in) :: Iact(Numgr)
integer, intent(in) :: Ipmax
integer, intent(in) :: Ismax
real(kind=wp), intent(in) :: Unit
integer, intent(in) :: Irch
real(kind=wp), intent(inout) :: Rchdwn
real(kind=wp), intent(inout) :: Rchin

public subroutine wolfe(Ndm, m, Pmat, Istrt, s, Ncor, Icor, Iwork, Liwrk, Work, Lwrk, r, Coef, Ptnr, Pmat1, Nparm, Numgr, Wcoef, Wpt, Wdist, Nmaj, Nmin, Jflag)

given m inequalities of the form a(k).x + b(k) <= 0.0 for k=1, ...,m, where a(k) and x are ndm dimensional vectors and b(k) are numbers, this subroutine returns the nearest point to the origin in the polytope defined by these inequalities (unless jflag > 0, which indicates failure). the user should put the mdm+1 dimensional vectors (a(k),b(k)) in the columns of pmat. the solution point will be returned in wpt, and will also be a linear combination of the a(k) vectors with (nonpositive) coefficients in the m dimensional vector wcoef. wcoef may not be accurate if refwl was used to refine wpt, which rarely happens. the number of vectors in the final corral will be returned in ncor with their indices in icor, and all entries of wcoef not corresponding to indices in icor will be zero. the distance will be returned in wdist, and the numbers of major and minor cycles in the cone subproblem will be returned in nmaj and nmin respectively. if the user sets istrt=0 the program will start from scratch, but the user can set istrt=1 (hot start) and specify ncor, icor, wcoef, and the factor s. (see later comments; set s=1.0 if no better value is available. set wcoef(j)=0.0 if icor(i) /= j for i=1,...,ncor.) (if inaccurate wcoef or s is used in a hot start attempt little will be lost, since ncor and icor are more important for a successful hot start than wcoef and s.) we must always have ncor <= ndm+1 in theory since the ncor ndm+1 dimensional vectors in a corral should be linearly independent, and in practice we will always require ncor <= ndm+1. if the user sets istrt=1 but the program fails, it will automatically try from scratch before giving up.

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: Ndm

The number of variables. it must be less than or equal to nparm.

integer, intent(in) :: m

The number of inequalities defining the polytope. it must be less than or equal to numgr.

real(kind=wp), intent(in) :: Pmat(Nparm+1,Numgr)

This is an array whose kth column contains the vector (a(k),b(k)) for k=1,...,m, where the m inequalities a(k).x + b(k) <= 0.0 define the polytope whose nearest point to the origin we seek. the first dimension of pmat in the driver program must be exactly nparm+1, while the second dimension of pmat in the driver program must be at least numgr. if we actually want the nearest point in the polytope to some point y other than the origin, we translate y to the origin before calling wolfe, that is, call wolfe to find the nearest point z to the origin in the polytope defined by a(k).z + (b(k) + a(k).y) <= 0.0, then compute x = y + z.

integer, intent(in) :: Istrt

Set this equal to zero unless a hot start is desired-- see next paragraph of comments for more details. if istrt is set equal to 1, then s, wcoef, ncor, and icor must also be assigned values initially.

real(kind=wp), intent(out) :: s

You may ignore this scale factor unless you want to use the hot start option.

integer, intent(out) :: Ncor

This is the number of vectors (i.e. colunns of pmat) in the final corral.

integer, intent(out) :: Icor(Nparm+1)

This array contains the ncor indices of the vectors in the final corral. its dimension in the driver program must be at least nparm+1.

integer :: Iwork(Liwrk)

Work array

integer, intent(in) :: Liwrk

This is the dimension of iwork. It must be at least 7nparm + 7numgr + 3.

real(kind=wp) :: Work(Lwrk)

Work array

integer, intent(in) :: Lwrk

This is the dimension of work. it must be at least 2nparm2 + 4numgrnparm + 11numgr + 27*nparm + 13. note that some storage could be saved by rewriting function subprogram iloc to take out all but the arrays needed (namely 1, 3, 4, 9, 28, 32, 34, 39 for work, 18, 23 for iwork) and scrunching work and iwork in iloc so the remaining arrays follow one after another.

real(kind=wp) :: r(Nparm+1)

Work array

real(kind=wp) :: Coef(Numgr)

Work array

real(kind=wp) :: Ptnr(Nparm+1)

Work array

real(kind=wp) :: Pmat1(Nparm+1,Numgr)

Work array

integer, intent(in) :: Nparm

This is basically a dimension parameter here. it must be greater than or equal to ndm.

integer, intent(in) :: Numgr

This is basically a dimension parameter here. it must be greater than or equal to m.

real(kind=wp), intent(out) :: Wcoef(Numgr)

This will give the coefficients of the vectors a(k) needed to form a linear combination equal to the solution in wpt. its dimension in the driver program must be at least numgr. wcoef may not be accurate if it was necessary to call refwl to refine wpt, which rarely happens.

real(kind=wp), intent(out) :: Wpt(Nparm)

This will give the coordinates of the point we are seeking, namely the nearest point in the polytope to the origin. its dimension in the driver program must be at least nparm.

real(kind=wp), intent(out) :: Wdist

This will be the (minimized) euclidean distance of wpt from the origin.

integer, intent(out) :: Nmaj

This will be the number of major cycles used in wolfe.

integer, intent(out) :: Nmin

This will be the number of major cycles used in wolfe.

integer, intent(out) :: Jflag

This is a flag variable which is 0 in case of a normal solution and is positive otherwise (in which case the returned solution may be no good).

private subroutine conenr(n, m, Pmat1, r, Istrt1, Ncor, Icor, Tol, Iwork, Liwrk, Work, Lwrk, Vec, Ptnrr, Picor, Nparm, Numgr, Coef, Ptnr, Dist, Nmaj, Nmin, Jflag)

Given m n-dimensional vectors p(j) as the first m columns of the matrix pmat1 and an n-vector r, this subroutine returns in ptnr the nearest point to r in the cone of points summation( coef(j)*p(j)), where coef(j) >= 0.0 for j=1,...,m (unless jflag > 0, which indicates failure). the number of vectors p(j) in the final corral is returned in ncor with their indices in icor, the distance is returned in dist, the number of major cycles (i.e. adding a vector) is returned in nmaj, and the number of minor cycles (i.e. removing a vector) is returned in nmin. if the user sets istrt1=0 the subroutne starts from scratch, but the user can set istrt1=1 and initially specify ncor, icor, and coef (noting that ncor must be <= n, and if j does not occur in icor, then coef(j) should be set to 0.0.)

Arguments

Type IntentOptional Attributes Name
integer :: n
integer :: m
real(kind=wp) :: Pmat1(Nparm+1,Numgr)
real(kind=wp) :: r(Nparm+1)
integer :: Istrt1
integer :: Ncor
integer :: Icor(Nparm+1)
real(kind=wp) :: Tol
integer :: Iwork(Liwrk)
integer, intent(in) :: Liwrk
real(kind=wp) :: Work(Lwrk)
integer, intent(in) :: Lwrk
real(kind=wp) :: Vec(Nparm+1)
real(kind=wp) :: Ptnrr(Nparm+1)
real(kind=wp) :: Picor(Nparm+1,Nparm+1)
integer, intent(in) :: Nparm
integer, intent(in) :: Numgr
real(kind=wp) :: Coef(Numgr)
real(kind=wp) :: Ptnr(Nparm+1)
real(kind=wp) :: Dist
integer :: Nmaj
integer :: Nmin
integer :: Jflag

private subroutine house(n, Ncor, Picor, r, Kpivot, Nparm, Aa, Beta, d, Save, b, Vec, Ihouse)

Given ncor n dimensional vectors as columns of the n by ncor matrix picor and an n dimensional vector r, this subroutine uses householder transformations to find the best least squares solution vec to the linear system of equations picor*vec = r, where vec is an ncor dimensional vector. if the rank of picor is (computationally) 0, the subroutine will return with the failure warning ihouse=1, otherwise it will return with ihouse=0. if the rank is > 0 but < ncor, then (ncor - rank) of the components of vec will be set to 0.0. the arays picor and r will not be changed by this subroutine. the subroutine will attempt up to numref iterative refinements of the solution, where the user can set numref as any nonnegative integer, but to get the most out of the iterative refinement process, the computation of the residual summ near the end of this subroutine should be done in higher precision than the other computations in the subroutine.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: n
integer, intent(in) :: Ncor
real(kind=wp) :: Picor(Nparm+1,Nparm+1)
real(kind=wp) :: r(Nparm+1)
integer :: Kpivot(Nparm+1)
integer, intent(in) :: Nparm
real(kind=wp) :: Aa(Nparm+1,Nparm+1)
real(kind=wp) :: Beta(Nparm+1)
real(kind=wp) :: d(Nparm+1)
real(kind=wp) :: Save(Nparm+1)
real(kind=wp) :: b(Nparm+1)
real(kind=wp) :: Vec(Nparm+1)
integer, intent(out) :: Ihouse

private subroutine refwl(Ndm, Ncor, Icor, Pmat, Pmat1, Nparm, Numgr, Ixrct, Save, Wpt)

This subroutine attempts to refine the ndm dimensional vector wpt produced by wolfe by directly solving the system summation(pmat(i,j)*wpt(i), i=1,...,ndm) = -pmat(ndm+1,j) for j = icor(l), l=1,...,ncor. nresl resolvents are chosen by total pivoting. if nresl < ndm then the remaining ndm-nresl elements of wpt are kept form the old wpt. itrlm steps of iterative refinement are attempted at the end.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: Ndm
integer, intent(in) :: Ncor
integer :: Icor(Nparm+1)
real(kind=wp) :: Pmat(Nparm+1,Numgr)
real(kind=wp) :: Pmat1(Nparm+1,Numgr)
integer, intent(in) :: Nparm
integer, intent(in) :: Numgr
integer :: Ixrct(2*Nparm)
real(kind=wp) :: Save(Nparm)
real(kind=wp) :: Wpt(Nparm+1)