      SUBROUTINE SQRR( M, N, A, LDA, TAU, WORK, LWORK, INFO )
*
*  -- LAPACK routine --
*     Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab,
*     Courant Institute, NAG Ltd., and Rice University
*     December 28, 1989
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, LWORK, M, N
*     ..
*     .. Array Arguments ..
      REAL               A( LDA, * ), TAU( * ), WORK( LWORK )
*     ..
*
*  Purpose
*  =======
*
*  SQRR computes a QR factorization of a real m by n matrix A
*  (m >= n):  A = Q * R.
*
*  This is the blocked right-looking version of the algorithm.
*
*  Arguments
*  =========
*
*  M       (input) INTEGER
*          The number of rows of the matrix A.  M >= 0.
*
*  N       (input) INTEGER
*          The number of columns of the matrix A.  M >= N >= 0.
*
*  A       (input/output) REAL array, dimension (LDA,N)
*          On entry, the m by n matrix A.
*          On exit, the upper triangle of the array contains the n by n
*          upper triangular matrix R; the elements below the diagonal,
*          with the array TAU, represent the orthogonal matrix Q as a
*          product of elementary reflectors.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A.  LDA >= M.
*
*  TAU     (output) REAL array, dimension (N)
*          Further details of the orthogonal matrix Q (see A).
*
*  WORK    (workspace) REAL array, dimension (LWORK)
*
*  LWORK   (input) INTEGER
*          The dimension of the array WORK.  LWORK >= N.
*          The block algorithm will not be used unless LWORK >= N*NB
*          where NB is the block size for this environment.
*
*  INFO    (output) INTEGER
*          = 0: successful exit
*          < 0: if INFO = -i, the i-th argument had an illegal value
*
*  Further Details
*  ===============
*
*  The matrix Q is represented as a product of elementary reflectors
*
*     Q = H(1) H(2) . . . H(n)
*
*  Each H(i) has the form
*
*     H = I - tau * v * v'
*
*  where tau is a real scalar, and v is a real vector with
*  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER            I, IB, NB
*     ..
*     .. External Subroutines ..
      EXTERNAL           ENVIR, SQRR2, SLARFB, SLARFT, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, MIN
*     ..
*     .. Executable Statements ..
*
*     Test the input arguments.
*
      INFO = 0
      IF( M.LT.0 ) THEN
         INFO = -1
      ELSE IF( N.LT.0 .OR. N.GT.M ) THEN
         INFO = -2
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -4
      ELSE IF( LWORK.LT.MAX( 1, N ) ) THEN
         INFO = -7
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SQRR', -INFO )
         RETURN
      END IF
*
*     Quick return if possible.
*
      IF( M.EQ.0 .OR. N.EQ.0 )
     $   RETURN
*
*     Determine the block size.
*
      CALL ENVIR( 'B', NB )
*
      IF( NB.LE.1 .OR. LWORK.LT.N*NB ) THEN
*
*        Use unblocked code.
*
         CALL SQRR2( M, N, A, LDA, TAU, INFO )
      ELSE
*
*        Use blocked code.
*
         DO 10 I = 1, N, NB
            IB = MIN( N-I+1, NB )
*
*           Compute the QR factorization of the current block.
*
            CALL SQRR2( M-I+1, IB, A( I, I ), LDA, TAU( I ), INFO )
            IF( I+IB.LE.N ) THEN
*
*              Form the upper triangular factor of the corresponding
*              block reflector H.
*
               CALL SLARFT( 'Lower', M-I+1, IB, A( I, I ), LDA,
     $                      TAU( I ), WORK, N )
*
*              Apply H' to A(i:m,i+ib:n) from the left.
*
               CALL SLARFB( 'Left', 'Lower', 'Transpose', M-I+1,
     $                      N-I-IB+1, IB, A( I, I ), LDA, WORK, N,
     $                      A( I, I+IB ), LDA, WORK( IB+1 ), N )
            END IF
   10    CONTINUE
      END IF
*
      RETURN
*
*     End of SQRR
*
      END
      SUBROUTINE SQRR2( M, N, A, LDA, TAU, INFO )
*
*  -- LAPACK routine --
*     Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab,
*     Courant Institute, NAG Ltd., and Rice University
*     December 28, 1989
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, M, N
*     ..
*     .. Array Arguments ..
      REAL               A( LDA, * ), TAU( * )
