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)
implicit none
class(conmax_solver), intent(inout) :: me
integer, intent(in) :: Ifun
integer, intent(in) :: Iptb
integer, intent(in) :: Liwrk
integer, intent(in) :: Lwrk
integer, intent(in) :: Nparm
integer, intent(in) :: Numgr
integer :: Indm
integer :: Ioptn
integer :: Iphse
integer :: Irk
integer :: Isucc
integer :: Iter
integer :: Ityp1
integer :: Ityp2
integer :: Itypm1
integer :: Itypm2
integer :: Nstep
real(wp) :: Enc1
real(wp) :: Enchg
real(wp) :: Projct
real(wp) :: Rchdwn
real(wp) :: Rchin
real(wp) :: Tolcon
integer :: Iact(Numgr)
integer :: Icntyp(Numgr)
integer :: Iwork(Liwrk)
real(wp) :: Actdif(Numgr)
real(wp) :: Confun(Numgr, Nparm + 1)
real(wp) :: Err1(Numgr + 3)
real(wp) :: Error(Numgr + 3)
real(wp) :: Fun(Ifun)
real(wp) :: Funtbl(Numgr, Nparm + 1)
real(wp) :: Param(Nparm)
real(wp) :: Parprj(Nparm)
real(wp) :: Parser(Nparm)
real(wp) :: Pmat(Nparm + 1, Numgr)
real(wp) :: Pttbl(Iptb, Indm)
real(wp) :: Work(Lwrk)
real(wp) :: Xrk(Nparm + 1)
real(wp) :: conup, emin, emin1, enorm, pe, &
prjbig, prjlim, prosea, qt, &
quots, s, ss, steplm, unit, wdist
integer :: i, icorct, ifrkpr, ilc06, &
ilc10, ilc15, ilc21, ilc22, ilc24, ilc27, ilc30, &
ilc31, ilc33, ilc35, ilc36, ilc37, ilc38, &
ilc40, ilc43, ilc48, ioptth, &
ipmax, ipt, ismax, iwarn, &
j, jflag, l, limfl, mactrk, &
ncor, nfail, nmaj, nmin, npar1, nsrch
real(wp), parameter :: qthi = (one + two)/four
real(wp), parameter :: qtlo = one/four
real(wp), parameter :: tol1 = ten*ten*spcmn
real(wp), parameter :: tol2 = ten*spcmn
real(wp), parameter :: prden = sqrt(sqrt(spcmn))
ioptth = (Ioptn - (Ioptn/100000)*100000)/10000
steplm = Tolcon/ten
ilc06 = iloc(6, Nparm, Numgr)
ilc10 = iloc(10, Nparm, Numgr)
ilc15 = iloc(15, Nparm, Numgr)
ilc21 = iloc(21, Nparm, Numgr)
ilc22 = iloc(22, Nparm, Numgr)
ilc24 = iloc(24, Nparm, Numgr)
ilc27 = iloc(27, Nparm, Numgr)
ilc30 = iloc(30, Nparm, Numgr)
ilc31 = iloc(31, Nparm, Numgr)
ilc33 = iloc(33, Nparm, Numgr)
ilc35 = iloc(35, Nparm, Numgr)
ilc36 = iloc(36, Nparm, Numgr)
ilc37 = iloc(37, Nparm, Numgr)
ilc38 = iloc(38, Nparm, Numgr)
ilc40 = iloc(40, Nparm, Numgr)
ilc43 = iloc(43, Nparm, Numgr)
ilc48 = iloc(48, Nparm, Numgr)
Isucc = 0
iwarn = 0
nfail = 0
conup = one
! LIMFL IS A SAFETY VALVE TO CATCH BLUNDERS; WE SET IT HIGH ENOUGH
! THAT IT WILL NOMALLY NOT COME INTO PLAY.
limfl = 20
enorm = Error(Numgr + 1)
npar1 = Nparm + 1
prjbig = one/spcmn
if (Ityp2 > 0) prjbig = enorm
! THE NEXT GROUP OF STATEMENTS SETS AN INITIAL VALUE FOR PROJCT.
if (Iter > 0) then
if (Irk <= 1) then
! HERE ITER > 0 AND IRK=1, AND WE BUILD ON THE PREVIOUS SUCCESSFUL
! RK ITERATION, WHICH REDUCED THE ERROR NORM. COMPUTE THE RATIO QT,
! WHICH WOULD BE 1.0 IF RUNGE-KUTTA WERE EXACT AND NO CORRECTION STEP
! WERE NEEDED.
qt = -Enc1/Projct
if (qt >= qthi) then
! HERE QT >= QTHI, SO WE INCREASE PROJCT BY A FACTOR OF 2.0.
Projct = two*Projct
! IF QTLO < QT < QTHI WE LEAVE PROJCT THE SAME, WHILE IF QT <=
! QTLO WE DIVIDE PROJCT BY 4.0.
else if (qt <= qtlo) then
Projct = Projct/four
end if
goto 100
end if
end if
! HERE ITER=0, OR ELSE ITER > 0 AND IRK=2, AND WE INITIALIZE PROJCT.
if (Iphse + 1 < 0) then
! HERE ITER=0 OR IRK=2, AND IPHSE=-2, SO WE ARE ATTEMPTING TO GAIN TYPE -2
! FEASIBILITY, AND WE SET THE INITIAL PROJCT TO ENOR3,
! WHICH WILL BE > TOLCON. NOTE THAT ENOR3 IS NOW IN ERROR(NUMGR+1).
Projct = enorm
else if (Iphse + 1 /= 0) then
! HERE ITER=0 OR IRK=2, AND IPHSE=0, SO WE ARE IN THE MAIN ITERATIONS,
! AND WE FIRST TRY PROJCT=1.0.
Projct = one
! CHECK TO SEE WHETHER ABS(ENORM) IS VERY
! LARGE RELATIVE TO THE INITIAL PROJCT. IF ABS(ENORM) >
! PROJCT/PRDEN, WE REPLACE THE INITIAL PROJCT BY PRDEN*ABS(ENORM)
! SO THAT IF WE ARE SUCCESSFUL IN REDUCING ENORM TO ENORM - PROJCT,
! THIS QUANTITY WILL DIFFER FROM ENORM IN AT LEAST SOME SIGNIFICANT
! DIGITS AND THE PROGRAM WILL HAVE A CHANCE TO CONTINUE.
pe = prden*abs(enorm)
if (pe > Projct) Projct = pe
! IF ITYP2 > 0 WE REDUCE THE INITIAL PROJCT TO ENORM (IF NECESSARY),
! WHICH WILL BE THE GREATEST DECREASE IN ENORM WE CAN HOPE FOR SINCE
! THERE WILL BE TYPE 2 CONSTRAINTS.
if (Ityp2 > 0) then
if (enorm < Projct) Projct = enorm
end if
end if
! WE DO NOT WISH FOR PROJCT TO BE SET BELOW 100.0*SPCMN
if (Projct < ten*ten*spcmn) Projct = ten*ten*spcmn
! WE DO NOT WANT PROJCT TO BE BIGGER THAN PRJBIG OR SMALLER THAN
! STEPLM.
100 if (Projct > prjbig) Projct = prjbig
if (Projct < steplm) then
iwarn = 1
Projct = steplm
end if
! CALL RKSACT TO PUT THE (SIGNED) INDICES OF THE ACTIVE CONSTRAINTS IN
! IACT AND THEIR NUMBER IN MACTRK.
call rksact(Ioptn, Numgr, Icntyp, Rchdwn, Rchin, conup, Projct, Error, &
mactrk, Actdif, Iact)
! SET UNIT FOR USE IN RCHMOD. UNIT WILL BE THE VALUE OF PROJCT WHEN
! RKSACT WAS LAST CALLED.
unit = Projct
! CALL PMTST TO SET UP PMAT.
call me%pmtst(Ioptn, Numgr, Nparm, Param, Icntyp, mactrk, Iact, Pttbl, Iptb, &
Indm, Actdif, Iphse, Iwork, Liwrk, Work, Lwrk, Confun, Pmat)
! COPY PMAT TRANSPOSE INTO FUNTBL FOR POSSIBLE LATER USE.
do j = 1, npar1
do i = 1, mactrk
Funtbl(i, j) = Pmat(j, i)
end do
end do
! STATEMENTS ABOVE THIS POINT WILL NOT BE EXECUTED AGAIN IN THIS CALL
! TO RKCON.
! INCREMENT NFAIL, WHICH COUNTS THE NUMBER OF WOLFE CALLS IN THIS CALL TO
! RKCON.
200 nfail = nfail + 1
! CALL WOLFE WITH ISTRT=0 TO SOLVE THE LEAST DISTANCE QP PROBLEM FROM
! SCRATCH.
call wolfe(Nparm, mactrk, Pmat, 0, s, ncor, Iwork(ilc15), Iwork, Liwrk, &
Work, Lwrk, Work(ilc33), Work(ilc06), Work(ilc31), &
Work(ilc30), Nparm, Numgr, Work(ilc40), Work(ilc36), wdist, &
nmaj, nmin, jflag)
! IF WOLFE FAILS, WE MAY TRY AGAIN WITH A SMALLER PROJCT.
if (jflag <= 0) then
! END OF GROUP OF STATEMENTS TO REDUCE PROJCT (IF POSSIBLE) TO HANDLE
! A FAILURE OF SOME KIND.
! DO AN RK STEP.
call me%rkpar(Ioptn, Numgr, Nparm, Icntyp, mactrk, Iact, Actdif, Projct, &
Param, Fun, Ifun, Pttbl, Iptb, Indm, Work(ilc36), Pmat, &
ncor, s, Itypm1, Itypm2, unit, Tolcon, Rchin, Nstep, Error, &
Iphse, Iwork, Liwrk, Work, Lwrk, Confun, Work(ilc37), &
Work(ilc38), Work(ilc43), Parprj, ifrkpr)
if (ifrkpr <= 0) then
! HERE RKPAR SUCCEEDED. IF THERE ARE ANY STANDARD CONSTRAINTS WE CALL
! CORRCT TO MOVE BACK INTO THE FEASIBLE REGION IF NECESSARY.
if (Itypm1 + Itypm2 <= 0) goto 500
call me%corrct(Ioptn, Nparm, Numgr, Fun, Ifun, Pttbl, Iptb, Indm, &
Icntyp, unit, Tolcon, Rchin, Error, mactrk, Iact, &
Projct, Iphse, Iwork, Liwrk, Work, Lwrk, Work(ilc27), &
Err1, Work(ilc10), Pmat, Confun, Work(ilc48), &
Iwork(ilc21), Parprj, icorct)
if (icorct <= 0) goto 500
end if
end if
! THE NEXT GROUP OF STATEMENTS IS TO REDUCE PROJCT (IF POSSIBLE) IN CASE
! OF A FAILURE OF SOME KIND.
300 if (nfail < limfl) then
! PREPARE TO TRY ANOTHER ITERATION IN RKCON BY
! REDUCING PROJCT, AND MAKING SURE PROJCT IS NOT TOO SMALL.
Projct = Projct/(four + four)
if (Projct >= steplm) goto 400
if (iwarn <= 0) then
iwarn = 1
Projct = steplm
goto 400
end if
end if
! HERE RKCON COULD NOT REDUCE THE ERROR NORM AND WE RETURN WITH THE
! WARNING ISUCC=1.
Isucc = 1
return
! NOW RESET ACTDIF FOR THIS PROJCT.
400 do l = 1, mactrk
i = abs(Iact(l))
if (Icntyp(i) < 0) then
if (Icntyp(i) + 1 < 0) then
! HERE WE HAVE AN ACTIVE TYPE -2 CONSTRAINT, AND WE SET ACTDIF(L)=
! MIN (CONUP, ERROR(I)/PROJCT).
Actdif(l) = Error(i)/Projct
if (Actdif(l) > conup) Actdif(l) = conup
else
! HERE WE HAVE AN ACTIVE TYPE -1 CONSTRAINT.
Actdif(l) = Error(i)/Projct
end if
else if (Icntyp(i) == 0) then
! ICNTYP(I)=0 SHOULD NOT OCCUR HERE SINCE CONSTRAINT I WAS DECLARED
! TO BE ACTIVE IN RKSACT, BUT WE ACCOUNT FOR IT ANYWAY AS A PRECAUTION.
Actdif(i) = zero
else if (Icntyp(i) <= 1) then
! HERE WE HAVE AN ACTIVE TYPE 1 CONSTRAINT.
Actdif(l) = one + (Error(i) - enorm)/Projct
else
! HERE WE HAVE AN ACTIVE TYPE 2 CONSTRAINT.
Actdif(l) = one + (abs(Error(i)) - enorm)/Projct
end if
end do
! COPY THE FIRST NPARM ROWS OF PMAT FROM OLD PMAT TRANSPOSE STORED
! IN FUNTBL, THEN APPEND ACTDIF AS THE LAST ROW.
do j = 1, mactrk
do i = 1, Nparm
Pmat(i, j) = Funtbl(j, i)
end do
Pmat(npar1, j) = Actdif(j)
end do
goto 200
! PUT THE SEARCH DIRECTION VECTOR PARPRJ - PARAM INTO XRK.
500 do j = 1, Nparm
Xrk(j) = Parprj(j) - Param(j)
end do
! CALL SEARSL TO DO A LINE SEARCH IN DIRECTION XRK AND PUT THE RESULTING
! VECTOR IN PARSER. START WITH A PROJECTION FACTOR PROSEA=1.0.
! PARPRJ WILL BE USED TEMPORARILY AS A WORK VECTOR IN SEARSL.
prosea = one
! WE NOW WISH TO DETERMINE PRJLIM = THE SMALLER OF 1.0/SPCMN AND
! THE LARGEST VALUE OF PROSEA FOR WHICH THE LINEAR STANDARD CONSTRAINTS
! ARE SATISFIED FOR THE PARAMETER VECTOR PARAM+PROSEA*XRK. THIS
! WILL GIVE AN UPPER BOUND FOR LINE SEARCHING. NOTE THAT IN
! THEORY WE SHOULD HAVE PRJLIM >= 1.0 SINCE THE LINEAR STANDARD
! CONSTRAINTS SHOULD BE SATISFIED FOR PROSEA=0.0 AND PROSEA=1.0, BUT
! ROUNDOFF ERROR COULD AFFECT THIS A LITTLE. IF THERE ARE NO
! LINEAR STANDARD CONSTRAINTS, WE SET PRJLIM=1.0/SPCMN.
prjlim = one/spcmn
!*****INSERT TO MAKE SEARCHING LESS VIOLENT.
! PRJLIM=TWO
!*****END INSERT
if (Itypm1 > 0) then
! HERE WE HAVE AT LEAST ONE TYPE -1 CONSTRAINT, AND IF IOPTTH=1 WE
! CALL DERST TO PUT ALL THE STANDARD CONSTRAINT VALUES AND GRADIENTS
! INTO CONFUN(.,.).
if (ioptth > 0) then
! WE SET IPT=-1 TO TELL DERST TO COMPUTE STANDARD CONSTRAINTS ONLY.
ipt = -1
call me%derst(Ioptn, Nparm, Numgr, Pttbl, Iptb, Indm, Param, ipt, &
Work(ilc24), Work(ilc35), Iwork(ilc22), Confun)
end if
do i = 1, Numgr
if (Icntyp(i) + 1 == 0) then
ipt = i
! HERE WE HAVE A TYPE -1 CONSTRAINT AND IF IOPTTH=0 WE CALL DERST
! TO PUT THE CONSTRAINT VALUE AND GRADIENT INTO CONFUN(IPT,.).
if (ioptth <= 0) call me%derst(Ioptn, Nparm, Numgr, Pttbl, Iptb, &
Indm, Param, ipt, Work(ilc24), Work(ilc35), &
Iwork(ilc22), Confun)
! WE WISH TO HAVE SUMMATION (CONFUN(IPT,J+1)*(PARAM(J)+PROSEA*XRK(J)))
! + C(IPT) <= 0.0 FOR IPT=1,...,NUMGR, ICNTYP(IPT) = -1,
! WHERE THE IPTTH CONSTRAINT APPLIED TO PARAM SAYS
! SUMMATION (CONFUN(IPT,J+1)*PARAM(J)) + C(IPT) <= 0.0, SO C(IPT) IS
! THE CONSTANT TERM IN THE LEFT SIDE OF LINEAR CONSTRAINT IPT.
! THUS FOR I=1PT,...,NUMGR, ICNTYP(IPT) = -1, WE WANT PRJLIM*SS <=
! SSS, WHERE SS = SUMMATION (CONFUN(IPT,J+1)*XRK(J)) AND SSS = -C(IPT) -
! SUMMATION (CONFUN(IPT,J+1)*PARAM(J)) = -CONFUN(IPT,1).
ss = zero
do j = 1, Nparm
ss = ss + Confun(i, j + 1)*Xrk(j)
end do
! IF SS < 10.0*SPCMN THIS CONSTRAINT WILL NOT PUT A SIGNIFICANT
! RESTRICTION ON PROSEA.
if (ss >= tol2) then
! HERE SS >= 10.0*SPCMN AND WE COMPARE SSS/SS AGIANST PRJLIM.
quots = -Confun(i, 1)/ss
if (prjlim > quots) prjlim = quots
end if
end if
end do
end if
call me%searsl(Ioptn, Numgr, Nparm, prjlim, tol1, Xrk, Fun, Ifun, Pttbl, Iptb, &
Indm, Param, Error, Rchdwn, mactrk, Iact, Iphse, unit, Tolcon, &
Rchin, Itypm1, Itypm2, Iwork, Liwrk, Work, Lwrk, Err1, Parprj, &
prosea, emin, emin1, Parser, nsrch)
! COMPUTE THE PRINCIPAL ERROR NORM CHANGE ENCHG. ALSO COMPUTE ENC1, THE
! CHANGE IN THE PRINCIPAL ERROR NORM WITHOUT THE LINE SEARCH.
Enchg = emin - enorm
Enc1 = emin1 - enorm
! IF WE OBTAINED MORE THAN A TOL1 REDUCTION IN ENORM WE UPDATE
! PARAM AND CALL ERCMP1 TO UPDATE ERROR, AND RETURN WITH ISUCC=0
! INDICATING SUCCESS.
! OTHERWISE WE CHECK TO SEE IF WE HAVE REACHED THE RKCON ITERATION
! LIMIT, AND IF SO WE RETURN WITH ISUCC=1, INDICATING FAILURE.
if (Enchg + tol1 >= 0) goto 300
do j = 1, Nparm
Param(j) = Parser(j)
end do
call me%ercmp1(Ioptn, Nparm, Numgr, Fun, Ifun, Pttbl, Iptb, Indm, Param, 1, &
Iphse, Iwork, Liwrk, Confun, Icntyp, ipmax, ismax, Error)
end subroutine rkcon