house Subroutine

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

        implicit none

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

        real(wp) :: aakk, amax, sqdk, store, sum, summ, test, testt
        integer :: i, ia, icount, ii, j, jj, k, kchnge, kk, &
                   kp, krank, kt

        real(wp), parameter :: tolsq = (ten*ten*spcmn)**2
        integer, parameter :: numref = 1 !! Set numref = the limit on the number of iterative refinement steps.
        integer, parameter :: nmref1 = numref + 1
        integer, parameter :: nmref2 = numref + 2

        ! COMPUTE MACHINE AND PRECISION DEPENDENT CONSTANTS.
        Ihouse = 0

        ! SET KRANK = MIN(N,NCOR).  THIS MAY BE REDUCED LATER.
        krank = Ncor
        if (n < Ncor) krank = n
        ! INITIALLY SET KPIVOT.  AFTER ALL COLUMN INTERCHANGES ARE DONE
        ! KPIVOT(J) WILL BE THE ORIGINAL POSITION OF THE COLUMN WHERE THE
        ! JTH PIVOT WAS DONE.  THIS COLUMN WILL BE MOVED TO COLUMN J.
        do j = 1, Ncor
            Kpivot(j) = j
        end do
        ! COPY R INTO B AND PICOR INTO AA, BUT IN THE PROCESS REPLACE ANY NUMBERS
        ! WITH ABSOLUTE VALUE LESS THAN SPCMN BY ZERO TO AVOID UNDERFLOWS.
        do i = 1, n
            if (abs(r(i)) < spcmn) then
                b(i) = zero
            else
                b(i) = r(i)
            end if
        end do
        do j = 1, Ncor
            do i = 1, n
                if (abs(Picor(i, j)) < spcmn) then
                    Aa(i, j) = zero
                else
                    Aa(i, j) = Picor(i, j)
                end if
            end do
        end do
        do k = 1, Ncor
            if (k > n) exit
            d(k) = zero
            kchnge = k
            do jj = k, Ncor
                sum = zero
                do ia = k, n
                    if (abs(Aa(ia, jj)) > spcmn) sum = sum + Aa(ia, jj)*Aa(ia, jj)
                end do
                if (d(k) < sum) then
                    kchnge = jj
                    d(k) = sum
                end if
            end do

            !  KCHNGE CONTAINS THE INDEX OF THE COLUMN OF GREATEST
            !  LENGTH BETWEEN K AND NCOR (FROM POSITION K TO THE BOTTOM).
            ! IF K=1 AND D(K) < TOLSQ WE RETURN WITH THE FAILURE WARNING
            ! IHOUSE=1.
            if (k <= 1) then
                if (d(k) < tolsq) then
                    Ihouse = 1
                    return
                end if
            end if

            if (kchnge /= k) then
                !  START COLUMN INTERCHANGE.
                do i = 1, n
                    store = Aa(i, kchnge)
                    Aa(i, kchnge) = Aa(i, k)
                    Aa(i, k) = store
                end do
                kk = Kpivot(k)
                Kpivot(k) = Kpivot(kchnge)
                Kpivot(kchnge) = kk
            end if
            if (k /= 1) then
                amax = abs(d(1))
                test = (real(n - k + 1, wp)*(ten*ten*spcmn)**2)*(amax*amax)
                if (abs(d(k)) <= test) then
                    ! HERE THE LENGTH OF THE BEST OF COLUMNS K THROUGH NCOR (FROM K DOWN)
                    ! WAS TOO SMALL, AND WE REDUCE KRANK TO K-1 AND LEAVE THIS LOOP.
                    d(k) = sqrt(d(k))
                    krank = k - 1
                    exit
                end if
            end if

            ! NOW COMPUTE THE SCALAR BETA(K) AND THE N-K+1 DIMENSIONAL VECTOR
            ! GNU(K) (TO BE PLACED IN AA(K,K),...,AA(N,K)) FOR I(K) - BETA(K)*
            ! GNU(K)*(GNU(K) TRANSPOSE), WHICH IS THE ACTIVE PART OF THE
            ! HOUSEHOLDER TRANSFORMATION PH(K) = DIAG(I(K-1), ACTIVE PART).  THIS
            ! IS A SYMMETRIC ORTHOGONAL MATRIX WHICH WHEN MULTIPLIED TIMES AA WILL
            ! ZERO OUT AA(K+1,K),...,AA(N,K) AND CHANGE AA(K,K) TO -SGN(OLD
            ! AA(K,K))*SQDK, WHERE SQDK = LENGTH OF OLD (AA(K,K),...,AA(N,K)) AND
            ! WE REDEFINE THE SGN FUNCTION TO HAVE VALUE 1.0 IF ITS ARGUMENT IS
            ! 0.0.  WE WILL HAVE BETA(K) = 1.0/(SQDK**2 + ABS(OLD AA(K,K))*SQDK)
            ! AND GNU(K) = (OLD AA(K,K) + SGN(OLD AA(K,K))*SQDK, OLD AA(K+1,K),...,
            ! OLD AA(N,K)).  WE WILL ALSO REPLACE D(K) BY THE NEW AA(K,K) (WHICH
            ! WILL NOT ACTUALLY BE WRITTEN INTO AA) FOR LATER USE.
            aakk = Aa(k, k)
            sqdk = sqrt(d(k))
            if (aakk < zero) then
                Beta(k) = one/(d(k) - aakk*sqdk)
                Aa(k, k) = -sqdk + aakk
                d(k) = sqdk
            else
                Beta(k) = one/(d(k) + aakk*sqdk)
                Aa(k, k) = sqdk + aakk
                d(k) = -sqdk
            end if
            kt = k + 1
            if (k /= Ncor) then
                ! HERE K < NCOR AND WE MULTIPLY COLUMNS K+1,...,NCOR OF AA BY THE
                ! HOUSEHOLDER TRANSFORMATION PH(K), WHICH WILL CHANGE ONLY POSITIONS
                ! K THROUGH THE BOTTOM OF THESE COLUMNS.  THIS IS DONE BY, FOR J =
                ! K+1,...,NCOR, REPLACING COLUMN J (FROM K DOWN) BY COLUMN J (FROM K DOWN)
                ! - GNU(K)*(GNU(K).COLUMN J (FROM K DOWN))*BETA(K).
                do j = kt, Ncor
                    Save(j) = zero
                    do ia = k, n
                        Save(j) = Save(j) + Aa(ia, k)*Aa(ia, j)
                    end do
                    do i = k, n
                        Aa(i, j) = Aa(i, j) - Aa(i, k)*Save(j)*Beta(k)
                    end do
                end do
            end if
        end do

        do i = 1, krank
            ! IF I <= MIN(KRANK,NCOR-1), DIVIDE ROW I OF AA FROM COLUMN I+1
            ! THROUGH COLUMN NCOR BY THE NEW AA(I,I) (WHICH IS NOT ACTUALLY
            ! WRITTEN INTO AA(I,I), BUT IS STORED IN D(I)).
            ii = i + 1
            if (i == Ncor) exit
            do j = ii, Ncor
                Aa(i, j) = Aa(i, j)/d(i)
            end do
        end do

        ! NOW ALL THE DIAGONAL ELEMENTS OF AA (ALTHOUGH NOT WRITTEN IN)
        ! ARE 1.0 AND ALL OFF DIAGONAL ELEMENTS OF AA ARE LESS THAN OR
        ! EQUAL TO 1.0.

        ! INITIALIZE THE ITERATIVE REFINEMENT COUNTER ICOUNT AND ZERO OUT VEC
        ! INITIALLY.  THE VEC VALUES NOT CORRESPONDING TO THE FIRST KRANK
        ! COLUMNS (MODULO EARLIER COLUMN INTERCHANGES) WILL REMAIN AT 0.0.
        icount = 1
        do i = 1, Ncor
            Vec(i) = zero
        end do

        iteration: do

            ! PREMULTIPLY B BY THE HOUSEHOLDER TRANSFORMATIONS PH(1),...,
            ! PH(KRANK).  RECALL THAT GNU(I) IS STILL IN AA(I,I),...,AA(N,I)
            ! FOR I=1,...,KRANK.
            do i = 1, krank
                sum = zero
                do ia = i, n
                    sum = sum + Aa(ia, i)*b(ia)
                end do
                sum = sum*Beta(i)
                do j = i, n
                    b(j) = b(j) - Aa(j, i)*sum
                end do
            end do

            ! NOW ONLY USE THE FIRST KRANK TERMS OF B, AS WE CANT DO ANYTHING ABOUT
            ! THE OTHERS, WHOSE SQUARE ROOT OF SUM OF SQUARES WILL GIVE THE LEAST
            ! SQUARES DISTANCE.
            ! DIVIDE B(I) BY D(I) FOR I=1,...,KRANK AS WE DID THIS TO ROW I OF AA.
            do i = 1, krank
                b(i) = b(i)/d(i)
            end do

            ! THE PROBLEM HAS NOW BEEN REDUCED TO SOLVING (UPPER LEFT KRANK BY
            ! KRANK PART OF AA)*(FIRST KRANK TERMS OF VEC, MODULO COLUMN
            ! INTERCHANGE UNSCRAMBLING) = (FIRST KRANK TERMS OF B).  ALTHOUGH THE
            ! DIAGONAL AND BELOW DIAGONAL TERMS OF THE COEFFICIENT MATRIX HAVE NOT
            ! BEEN WRITTEN IN, THE SYSTEM IS UPPER TRIANGULAR WITH DIAGONAL ELEMENTS
            ! ALL EQUAL TO 1.0, SO WE SOLVE BY BACK SUBSTITUTION.  WE FIRST PUT
            ! THE SOLUTION TO THIS SYSTEM IN B(1),...,B(KRANK) AND SORT IT OUT
            ! LATER.  IF ICOUNT > 1 THE SOLUTION IS AN ITERATIVE CORRECTION TO
            ! VEC RATHER THAN VEC ITSELF.
            do ii = 1, krank
                i = krank + 1 - ii
                kk = i - 1
                if (i /= 1) then
                    ! HERE WE ALREADY HAVE B(I) (WHERE I  > 1) AND WE SUBTRACT AA(J,I)*
                    ! B(I) FROM B(J) FOR J = 1,...,I-1.
                    do j = 1, kk
                        b(j) = b(j) - Aa(j, i)*b(i)
                    end do
                end if
            end do

            !  TEST FOR CONVERGENCE.
            !  FIRST TEST, TOO MANY ITERATIONS.
            !  SECOND TEST, SEE IF VEC IS DECREASING.

            ! COMPUTE THE LENGTH SQUARED OF THE FIRST TOP 1 THROUGH KRANK PART OF
            ! B, WHICH WILL BE THE RESIDUAL VECTOR IF ICOUNT > 1.
            sum = zero
            do i = 1, krank
                if (abs(b(i)) > spcmn) sum = sum + b(i)*b(i)
            end do
            if (icount == 1) then
                testt = sum
            else if (sum > test/two) then
                icount = nmref2
            end if
            test = sum

            ! COMPUTE THE VEC VALUES, WHICH WILL BE ACTUAL VEC VALUES IF ICOUNT=1
            ! AND CORRECTIONS TO VEC VALUES IF ICOUNT > 1.  WE GET THESE BY
            ! UNSCRAMBLING THE B VALUES AND ADDING THEM TO THE APPROPRIATE OLD VEC
            ! VALUES (WHICH WILL BE 0.0 IF ICOUNT=1).
            do i = 1, krank
                kp = Kpivot(i)
                Vec(kp) = b(i) + Vec(kp)
            end do

            ! CALCULATE THE RESIDUAL R - ACOEF*VEC.  RECALL THAT ACOEF AND R
            ! CONTAIN THE ORIGINAL COEFFICIENT AND RIGHT SIDE ARRAYS RESPECTIVELY.
            ! TO GET THE MOST OUT OF ITERATIVE REFINEMENT THIS COMPUTATION SHOULD
            ! PROBABLY BE DONE IN HIGHER PRECISION, IN WHICH CASE IT MAY BE
            ! FRUITFUL TO ALSO SET NUMREF LARGER AT THE BEGINNING OF THIS
            ! SUBROUTINE.
            do i = 1, n
                summ = zero
                do j = 1, Ncor
                    if (abs(Picor(i, j)) >= spcmn) summ = summ + Picor(i, j)*Vec(j)
                end do
                b(i) = r(i) - summ
            end do

            !  THIRD TEST, WAS THE CORRECTION SIGNIFICANT.
            if (test >= spcmn*testt) then
                if (icount /= nmref1) then
                    if (icount < nmref2) then
                        icount = icount + 1
                        cycle iteration
                    end if
                end if
            end if

            exit iteration  ! done

        end do iteration

    end subroutine house