conenr Subroutine

    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)

        implicit none

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

        real(wp) :: amax, amin, cjj, diff, dmax, &
                    dp, dsq, omt, pdotj, quot, theta, tst
        integer :: i, icoro, ihouse, ii, ilc01, ilc03, ilc04, &
                   ilc09, ilc23, ilc34, itst1, j, jj, jmax, jmin, &
                   kntsl, l, limsl, ll, mincf, mp1, &
                   ncoro, ndm

        real(wp), parameter :: tolel = ten*ten*spcmn
        real(wp), parameter :: z1 = ten*tolel
        real(wp), parameter :: z2 = ten*z1
        real(wp), parameter :: z3 = ten*z1

        ilc01 = iloc(1, Nparm, Numgr)
        ilc03 = iloc(3, Nparm, Numgr)
        ilc04 = iloc(4, Nparm, Numgr)
        ilc09 = iloc(9, Nparm, Numgr)
        ilc23 = iloc(23, Nparm, Numgr)
        ilc34 = iloc(34, Nparm, Numgr)
        kntsl = 0
        limsl = 100
        mp1 = m + 1
        ndm = n - 1
        Nmaj = 0
        Nmin = 0
        Jflag = 0
        itst1 = 0
        ncoro = -1

        if (Istrt1 > 0) goto 200

        ! HERE ISTRT1=0 SO WE START FROM SCRATCH.  FIND THE INDEX JMAX FOR
        ! WHICH (P(J).R)/SQRT(P(J).P(J)) IS MAXIMIZED FOR P(J).P(J) > Z1.
100     amax = zero
        jmax = 0
        do j = 1, m
            do i = 1, n
                Vec(i) = Pmat1(i, j)
            end do
            pdotj = dotprd(n, Vec, Vec, Nparm)
            if (pdotj > z1) then
                quot = dotprd(n, Vec, r, Nparm)/sqrt(pdotj)
                if (quot > amax) then
                    amax = quot
                    jmax = j
                end if
            end if
        end do
        if (jmax > 0) then
            ! IF AMAX IS NOT SIGINFICANTLY POSITIVE WE PROCEED AS IF IT WERE ZERO.
            if (amax*sqrt(ndm + one) > tolel) then
                ! HERE WE FOUND THE RAY CLOSEST TO R AND WE COMPLETE THE
                ! INITIALIZATION BY SETTING NCOR=1, ICOR(1)=JMAX, AND COEF(JMAX)=1.0
                ! (WITH ALL OTHER ENTRIES OF COEF EQUAL TO ZERO).
                Ncor = 1
                Icor(1) = jmax
                do i = 1, m
                    Coef(i) = zero
                end do
                Coef(jmax) = one
                goto 200
            end if
        end if

        ! HERE THERE WERE NO VECTORS P(J) WHICH HAVE BOTH LENGTH SQUARED
        ! GREATER THAN Z1 AND ANGLE WITH R SIGNIFICANTLY LESS THAN 90 DEGREES,
        ! AND WE SET NCOR=0, PTNR=THE ZERO VECTOR, COEF=THE ZERO VECTOR, DIST=
        ! THE LENGTH OF R, AND WE RETURN.
        Ncor = 0
        do i = 1, n
            Ptnr(i) = zero
        end do
        do j = 1, m
            Coef(j) = zero
        end do
        Dist = sqrt(dotprd(n, r, r, Nparm))
        return

        ! SET PTNR TO THE CURRENT NEAREST POINT.  FIRST ZERO IT OUT.
