setu1 Subroutine

    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)

        implicit none

        class(conmax_solver), intent(inout) :: me
        real(wp) :: actlim, bndfud, Bndlgt, Cofbnd, Confun, enorm, &
                    Error, Fun, Funtbl, grdlgt, Param, &
                    Pttbl, Rchdwn, Rchin, rchind, rt, stfudg, sum
        real(wp) :: v, Work
        integer :: i, Iact, Iact1, Icntyp, Ifun, ii, ilc22, ilc24, &
                   Indm, Ioptn, ioptth, Iphse, ipt, Iptb, &
                   Iwork, Iyrct, j, jj, k
        integer :: kk, l, Liwrk, Lwrk, m, mact, Mact1, mm1, mp1, &
                   npar1, npar2, Nparm, Numgr, Numin

        dimension Pttbl(Iptb, Indm), Fun(Ifun), Funtbl(Numgr, Nparm + 1), &
            Cofbnd(Nparm), Param(Nparm), Error(Numgr + 3), &
            v(Numgr + 2*Nparm + 1, Nparm + 2), Iact(Numgr), Iact1(Numgr), &
            Iyrct(Numgr + 2*Nparm), Icntyp(Numgr), &
            Confun(Numgr, Nparm + 1), Iwork(Liwrk), Work(Lwrk)

        ilc22 = iloc(22, Nparm, Numgr)
        ilc24 = iloc(24, Nparm, Numgr)
        npar1 = Nparm + 1
        npar2 = Nparm + 2
        ioptth = (Ioptn - (Ioptn/100000)*100000)/10000

        ! THE LINEARIZED PROBLEM REPLACES THE APPROXIMATING FUNCTION BY ITS
        ! FIRST ORDER TAYLOR SERIES, SO FUN(I)-(APPROXIMATING FUNCTION)(I) IS
        ! REPLACED BY ERROR(I)-(SUMMATION OF COEFFICIENT CHANGES TIMES PARTIAL
        ! DERIVATIVES OF THE APPROXIMATING FUNCTION WITH RESPECT TO THOSE
        ! COEFFICIENTS) IF ICNTYP(I)=2, AND IF ICNTYP(I)=1 WE REPLACE THE LEFT
        ! SIDE OF CONSTRAINT I BY ERROR(I)+(SUMMATION OF COEFFICIENT CHANGES TIMES
        ! PARTIAL DERIVATIVES OF THE LEFT SIDE OF CONSTRAINT I).
        ! V AND M ARE THE OUTPUT QUANTITIES.  M WILL KEEP TRACK OF THE NUMBER
        ! OF CONSTRAINTS IN THE LP PROBLEM TO BE SOLVED BY SLNPRO.
        m = 0
        enorm = Error(Numgr + 1)
        stfudg = one/ten

        ! COMPUTE THE LENGTH OF THE LONGEST X VECTOR SATISFYING THE COEFFICIENT
        ! CHANGE BOUNDS.
        sum = zero
        do j = 1, Nparm
            sum = sum + (Cofbnd(j))**2
        end do
        Bndlgt = sqrt(sum)
        bndfud = stfudg*Bndlgt

        ! WE WILL SAY A PRIMARY CONSTRAINT IS ACTIVE IF ERROR(I) (OR ABS(ERROR(I
        ! IF ICNTYP(I)=2) >= ENORM-RCHDWN*BNDLGT.
        actlim = enorm - Rchdwn*Bndlgt

        ! WE WILL SAY A TYPE -2 CONSTRAINT IS ACTIVE IF ERROR(I) >= -RCHIND.
        rchind = Rchin*Bndlgt

        if (Numin <= 0) then

            ! HERE NUMIN=0, SO WE WILL FIRST COMPUTE A NEW SET OF ACTIVE INDICES,
            ! THEN PUT THE FUNCTION VALUES AND GRADIENTS FOR THESE INDICES IN
            ! FUNTBL, WHERE THEY WILL REMAIN THROUGHOUT THIS CALL TO SLPCON.
            do i = 1, Numgr
                if (Icntyp(i) < 0) then
                    ! HERE ICNTYP(I) < 0 AND WE WILL DECLARE THE CONSTRAINT TO BE ACTIVE IF
                    ! AND ONLY IF ICNTYP(I)=-1, OR ICNTYP(I)=-2 AND ERROR(I) >= -RCHIND.
                    if (Icntyp(i) + 1 < 0) then
                        if (Error(i) + rchind < 0) cycle
                    end if
                else if (Icntyp(i) == 0) then
                    cycle
                else if (Icntyp(i) <= 1) then
                    ! HERE ICNTYP(I)=1 AND WE WILL DECLARE THE CONSTRAINT TO BE +ACTIVE IF AND
                    ! ONLY IF ERROR(I) >= ACTLIM.
                    if (Error(i) < actlim) cycle
                    ! HERE ICNTYP(I)=2 AND WE WILL DECLARE THE CONSTRAINT TO BE +ACTIVE IF AND
                    ! ONLY IF ERROR(I) >= ACTLIM OR -ACTIVE IF AND ONLY IF ERROR(I)  <=
                    ! -ACTLIM.
                else if (Error(i) < actlim) then
                    if (Error(i) + actlim <= 0) then
                        ! DECLARE CONSTRAINT I TO BE -ACTIVE.
                        m = m + 1
                        Iact(m) = -i
                    end if
                    cycle
                end if
                ! DECLARE CONSTRAINT I TO BE (+)ACTIVE.
                m = m + 1
                Iact(m) = i
            end do
            mact = m

            ! NOW PUT ACTIVE VALUES AND GRADIENTS IN FUNTBL.
            if (ioptth <= 0) then
                ! HERE IOPTTH=0 AND WE CALL DERST FOR EACH ACTIVE CONSTRAINT.
                do l = 1, mact
                    i = abs(Iact(l))
                    ipt = i
                    ! CALL DERST TO COMPUTE BOTH FUNCTION AND GRADIENT VALUES.
                    call me%derst(Ioptn, Nparm, Numgr, Pttbl, Iptb, Indm, Param, ipt, &
                                  Work(ilc24), v, Iwork(ilc22), Confun)
                    ! COPY THE VALUES FOR CONSTRAINT I INTO FUNTBL.
                    do j = 1, npar1
                        Funtbl(i, j) = Confun(i, j)
                    end do
                end do
                goto 300
            else
                ! HERE IOPTTH=1 AND ONLY ONE DERST CALL IS NEEDED.
                ! IF IPHSE < 0 OR NO ICNTYP(L) IS POSITIVE, SET IPT=-1 TO TELL DERST
                ! TO COMPUTE STANDARD CONSTRAINTS ONLY, WHILE OTHERWISE SET IPT=0 TO
                ! TELL DERST TO COMPUTE ALL CONSTRAINTS.
                if (Iphse >= 0) then
                    do l = 1, Numgr
                        if (Icntyp(l) > 0) goto 100
                    end do
                end if
                ipt = -1
                goto 200
            end if
100         ipt = 0
        else
            ! HERE NUMIN IS NOT 0, AND WE WILL KEEP THE OLD ACTIVE CONSTRAINT SET
            ! AND FOREGO RECOMPUTING FUNCTION VALUES AND GRADIENTS.
            mact = Mact1
            m = mact
            do l = 1, mact
                Iact(l) = Iact1(l)
            end do
            goto 300
        end if
200     call me%derst(Ioptn, Nparm, Numgr, Pttbl, Iptb, Indm, Param, ipt, Work(ilc24), &
                      v, Iwork(ilc22), Confun)
        ! COPY THE ACTIVE FUNCTION AND GRADIENT VALUES INTO FUNTBL.
        do l = 1, mact
            i = abs(Iact(l))
            do j = 1, npar1
                Funtbl(i, j) = Confun(i, j)
            end do
        end do

        ! NOW SET UP THE ACTIVE CONSTRAINTS IN V FOR SLNPRO.
300     do l = 1, mact
            i = abs(Iact(l))
            if (Icntyp(i) < 0) then
                if (Icntyp(i) + 1 < 0) then
                    ! HERE ICNTYP(I)=-2 AND WE FIRST COMPUTE THE LENGTH OF THE GRADIENT.
                    sum = zero
                    do j = 1, Nparm
                        sum = sum + (Funtbl(i, j + 1))**2
                    end do
                    grdlgt = sqrt(sum)
                    ! NOW SET UP A CONSTRAINT OF THE FORM GRADIENT.CHANGE  <=
                    ! -MIN(1.0,CONSTRAINT VALUE)*BNDFUD*GRDLGT, SO IF GRDLGT > 0.0 WE
                    ! HAVE (-GRADIENT/GRDLGT).(CHANGE/BNDLGT) >= STFUDG*MIN(1.0,
                    ! CONSTRAINT VALUE).
                    do j = 1, Nparm
                        v(l, j) = Funtbl(i, j + 1)
                    end do
                    v(l, npar1) = zero
                    rt = Error(i)
                    if (rt > one) rt = one
                    v(l, npar2) = -rt*bndfud*grdlgt
                else
                    ! HERE ICNTYP(I)=-1 AND WE SET UP A CONSTRAINT OF THE FORM
                    ! GRADIENT.CHANGE <= -CONSTRAINT VALUE.
                    do j = 1, Nparm
                        v(l, j) = Funtbl(i, j + 1)
                    end do
                    v(l, npar1) = zero
                    v(l, npar2) = -Error(i)
                end if
            else if (Icntyp(i) /= 0) then
                if (Icntyp(i) <= 1) then
                    ! HERE ICNTYP(I)=1 AND WE SET UP A CONSTRAINT OF THE FORM
                    ! GRADIENT.CHANGE - W <= -CONSTRAINT VALUE.
                    do j = 1, Nparm
                        v(l, j) = Funtbl(i, j + 1)
                    end do
                    v(l, npar1) = -one
                    v(l, npar2) = -Error(i)
                else if (Iact(l) <= 0) then
                    ! HERE ICNTYP(I)=2 AND IACT(L) < 0, AND WE SET UP A CONSTRAINT OF THE
                    ! FORM GRADIENT.CHANGE - W <= FUN - CONSTRAINT VALUE.
                    do j = 1, Nparm
                        v(l, j) = Funtbl(i, j + 1)
                    end do
                    v(l, npar1) = -one
                    v(l, npar2) = Error(i)
                else
                    ! HERE ICNTYP(I)=2 AND IACT(L) > 0, AND WE SET UP A CONSTRAINT OF THE
                    ! FORM -GRADIENT.CHANGE - W <= -(FUN - CONSTRAINT VALUE).
                    do j = 1, Nparm
                        v(l, j) = -Funtbl(i, j + 1)
                    end do
                    v(l, npar1) = -one
                    v(l, npar2) = -Error(i)
                end if
            end if
        end do
        ! SET THE CONSTRAINTS OF THE FORM -X(J) <= COFBND(J) AND
        ! X(J) <= COFBND(J).
        do j = 1, Nparm
            m = m + 2
            mm1 = m - 1
            do k = 1, npar1
                v(mm1, k) = zero
                v(m, k) = zero
            end do
            v(mm1, j) = -one
            v(m, j) = one
            v(mm1, npar2) = Cofbnd(j)
            v(m, npar2) = Cofbnd(j)
        end do

        ! NOW SET THE BOTTOM ROW.  TO MINIMIZE W = X(NPARM+1) WE MAXIMIZE -W.
        mp1 = m + 1
        do j = 1, npar2
            v(mp1, j) = zero
        end do
        v(mp1, npar1) = one

        ! THIS SECTION ADJUSTS IYRCT TO EITHER TELL SLNPRO TO DO THE INITIAL
        ! EXCHANGES STRICTLY ACCORDING TO A PIVOTING STRATEGY (BY SETTING
        ! IYRCT(1)=-1) OR TO SPECIFY AN INITIAL VERTEX FOR SLNPRO, NAMELY THE
        ! VERTEX CORRESPONDING TO THE LAST LINEAR PROGRAMMING SOLUTION.
        ! IF IYRCT(1) IS -1 ALREADY WE DO NOT ATTEMPT TO SPECIFY A VERTEX, BUT
        ! WE STORE MACT IN MACT1 AND IACT IN IACT1 FOR POSSIBLE LATER USE.
        if (Iyrct(1) >= 0) then
            ! HERE IYRCT(1) /= -1, AND WE CONSIDER THE PRESENT ENTRIES IN IYRCT
            ! ONE BY ONE.
            do j = 1, npar1
                jj = Iyrct(j)
                if (jj <= Mact1) then
                    ! HERE ENTRY J OF IYRCT CORRESPONDS TO A FORMER ACTIVE CONSTRAINT AT
                    ! SOME POINT abs(KK), WHERE THE SIGN OF KK WILL INDICATE WHETHER THE
                    ! CONSTRAINT WAS +ACTIVE OR -ACTIVE.
                    kk = Iact1(jj)
                    ! WE NOW CHECK TO SEE IF THIS FORMER ACTIVE CONSTRAINT IS STILL
                    ! ACTIVE WITH THE SAME SIGN.  IF SO, WE RESET IYRCT(J) TO THE PRESENT
                    ! NUMBER OF THIS CONSTRAINT, AND IF NOT (WHICH WILL OCCUR IFF THE K
                    ! LOOP BELOW IS COMPLETED), WE WILL NOT TRY TO DETERMINE A VERTEX, SO
                    ! WE WILL SET IYRCT(1)=-1 AND LEAVE THE J LOOP.
                    do k = 1, mact
                        if (kk == Iact(k)) then
                            Iyrct(j) = k
                            goto 350
                        end if
                    end do
                    Iyrct(1) = -1
                    goto 500
                else
                    ! HERE ENTRY J OF IYRCT CORRESPONDS TO A CONSTRAINT BEYOND THE ACTIVE
                    ! POINT CONSTRAINTS, AND WE ADJUST IYRCT(J) BY THE DIFFERENCE OF THE
                    ! PRESENT AND FORMER NUMBER OF ACTIVE CONSTRAINTS.
                    Iyrct(j) = Iyrct(j) + mact - Mact1
                end if
350         end do
            ! WE HAVE NOW FILLED IN IYRCT(1),...,IYRCT(NPARM+1) WITH DISTINCT
            ! POSITIVE INTEGERS BETWEEN 1 AND M, AND WE FILL IN THE REST OF IYRCT
            ! SO THAT IYRCT WILL CONTAIN A PERMUTATION OF 1,...,M.  TO BE CONSISTENT
            ! WITH SLNPRO WE PUT IYRCT(NPARM+2),...,IYRCT(M) IN DECREASING ORDER.
            l = npar1
            do i = 1, m
                ii = m - i + 1
                ! SKIP II IF IT IS ALREADY IN IYRCT.
                do j = 1, npar1
                    if (ii == Iyrct(j)) goto 400
                end do
                l = l + 1
                Iyrct(l) = ii
400         end do
        end if

        ! SAVE MACT IN MACT1 AND IACT IN IACT1 AND RETURN.
500     Mact1 = mact
        do j = 1, mact
            Iact1(j) = Iact(j)
        end do

    end subroutine setu1