*     ..
*
*  Purpose
*  =======
*
*  SQRR2 computes a QR factorization of a real m by n matrix A
*  (m >= n):  A = Q * R.
*
*  This is the unblocked right-looking version of the algorithm.
*
*  Arguments
*  =========
*
*  M       (input) INTEGER
*          The number of rows of the matrix A.  M >= 0.
*
*  N       (input) INTEGER
*          The number of columns of the matrix A.  M >= N >= 0.
*
*  A       (input/output) REAL array, dimension (LDA,N)
*          On entry, the m by n matrix A.
*          On exit, the upper triangle of the array contains the n by n
*          upper triangular matrix R; the elements below the diagonal,
*          with the array TAU, represent the orthogonal matrix Q as a
*          product of elementary reflectors.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A.  LDA >= M.
*
*  TAU     (output) REAL array, dimension (N)
*          Further details of the orthogonal matrix Q (see A).
*
*  INFO    (output) INTEGER
*          = 0: successful exit
*          < 0: if INFO = -i, the i-th argument had an illegal value
*
*  Further Details
*  ===============
*
*  The matrix Q is represented as a product of elementary reflectors
*
*     Q = H(1) H(2) . . . H(n)
*
*  Each H(i) has the form
*
*     H = I - tau * v * v'
*
*  where tau is a real scalar, and v is a real vector with
*  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ONE
      PARAMETER          ( ONE = 1.0E+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I
      REAL               AII
*     ..
*     .. External Subroutines ..
      EXTERNAL           SLARF, SLARFG, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX
*     ..
*     .. Executable Statements ..
*
*     Test the input arguments.
*
      INFO = 0
      IF( M.LT.0 ) THEN
         INFO = -1
      ELSE IF( N.LT.0 .OR. N.GT.M ) THEN
         INFO = -2
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -4
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SQRR2', -INFO )
         RETURN
      END IF
*
*     Quick return if possible.
*
      IF( M.EQ.0 .OR. N.EQ.0 )
     $   RETURN
*
*     Compute factorization.
*
      DO 10 I = 1, N
*
*        Generate elementary reflector H(i).
*
         IF( I.LT.M ) THEN
            CALL SLARFG( M-I+1, A( I, I ), A( I+1, I ), 1, TAU( I ) )
         ELSE
            CALL SLARFG( 1, A( M, M ), A( M, M ), 1, TAU( M ) )
         END IF
         IF( I.LT.N ) THEN
*
*           Apply H(i) to A(i:m,i+1:n) from the left.
*
            AII = A( I, I )
            A( I, I ) = ONE
            CALL SLARF( 'Left', M-I+1, N-I, A( I, I ), 1, TAU( I ),
     $                  A( I, I+1 ), LDA, TAU( I+1 ) )
            A( I, I ) = AII
         END IF
*
   10 CONTINUE
      RETURN
*
*     End of SQRR2
*
      END
      SUBROUTINE SQRL( M, N, A, LDA, TAU, WORK, LWORK, INFO )
*
*  -- LAPACK routine --
*     Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab,
*     Courant Institute, NAG Ltd., and Rice University
*     January 11, 1990
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, LDS, LS, M, N, NBS
*     ..
*     .. Array Arguments ..
      REAL               A( LDA, * ), TAU( * )
*     ..
*
*  Purpose
*  =======
*
*  SQRL computes a QR factorization of a real m by n matrix A
*  (m >= n):  A = Q * R.
*
*  This is the blocked left-looking version of the algorithm.
*
*  Arguments
*  =========
*
*  M       (input) INTEGER
*          The number of rows of the matrix A.  M >= 0.
*
*  N       (input) INTEGER
*          The number of columns of the matrix A.  M >= N >= 0.
*
*  A       (input/output) REAL array, dimension (LDA,N)
*          On entry, the m by n matrix A.
*          On exit, the upper triangle of the array contains the n by n
*          upper triangular matrix R; the elements below the diagonal,
*          with the array TAU, represent the orthogonal matrix Q as a
*          product of elementary reflectors.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A.  LDA >= M.
*
*  TAU     (output) REAL array, dimension (N)
*          Further details of the orthogonal matrix Q (see A).
*
*  WORK    (workspace) REAL array, dimension (LWORK)
*
*  LWORK   (input) INTEGER
*          The dimension of the array WORK.  LWORK >= N.
*          The block algorithm will not be used unless LWORK >= N*NB
*          where NB is the block size for this environment.
*
*  INFO    (output) INTEGER
*          = 0: successful exit
*          < 0: if INFO = -i, the i-th argument had an illegal value
*
*  Further Details
*  ===============
*
*  The matrix Q is represented as a product of elementary reflectors
*
*     Q = H(1) H(2) . . . H(n)
*
*  Each H(i) has the form
*
*     H = I - tau * v * v'
*
*  where tau is a real scalar, and v is a real vector with
*  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO
      PARAMETER          ( ZERO = 0.0E0 )
*     ..
*     .. Local Scalars ..
      INTEGER            J, K, KB, NB
*     ..
*     .. External Subroutines ..
      EXTERNAL           ENVIR, SLARFB, SLARFT, SQRL2, SQRR2, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, MIN
*     ..
*     .. Executable Statements ..
*
*
*     Test the input parameters.
*
      INFO = 0
      IF( M.LT.0 ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -4
      ELSE IF( LS.LT.N ) THEN
         INFO = -6
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SQRL  ', -INFO )
         RETURN
      END IF
*
*     Quick return if possible.
*
      IF( M.EQ.0 .OR. N.EQ.0 ) THEN
         LDS = 1
         NBS = 1
         TAU( 1 ) = ZERO
         RETURN
      END IF
*
*     Determine the correct block size for this environment.
*
      CALL ENVIR( 'Block', NB )
      NB = MIN( NB, LS/N )
*
      LDS = N
      NBS = NB
      IF( NB.EQ.1 ) THEN
         CALL SQRL2( M, N, A, LDA, TAU, INFO )
      ELSE
*
         DO 20 K = 1, N, NB
            KB = MIN( N-K+1, NB )
*
*           Apply the previous transformations to the current block.
*
            DO 10 J = 1, K - NB, NB
               CALL SLARFB( 'Left', 'Transpose', M-J+1, KB, NB,
     $                    A( J, J ), LDA, TAU( J ), LDS, A( J, K ), LDA,
     $                      TAU( K ), LDS )
   10       CONTINUE
*
*           Find the QR factorization of the current block.
*
            CALL SQRR2( M-K+1, KB, A( K, K ), LDA, TAU( K ), INFO )
*
*           Accumulate the block Householder matrix.
*
            CALL SLARFT( M-K+1, KB, A( K, K ), LDA, TAU( K ), LDS )
   20    CONTINUE
      END IF
      RETURN
*
*     End of SQRL
*
      END
      SUBROUTINE SQRL2( M, N, A, LDA, S, INFO )
*
*  -- LAPACK routine --
*     Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab,
*     Courant Institute, NAG Ltd., and Rice University
*     November 29, 1989
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, M, N
*     ..
*     .. Array Arguments ..
      REAL               A( LDA, * ), S( * )
*     ..
*
*  Purpose
*  =======
*
*     Computes a QR factorization of an m-by-n matrix A (m.ge.n).
*     This version assumes that elementary Householder matrices are
*     represented in the form  I - alpha*v*v'  where v(1) = 1.
*
*     This is the left-looking unblocked version of the algorithm.
*
*  Arguments
*  =========
*
*  M      - INTEGER.
*           On entry, M specifies the number of rows of the matrix
*           A .  M must be at least zero.
*           Unchanged on exit.
*
*  N      - INTEGER.
*           On entry, N specifies the number of columns of the matrix
*           A . N must be at least zero.
*           Unchanged on exit.
*
*  A      - REAL             array of DIMENSION ( LDA, N ).
*           On entry, A specifies the array which contains the matrix
*           being factored.
*           On exit, the array A is overwritten by the
*           QR factorization. The factorization can be written as
*           A = Q*R where Q is an orthogonal matrix and R is an upper
*           triangular matrix.
*
*  LDA    - INTEGER.
*           On entry, LDA specifies the first dimension of A as declared
*           in the calling (sub) program.
*           LDA must be at least  max( 1, M ).
*           Unchanged on exit.
*
*  S      - REAL             array of DIMENSION ( N ).
*           S is used as a temporary array.
*
*  INFO   - INTEGER.
*           On exit, a value of 0 indicates a normal return;
*           a negative value, say -K, indicates the Kth argument has an
*           illegal value.
*
*     .. Parameters ..
      REAL               ONE, ZERO
      PARAMETER          ( ONE = 1.0E0, ZERO = 0.0E0 )
*     ..
*     .. Local Scalars ..
      INTEGER            J, K
      REAL               ALPHA, DJ, T
*     ..
*     .. External Functions ..
      REAL               SDOT
      EXTERNAL           SDOT
*     ..
*     .. External Subroutines ..
      EXTERNAL           SAXPY, SLARFG
*     ..
*     .. Executable Statements ..
*
*
      DO 20 K = 1, N
*
*        Apply the previous transformations to the current column.
*
         DO 10 J = 1, K - 1
            IF( S( J ).NE.ZERO ) THEN
               DJ = A( J, J )
               A( J, J ) = ONE
               ALPHA = -S( J )*SDOT( M-J+1, A( J, J ), 1, A( J, K ), 1 )
               CALL SAXPY( M-J+1, ALPHA, A( J, J ), 1, A( J, K ), 1 )
               A( J, J ) = DJ
            END IF
   10    CONTINUE
*
*        Generate reflection
*
         T = ZERO
         IF( K+1.LE.M ) THEN
            CALL SLARFG( M-K+1, A( K, K ), A( K+1, K ), 1, T )
         END IF
         S( K ) = T
   20 CONTINUE
*
      RETURN
*
*     End of SQRL2
*
      END