200     do i = 1, n
            Ptnr(i) = zero
        end do
        if (Ncor > 0) then
            ! HERE NCOR > 0 AND WE SET PTNR=SUMMATION(COEF(J)*P(J)).
            do j = 1, Ncor
                jj = Icor(j)
                cjj = Coef(jj)
                do i = 1, n
                    Ptnr(i) = Ptnr(i) + cjj*Pmat1(i, jj)
                end do
            end do
        end if

        ! PUT PTNR-R INTO PTNRR AND COMPUTE THE DISTANCE FROM PTNR TO R.
        do i = 1, n
            Ptnrr(i) = Ptnr(i) - r(i)
        end do
        dsq = dotprd(n, Ptnrr, Ptnrr, Nparm)
        Dist = sqrt(dsq)

        ! NOW CHECK OPTIMALITY.
        ! FIRST SEE WHETHER THE HYPERPLANE THROUGH PTNR PERPENDICULAR TO
        ! PTNR-R PASSES THROUGH THE ORIGIN.  IF NCOR=0 THIS WILL
        ! AUTOMATICALLY BE TRUE SINCE THEN PTNR IS THE ORIGIN.  IF IT IS NOT
        ! TRUE WE GO DOWN TO SOLVE FOR A NEW NEAREST POINT IN THE SUBSPACE
        ! DETERMINED BY THE CURRENT ICOR.
        if (Ncor > 0) then
            tst = dotprd(n, Ptnr, Ptnrr, Nparm)
            if (abs(tst) >= z1) goto 300
        end if
        ! HERE THE HYPERPLANE ROUGHLY PASSES THROUGH THE ORIGIN, AND WE
        ! CHECK WHETHER ALL P(J) VECTORS ARE ROUGHLY SEPARATED FROM R BY IT.
        ! PUT THE MINIMUM OF (PTNR-R).(P(J)-R) IN AMIN AND THE INDEX WHERE IT
        ! IS ACHIEVED IN JMIN.
        do i = 1, n
            Vec(i) = Pmat1(i, 1) - r(i)
        end do
        jmin = 1
        amin = dotprd(n, Ptnrr, Vec, Nparm)
        if (m > 1) then
            do j = 2, m
                do i = 1, n
                    Vec(i) = Pmat1(i, j) - r(i)
                end do
                dp = dotprd(n, Ptnrr, Vec, Nparm)
                if (dp < amin) then
                    amin = dp
                    jmin = j
                end if
            end do
        end if

        ! FOR TESTING PURPOSES COMPUTE THE MAXIMUM OF THE SQUARES OF THE
        ! LENGTHS OF THE DISTANCES CONSIDERED.
        do i = 1, n
            Vec(i) = Pmat1(i, jmin) - r(i)
        end do
        dmax = dotprd(n, Vec, Vec, Nparm)
        if (Ncor > 0) then
            do j = 1, Ncor
                jj = Icor(j)
                do i = 1, n
                    Vec(i) = Pmat1(i, jj) - r(i)
                end do
                dp = dotprd(n, Vec, Vec, Nparm)
                if (dp > dmax) dmax = dp
            end do
        end if
        ! DO THE TEST.  IF IT IS SUCCESSFUL, THEN WE HAVE (APPROXIMATE)
        ! OPTIMALITY AND WE RETURN.
        if (amin - dsq + z1*dmax < 0) then
            ! HERE PTNR IS NOT OPTIMAL.  AS A CHECK AGAINST BLUNDERS WE MAKE SURE
            ! NCOR < N AND JMIN IS NOT IN ICOR.
            if (Ncor > 0) then
                if (Ncor < n) then
                    do l = 1, Ncor
                        if (jmin == Icor(l)) goto 220
                    end do
                    goto 250
                end if
                ! HERE WE HAVE BLUNDERED SO WE SET JFLAG=1 AS A WARNING, COMPUTE DIST,
                ! AND RETURN.  FIRST TRY FROM SCRATCH IF THIS HAS NOT BEEN DONE.
220             if (Istrt1 + Jflag <= 0) then
                    Jflag = 1
                    ! write(nwrit,'(A)') '*****JFLAG IS 1 IN CONENR'
                    return
                else
                    Jflag = -1
                    kntsl = 0
                    goto 100
                end if
            end if
            ! HERE PTNR IS NOT OPTIMAL, NCOR < N, AND JMIN IS NOT IN ICOR.
            ! WE INCREMENT THE MAJOR CYCLE COUNTER AND ADD P(JMIN).
