*> \brief \b SGELQT3 * * Definition: * =========== * * RECURSIVE SUBROUTINE SGELQT3( M, N, A, LDA, T, LDT, INFO ) * * .. Scalar Arguments .. * INTEGER INFO, LDA, M, N, LDT * .. * .. Array Arguments .. * REAL A( LDA, * ), T( LDT, * ) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> DGELQT3 recursively computes a LQ factorization of a real M-by-N *> matrix A, using the compact WY representation of Q. *> *> Based on the algorithm of Elmroth and Gustavson, *> IBM J. Res. Develop. Vol 44 No. 4 July 2000. *> \endverbatim * * Arguments: * ========== * *> \param[in] M *> \verbatim *> M is INTEGER *> The number of rows of the matrix A. M =< N. *> \endverbatim *> *> \param[in] N *> \verbatim *> N is INTEGER *> The number of columns of the matrix A. N >= 0. *> \endverbatim *> *> \param[in,out] A *> \verbatim *> A is REAL array, dimension (LDA,N) *> On entry, the real M-by-N matrix A. On exit, the elements on and *> below the diagonal contain the N-by-N lower triangular matrix L; the *> elements above the diagonal are the rows of V. See below for *> further details. *> \endverbatim *> *> \param[in] LDA *> \verbatim *> LDA is INTEGER *> The leading dimension of the array A. LDA >= max(1,M). *> \endverbatim *> *> \param[out] T *> \verbatim *> T is REAL array, dimension (LDT,N) *> The N-by-N upper triangular factor of the block reflector. *> The elements on and above the diagonal contain the block *> reflector T; the elements below the diagonal are not used. *> See below for further details. *> \endverbatim *> *> \param[in] LDT *> \verbatim *> LDT is INTEGER *> The leading dimension of the array T. LDT >= max(1,N). *> \endverbatim *> *> \param[out] INFO *> \verbatim *> INFO is INTEGER *> = 0: successful exit *> < 0: if INFO = -i, the i-th argument had an illegal value *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \ingroup doubleGEcomputational * *> \par Further Details: * ===================== *> *> \verbatim *> *> The matrix V stores the elementary reflectors H(i) in the i-th row *> above the diagonal. For example, if M=5 and N=3, the matrix V is *> *> V = ( 1 v1 v1 v1 v1 ) *> ( 1 v2 v2 v2 ) *> ( 1 v3 v3 v3 ) *> *> *> where the vi's represent the vectors which define H(i), which are returned *> in the matrix A. The 1's along the diagonal of V are not stored in A. The *> block reflector H is then given by *> *> H = I - V * T * V**T *> *> where V**T is the transpose of V. *> *> For details of the algorithm, see Elmroth and Gustavson (cited above). *> \endverbatim *> * ===================================================================== RECURSIVE SUBROUTINE SGELQT3( M, N, A, LDA, T, LDT, INFO ) * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N, LDT * .. * .. Array Arguments .. REAL A( LDA, * ), T( LDT, * ) * .. * * ===================================================================== * * .. Parameters .. REAL ONE PARAMETER ( ONE = 1.0E+00 ) * .. * .. Local Scalars .. INTEGER I, I1, J, J1, M1, M2, IINFO * .. * .. External Subroutines .. EXTERNAL SLARFG, STRMM, SGEMM, XERBLA * .. * .. Executable Statements .. * INFO = 0 IF( M .LT. 0 ) THEN INFO = -1 ELSE IF( N .LT. M ) THEN INFO = -2 ELSE IF( LDA .LT. MAX( 1, M ) ) THEN INFO = -4 ELSE IF( LDT .LT. MAX( 1, M ) ) THEN INFO = -6 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'SGELQT3', -INFO ) RETURN END IF * IF( M.EQ.1 ) THEN * * Compute Householder transform when N=1 * CALL SLARFG( N, A, A( 1, MIN( 2, N ) ), LDA, T ) * ELSE * * Otherwise, split A into blocks... * M1 = M/2 M2 = M-M1 I1 = MIN( M1+1, M ) J1 = MIN( M+1, N ) * * Compute A(1:M1,1:N) <- (Y1,R1,T1), where Q1 = I - Y1 T1 Y1^H * CALL SGELQT3( M1, N, A, LDA, T, LDT, IINFO ) * * Compute A(J1:M,1:N) = Q1^H A(J1:M,1:N) [workspace: T(1:N1,J1:N)] * DO I=1,M2 DO J=1,M1 T( I+M1, J ) = A( I+M1, J ) END DO END DO CALL STRMM( 'R', 'U', 'T', 'U', M2, M1, ONE, & A, LDA, T( I1, 1 ), LDT ) * CALL SGEMM( 'N', 'T', M2, M1, N-M1, ONE, A( I1, I1 ), LDA, & A( 1, I1 ), LDA, ONE, T( I1, 1 ), LDT) * CALL STRMM( 'R', 'U', 'N', 'N', M2, M1, ONE, & T, LDT, T( I1, 1 ), LDT ) * CALL SGEMM( 'N', 'N', M2, N-M1, M1, -ONE, T( I1, 1 ), LDT, & A( 1, I1 ), LDA, ONE, A( I1, I1 ), LDA ) * CALL STRMM( 'R', 'U', 'N', 'U', M2, M1 , ONE, & A, LDA, T( I1, 1 ), LDT ) * DO I=1,M2 DO J=1,M1 A( I+M1, J ) = A( I+M1, J ) - T( I+M1, J ) T( I+M1, J )=0 END DO END DO * * Compute A(J1:M,J1:N) <- (Y2,R2,T2) where Q2 = I - Y2 T2 Y2^H * CALL SGELQT3( M2, N-M1, A( I1, I1 ), LDA, & T( I1, I1 ), LDT, IINFO ) * * Compute T3 = T(J1:N1,1:N) = -T1 Y1^H Y2 T2 * DO I=1,M2 DO J=1,M1 T( J, I+M1 ) = (A( J, I+M1 )) END DO END DO * CALL STRMM( 'R', 'U', 'T', 'U', M1, M2, ONE, & A( I1, I1 ), LDA, T( 1, I1 ), LDT ) * CALL SGEMM( 'N', 'T', M1, M2, N-M, ONE, A( 1, J1 ), LDA, & A( I1, J1 ), LDA, ONE, T( 1, I1 ), LDT ) * CALL STRMM( 'L', 'U', 'N', 'N', M1, M2, -ONE, T, LDT, & T( 1, I1 ), LDT ) * CALL STRMM( 'R', 'U', 'N', 'N', M1, M2, ONE, & T( I1, I1 ), LDT, T( 1, I1 ), LDT ) * * * * Y = (Y1,Y2); L = [ L1 0 ]; T = [T1 T3] * [ A(1:N1,J1:N) L2 ] [ 0 T2] * END IF * RETURN * * End of SGELQT3 * END