refwl Subroutine

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

        implicit none

        integer, intent(in) :: Nparm
        integer, intent(in) :: Numgr
        integer, intent(in) :: Ncor
        integer, intent(in) :: Ndm
        integer :: Icor(Nparm + 1)
        integer :: Ixrct(2*Nparm)
        real(wp) :: Pmat1(Nparm + 1, Numgr)
        real(wp) :: Pmat(Nparm + 1, Numgr)
        real(wp) :: Save(Nparm)
        real(wp) :: Wpt(Nparm + 1)

        real(wp) :: aa, amax, fact, wrst, wrsto
        integer :: i, imax, itrct, itrlm, j, jmax, jstrt, &
                   k, kcol, kk, kp1, l, maxrs, n, nresl

        real(wp), parameter :: tole = spcmn

        ! COMPUTE MACHINE AND PRECISION DEPENDENT CONSTANTS.
        itrlm = 2
        itrct = 0
        nresl = 0
        n = Ndm + 1

        ! IF NCOR=0 WE HAVE NOTHING TO DO SO WE RETURN.
        if (Ncor > 0) then

            ! COPY COLUMN ICOR(L) OF PMAT WITH THE SIGN OF THE LAST ELEMENT REVERSED
            ! INTO COLUMN L OF THE WORK MATRIX PMAT1 FOR L=1,...,NCOR.
            do l = 1, Ncor
                j = Icor(l)
                do i = 1, Ndm
                    Pmat1(i, l) = Pmat(i, j)
                end do
                Pmat1(n, l) = -Pmat(n, j)
            end do

            ! NOW COLUMN REDUCE PMAT1.  NOTE THAT PMAT1 IS THE TRANSPOSE OF THE USUAL
            ! AUGMENTED MATRIX FOR SOLVING A LINEAR SYSTEM OF EQUATONS.
            ! THERE WILL BE AT MOST MAXRS = MIN(NDM,NCOR) RESOLVENTS.
            maxrs = Ncor
            if (Ndm < maxrs) maxrs = Ndm
            column_reduce: do k = 1, maxrs
                ! SEARCH FOR THE INDICES IMAX AND JMAX WITH 1 <= IMAX <= NDM, 1 <=
                ! JMAX <= NCOR, PMAT1(IMAX,JMAX) IS NOT IN THE ROW OR COLUMN OF ANY
                ! OTHER RESOLVENT (I.E. PIVOT), AND ABS(PMAT1(IMAX,JMAX)) IS MAXIMIZED.
                ! WE USE THE VECTOR IXRCT TO SAVE THE RESOLVENT POSITIONS TO SAVE SPACE.
                jstrt = 0
                outer: do j = 1, Ncor
                    if (nresl > 0) then
                        do l = 1, nresl
                            if (j == Ixrct(2*l)) cycle outer
                        end do
                    end if
                    ! HERE THERE IS NO EARLIER RESOLVENT IN COLUMN J.
                    inner: do i = 1, Ndm
                        if (nresl > 0) then
                            do l = 1, nresl
                                if (i == Ixrct(2*l - 1)) cycle inner
                            end do
                        end if
                        ! HERE THERE IS NO EARLIER RESOLVENT IN ROW I.
                        aa = abs(Pmat1(i, j))
                        if (jstrt <= 0) then
                            jstrt = 1
                        else if (aa <= amax) then
                            cycle inner
                        end if
                        amax = aa
                        imax = i
                        jmax = j
                    end do inner
                end do outer

                ! IF THE ABSOLUTE VALUE OF THIS RESOLVENT IS VERY SMALL WE DO NOT ATTEMPT
                ! ANY FURTHER COLUMN OPERATIONS.
                if (amax < tole) exit column_reduce

                ! INCREMENT NRESL AND PUT THE LOCATION OF THE NRESLTH RESOLVENT IN
                ! (IXRCT(2*L-1),IXRCT(2*L)).
                nresl = nresl + 1
                Ixrct(2*nresl - 1) = imax
                Ixrct(2*nresl) = jmax

                ! NOW ELIMINATE WPT(IMAX) FROM THOSE COLUMNS WHICH DO NOT CONTAIN ANY OF
                ! THE RESOLVENTS FOUND SO FAR (INCLUDING THE PRESENT RESOLVENT).
                outer2: do j = 1, Ncor
                    do l = 1, nresl
                        if (j == Ixrct(2*l)) cycle outer2
                    end do
                    ! HERE COLUMN J DOES NOT CONTAIN ANY OF THE RESOLVENTS FOUND SO FAR, AND
                    ! WE COMPUTE THE FACTOR FOR THE COLUMN OPERATION NEEDED TO ZERO OUT
                    ! PMAT1(IMAX,J) (ALTHOUGH WE DO NOT ACTUALLY WRITE IN THE ZERO).
                    fact = Pmat1(imax, j)/Pmat1(imax, jmax)
                    ! NOW DO THE OPERATION IN COLUMN J FOR ALL ROWS NOT CONTAINING A
                    ! RESOLVENT.  THE ELEMENTS IN THIS COLUMN IN THE ROWS WHICH CONTAIN AN
                    ! EARLIER (OR PRESENT) RESOLVENT WILL NOT BE NEEDED LATER.
                    inner2: do i = 1, n
                        do l = 1, nresl
                            if (i == Ixrct(2*l - 1)) cycle inner2
                        end do
                        Pmat1(i, j) = Pmat1(i, j) - fact*Pmat1(i, jmax)
                    end do inner2
                end do outer2
            end do column_reduce
            ! END OF COLUMN REDUCTION OF PMAT1.

            ! IF NRESL=0 THEN ALL THE ELEMENTS IN PMAT1 FOR 1 <= I <= NDM AND
            ! 1 <= J <= NCOR WERE VERY SMALL IN ABSOLUTE VALUE, AND THERE IS
            ! NOTHING WE CAN DO, SO WE RETURN.
            if (nresl <= 0) return

        else
            return
        end if

        do
            ! NOW DO BACK SUBSTITUTION TO COMPUTE, FOR K=NRESL,...,1,
            ! WPT(IXRCT(2*K-1)) = (PMAT1(NDM+1,IXRCT(2*K)) - SUMMATION(
            ! PMAT1(I,IXRCT(2*K))*WPT(I), FOR I = 1,...,NDM, I /= IXRCT(2*L-1)
            ! FOR ANY L=1,...,K))/PMAT1(IXRCT(2*K-1),IXRCT(2*K)).  IF WE ARE IN AN
            ! ITERATIVE REFINEMENT STEP WE WISH TO CONSIDER WPT(I) (WHICH IS THEN
            ! JUST A CORRECTION TO WPT(I)) = 0.0 IF I CORRESPONDS TO NO RESOLVENT
            ! (SINCE THE VALUE OF SUCH WPT(I) IN SAVE SHOULD NOT CHANGE) SO WE OMIT
            ! THE CORRESPONDING TERMS IN THE SUMMATION ABOVE.
            do kk = 1, nresl
                k = nresl - kk + 1
                imax = Ixrct(2*k - 1)
                jmax = Ixrct(2*k)
                Wpt(imax) = Pmat1(n, jmax)
                iloop: do i = 1, Ndm
                    do l = 1, k
                        if (i == Ixrct(2*l - 1)) cycle iloop
                    end do
                    ! HERE ROW I CONTAINS NO EARLIER (OR PRESENT) RESOLVENTS.
                    if (itrct > 0) then
                        if (k < nresl) then
                            ! HERE WE ARE DOING ITERATIVE REFINEMENT, K < NRESL, AND I /=
                            ! IXRCT(2*L-1) FOR L=1,...,K.  WE WILL USE THE TERM CORRESPONDING TO
                            ! WPT(I) IFF I = IXRCT(2*L-1) FOR SOME L = K+1,...,NRESL.
                            kp1 = k + 1
                            do l = kp1, nresl
                                if (i == Ixrct(2*l - 1)) then
                                    Wpt(imax) = Wpt(imax) - Pmat1(i, jmax)*Wpt(i)
                                    cycle iloop
                                end if
                            end do
                        end if
                    else
                        Wpt(imax) = Wpt(imax) - Pmat1(i, jmax)*Wpt(i)
                    end if
                end do iloop
                Wpt(imax) = Wpt(imax)/Pmat1(imax, jmax)
            end do
            ! END OF BACK SUBSTITUTION.

            ! IF ITRCT IS POSITIVE THEN WPT WILL CONTAIN ONLY AN ITERATIVE
            ! REFINEMENT CORRECTION IN THOSE POSITIONS CORRESPONDING TO RESOLVENTS
            ! AND WE ADD THIS TO SAVE TO GET THE TRUE WPT.
            if (itrct > 0) then
                do i = 1, Ndm
                    do l = 1, nresl
                        if (i == Ixrct(2*l - 1)) then
                            Wpt(i) = Wpt(i) + Save(i)
                            exit
                        end if
                    end do
                end do
            end if

            ! NOW COMPUTE THE RESIDUAL AND PUT IT INTO PMAT1(NDM+1,.).
            do k = 1, Ncor
                ! COMPUTE THE COLUMN INDEX KCOL IN PMAT CORRESPONDING TO COLUMN K IN
                ! PMAT1.
                kcol = Icor(k)
                Pmat1(n, k) = -Pmat(n, kcol)
                do i = 1, Ndm
                    Pmat1(n, k) = Pmat1(n, k) - Pmat(i, kcol)*Wpt(i)
                end do
            end do

            ! COMPUTE THE WORST ABSOLUTE VALUE OF THE RESIDUAL ELEMENTS.
            do k = 1, Ncor
                aa = abs(Pmat1(n, k))
                if (k > 1) then
                    if (aa <= wrst) cycle
                end if
                wrst = aa
            end do

            if (itrct <= 0) then
                wrsto = wrst
            else if (wrst > wrsto) then
                ! HERE ITRCT > 0 AND WRST > WRSTO, SO WE GO BACK TO THE PREVIOUS
                ! WPT AND RETURN.
                wrst = wrsto
                do i = 1, Ndm
                    Wpt(i) = Save(i)
                end do
                return
            end if

            if (itrct >= itrlm) return

            ! HERE ITRCT < ITRLM AND WE INCREMENT ITRCT AND SET UP FOR THE ITRCTTH
            ! ITERATIVE REFINEMENT STEP.
            itrct = itrct + 1
            ! COPY WPT INTO SAVE.
            do i = 1, Ndm
                Save(i) = Wpt(i)
            end do

        end do

    end subroutine refwl