250         Nmaj = Nmaj + 1
            Ncor = Ncor + 1
            Icor(Ncor) = jmin
            Coef(jmin) = zero
        else
            return
        end if

        ! CHECK TO SEE IF WE HAVE SOLVED THE SYSTEM BELOW LIMSL TIMES ALREADY,
        ! AND IF SO, SET JFLAG=6 AS A WARNING AND RETURN (BUT
        ! TRY FROM SCRATCH BEFORE GIVING UP IF THIS HAS NOT ALREADY BEEN DONE).
300     if (kntsl < limsl) then

            ! CHECK TO SEE IF NCOR AND THE LAST ELEMENT IN ICOR ARE UNCHANGED FROM THE
            ! PREVIOUS HOUSE CALL (HA HA), WHICH INDICATES FAILURE.  NOTE THAT HERE WE
            ! MUST HAVE NCOR > 0.
            if (Ncor /= ncoro) then
                ncoro = Ncor
            else if (Icor(Ncor) == icoro) then
                ! HERE WE HAVE CYCLING AND WE SET JFLAG=2 AS A WARNING AND RETURN.  FIRST
                ! TRY FROM SCRATCH IF THIS HAS NOT BEEN DONE.
                if (Istrt1 + Jflag <= 0) then
                    Jflag = 2
                    return
                else
                    Jflag = -1
                    kntsl = 0
                    goto 100
                end if
            end if
            icoro = Icor(Ncor)
            kntsl = kntsl + 1

            ! NOW WE SOLVE THE SYSTEM PICOR*VEC = R IN THE LEAST SQUARES
            ! SENSE FOR THE COEFFICIENT VECTOR VEC (RELATIVE TO
            ! ICOR) OF THE NEAREST POINT TO R IN THE SUBSPACE SPANNED BY
            ! P(ICOR(1)),...,P(ICOR(NCOR)), WHERE P(ICOR) IS THE N X NCOR MATRIX
            ! WHOSE COLUMNS ARE THESE VECTORS.
            ! NOW FILL IN PICOR AND CALL HOUSE TO COMPUTE VEC.
            do j = 1, Ncor
                jj = Icor(j)
                do i = 1, n
                    Picor(i, j) = Pmat1(i, jj)
                end do
            end do

            call house(n, Ncor, Picor, r, Iwork(ilc23), Nparm, Work(ilc01), &
                       Work(ilc04), Work(ilc09), Work(ilc34), Work(ilc03), Vec, &
                       ihouse)

            ! IF HOUSE FAILS (INDICATED BY IHOUSE=1) WE SET JFLAG=3 AS A
            ! WARNING AND RETURN.  FIRST TRY FROM SCRATCH IF THIS HAS NOT BEEN DONE.
            if (ihouse <= 0) then
                ! CHECK TO SEE IF ALL THE COEFFICIENTS IN VEC ARE > Z2, AND IF SO,
                ! PUT VEC INTO COEF AND GO BACK TO COMPUTE PTNR.  THE COEFFICIENTS IN
                ! COEF NOT CORRESPONDING TO THOSE IN VEC WILL REMAIN EQUAL TO ZERO.
                do i = 1, Ncor
                    if (Vec(i) <= z2) goto 350
                end do
                do i = 1, Ncor
                    ii = Icor(i)
                    Coef(ii) = Vec(i)
                end do
                goto 200
            else if (Istrt1 + Jflag <= 0) then
                Jflag = 3
                return
            else
                Jflag = -1
                kntsl = 0
                goto 100
            end if

            ! HERE SOME ELEMENT OF VEC IS <= Z2.  COMPUTE THETA=MIN(1.0, MIN(
            ! COEF(ICOR(I))/(COEF(ICOR(I))-VEC(I)): COEF(ICOR(I))-VEC(I) > Z3)).
