      SUBROUTINE ORTHO(M, N, QR, MS, S, IP, EXIT)
C***********************************************************************
C      IDENTIFICATION
C        ORTHO - ORTHOGONAL TRANSFORMATION OF A GIVEN GENERAL M BY N
C                MATRIX A TO UPPER TRIANGULAR FORM
C        FORTRAN SUBROUTINE SUBPROGRAM
C        AEROSPACE RESEARCH LABORATORIES
C        WRIGHT-PATTERSON AFB, OHIO  45433
C      PURPOSE
C        ORTHO COMPUTES AN IMPLICIT ORTHOGONAL MATRIX Q AND AN EXPLICIT
C        UPPER TRIANGULAR MATRIX R AND A PERMUTATION MATRIX P SATISFYING
C        QR = PA GIVEN AN M BY N REAL MATRIX A.  ORTHO IS INTENDED FOR
C        USE WITH THE SUBROUTINE ORSOL TO PRODUCE THE LEAST SQUARES
C        SOLUTION OF THE EQUATION AX = B.
C      CONTROL
C
C        DIMENSION QR(MS,N), S(N), IP(N)
C        LOGICAL EXIT
C        .
C        .
C        .
C        CALL ORTHO(M, N, QR, MS, S, IP, EXIT)
C
C        WHERE
C        M    IS AN INTEGER INPUT VARIABLE, THE NUMBER OF ROWS OF A.
C        N    IS AN INTEGER INPUT VARIABLE, THE NUMBER OF COLUMNS OF A,
C             (1 .LT. N .LE. M).
C        QR   AS A REAL INPUT ARRAY IS MATRIX A TO BE TRIANGULARIZED.
C        QR   AS A REAL OUTPUT ARRAY IS THE UPPER TRIANGULAR FACTOR R IN
C             QR(I,J), I .LE. J, AND THE RELEVANT PARTS OF Q IN QR(I,J),
C             I .GT. J.
C        MS   IS AN INTEGER INPUT VARIABLE, THE LEADING DIMENSION OF
C             QR IN THE CALLING PROGRAM.
C        S    IS A REAL OUTPUT ARRAY, THE RELEVANT PARTS OF Q.
C        IP   IS AN INTEGER OUTPUT ARRAY CONTAINING IN IP(I), I=1,...,N,
C             THE IMAGES OF THE PERMUTATION CORRESPONDING TO THE PERMU-
C             TATION MATRIX P.
C        EXIT IS SET TO THE VALUE .TRUE. IF THE RANK OF A IS EQUAL TO N
C             AND .FALSE. OTHERWISE.
C      METHOD
C        THE MATRIX A IN THE ARRAY QR IS REDUCED TO UPPER TRIANGULAR
C        FORM USING ORTHOGONAL TRANSFORMATION WITH PARTIAL PIVOTING.
C      REFERENCES
C        (1) PETER BUSINGER AND G.H. GOLUB, LINEAR LEAST SQUARES SOLU-
C            TIONS BY HOUSEHOLDER TRANSFORMATIONS, NUMER. MATH. 7(1965),
C            269-276.
C        (2) N.K. TSAO AND P.J. NIKOLAI, PROCEDURES USING ORTHOGONAL
C            TRANSFORMATIONS FOR LINEAR LEAST SQUARES PROBLEMS, ARL
C            TECHNICAL REPORT ARL TR 74-0124(1974).
C***********************************************************************
      DIMENSION QR(MS,N), S(N), IP(N)
      LOGICAL EXIT
      EXIT = .TRUE.
      NN = N
      IF(N .EQ. M) NN = N - 1
      DO 8 J = 1, NN
        IP(J) = J
        JP = J + 1
        KJ = J
C                                       SEARCH FOR PIVOT IN THE J-TH
C                                       COLUMN AND INTERCHANGE ROWS.
        DO 1 K = JP, M
    1     IF(ABS(QR(K,J)) .GT. ABS(QR(KJ,J))) KJ = K
        IF(QR(KJ,J) .EQ. 0.) GO TO 9
        IF(KJ .EQ. J) GO TO 30
        IP(J) = KJ
        DO 2 I = J, N
          SAV = QR(J,I)
          QR(J,I) = QR(KJ,I)
    2     QR(KJ,I) = SAV
C                                       NORMALIZE THE PIVOTING COLUMN
C                                       AND FIND ITS NORM.
   30   AJJ = QR(J,J)
        DO 3 I = JP, M
    3     QR(I,J) = QR(I,J)/AJJ
        SAV = 1.
        DO 4 I = JP, M
    4     SAV = SAV + QR(I,J)*QR(I,J)
        S(J) = -SQRT(SAV)
        QR(J,J) = S(J)*AJJ
        IF(JP .GT. N) GO TO 8
C                                       PREMULTIPLY QR WITH THE J-TH
C                                       ORTHOGONAL MATRIX.
        Y = 1. - S(J)
        DO 7 K = JP, N
          SAV = QR(J,K)
          DO 5 I = JP, M
    5       SAV = SAV + QR(I,J)*QR(I,K)
          SS = QR(J,K)
          QR(J,K) = SAV/S(J)
          SS = (SS - QR(J,K))/Y
          DO 6 I = JP, M
    6       QR(I,K) = QR(I,K) - QR(I,J)*SS
    7   CONTINUE
    8 CONTINUE
      RETURN
    9 EXIT = .FALSE.
      RETURN
      END