350         theta = one
            do l = 1, Ncor
                ll = Icor(l)
                diff = Coef(ll) - Vec(l)
                if (diff > z3) then
                    quot = Coef(ll)/diff
                    if (quot < theta) theta = quot
                end if
            end do
            ! COMPUTE THE NEW COEF AS (1.0-THETA)*COEF+THETA*VEC.
            omt = one - theta
            do l = 1, Ncor
                ll = Icor(l)
                Coef(ll) = omt*Coef(ll) + theta*Vec(l)
            end do
            ! COMPUTE THE INDEX MINCF (RELATIVE TO ICOR) OF THE SMALLEST ELEMENT OF
            ! COEF AND SET ALL ELEMENTS OF COEF WHICH ARE <= Z2 TO ZERO.
            mincf = 0
            amin = z2
            do i = 1, Ncor
                ii = Icor(i)
                if (Coef(ii) <= z2) then
                    if (Coef(ii) <= amin) then
                        amin = Coef(ii)
                        mincf = i
                    end if
                    Coef(ii) = zero
                end if
            end do

            if (mincf <= 0) then
                ! HERE MINCF=0 AND AN UNLIKELY BLUNDER HAS OCCURRED.  THIS MUST BE DUE TO
                ! ROUNDOFF ERROR SINCE IN THEORY (NEW) COEF(ICOR(I)) MUST BE <= Z2
                ! FOR SOME I=1,...,NCOR, WHICH MAKES MINCF > 0 IN THE LAST LOOP.
                ! TO SEE THIS, FIRST NOTE THAT FOR SOME IBAR=1,...,NCOR, VEC(IBAR) MUST
                ! BE <= Z2 SINCE OTHERWISE WE WOULD NOT BE HERE.  BY ITS DEFINITION,
                ! THETA MUST BE <= 1.0.  IF THETA = 1.0, THEN (NEW) COEF(ICOR(IBAR))
                ! = (1.0 - THETA)*(OLD) COEF(ICOR(IBAR)) + THETA*VEC(IBAR) = VEC(IBAR)
                ! <= Z2.  IF ON THE OTHER HAND THETA < 1.0, THEN FOR SOME ISTAR=1,
                ! ...,ICOR WE HAVE (OLD) COEF(ICOR(ISTAR)) - VEC(ISTAR) >= Z3 AND
                ! THETA = (OLD) COEF(ICOR(ISTAR))/((OLD) COEF(ICOR(ISTAR)) - VEC(ISTAR)),
                ! SO (NEW) COEF(ICOR(ISTAR)) = (1.0 - THETA)*(OLD) COEF(ICOR(ISTAR)) +
                ! THETA*VEC(ISTAR) = (-VEC(ISTAR)*(OLD) COEF(ICOR(ISTAR)) + (OLD)
                ! COEF(ICOR(ISTAR))*VEC(ISTAR))/((OLD) COEF(ICOR(ISTAR)) - VEC(ISTAR))
                ! = 0.0.  NOTE THAT WE HAVE Z2 >= 0.0 AND Z3 > 0.0.
                ! TO CORRECT THIS BLUNDER WE SET MINCF = AN INDEX I FOR WHICH (NEW)
                ! COEF(ICOR(I)) IS MINIMIZED AND SET COEF(ICOR(I)) = 0.0.
                do i = 1, Ncor
                    ii = Icor(i)
                    if (i > 1) then
                        if (Coef(ii) >= amin) cycle
                    end if
                    amin = Coef(ii)
                    mincf = i
                end do
                ii = Icor(mincf)
                Coef(ii) = zero
            end if
            ! INCREMENT THE MINOR ITERATION COUNTER NMIN, REMOVE ICOR(MINCF),
            ! AND DECREMENT NCOR.
            Nmin = Nmin + 1
            do l = 1, Ncor
                if (l > mincf) Icor(l - 1) = Icor(l)
            end do
            Ncor = Ncor - 1
            ! GO BACK TO COMPUTE PTNR.
            goto 200
        else if (Istrt1 + Jflag <= 0) then
            Jflag = 6
            ! write(nwrit,'(A)') '*****JFLAG IS 6 IN CONENR'
            return
        else
            Jflag = -1
            kntsl = 0
            goto 100
        end if

    end subroutine conenr