SUBROUTINE SORGQR( M, N, K, A, LDA, TAU, WORK, LWORK, INFO ) * * -- LAPACK routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER INFO, K, LDA, LWORK, M, N * .. * .. Array Arguments .. REAL A( LDA, * ), TAU( * ), WORK( * ) * .. * * Purpose * ======= * * SORGQR generates an M-by-N real matrix Q with orthonormal columns, * which is defined as the first N columns of a product of K elementary * reflectors of order M * * Q = H(1) H(2) . . . H(k) * * as returned by SGEQRF. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix Q. M >= 0. * * N (input) INTEGER * The number of columns of the matrix Q. M >= N >= 0. * * K (input) INTEGER * The number of elementary reflectors whose product defines the * matrix Q. N >= K >= 0. * * A (input/output) REAL array, dimension (LDA,N) * On entry, the i-th column must contain the vector which * defines the elementary reflector H(i), for i = 1,2,...,k, as * returned by SGEQRF in the first k columns of its array * argument A. * On exit, the M-by-N matrix Q. * * LDA (input) INTEGER * The first dimension of the array A. LDA >= max(1,M). * * TAU (input) REAL array, dimension (K) * TAU(i) must contain the scalar factor of the elementary * reflector H(i), as returned by SGEQRF. * * WORK (workspace/output) REAL array, dimension (MAX(1,LWORK)) * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. * * LWORK (input) INTEGER * The dimension of the array WORK. LWORK >= max(1,N). * For optimum performance LWORK >= N*NB, where NB is the * optimal blocksize. * * If LWORK = -1, then a workspace query is assumed; the routine * only calculates the optimal size of the WORK array, returns * this value as the first entry of the WORK array, and no error * message related to LWORK is issued by XERBLA. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument has an illegal value * * ===================================================================== * * .. Parameters .. REAL ZERO PARAMETER ( ZERO = 0.0E+0 ) * .. * .. Local Scalars .. LOGICAL LQUERY INTEGER I, IB, IINFO, IWS, J, KI, KK, L, LDWORK, \$ LWKOPT, NB, NBMIN, NX * .. * .. External Subroutines .. EXTERNAL SLARFB, SLARFT, SORG2R, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. External Functions .. INTEGER ILAENV EXTERNAL ILAENV * .. * .. Executable Statements .. * * Test the input arguments * INFO = 0 NB = ILAENV( 1, 'SORGQR', ' ', M, N, K, -1 ) LWKOPT = MAX( 1, N )*NB WORK( 1 ) = LWKOPT LQUERY = ( LWORK.EQ.-1 ) IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 .OR. N.GT.M ) THEN INFO = -2 ELSE IF( K.LT.0 .OR. K.GT.N ) THEN INFO = -3 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -5 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN INFO = -8 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'SORGQR', -INFO ) RETURN ELSE IF( LQUERY ) THEN RETURN END IF * * Quick return if possible * IF( N.LE.0 ) THEN WORK( 1 ) = 1 RETURN END IF * NBMIN = 2 NX = 0 IWS = N IF( NB.GT.1 .AND. NB.LT.K ) THEN * * Determine when to cross over from blocked to unblocked code. * NX = MAX( 0, ILAENV( 3, 'SORGQR', ' ', M, N, K, -1 ) ) IF( NX.LT.K ) THEN * * Determine if workspace is large enough for blocked code. * LDWORK = N IWS = LDWORK*NB IF( LWORK.LT.IWS ) THEN * * Not enough workspace to use optimal NB: reduce NB and * determine the minimum value of NB. * NB = LWORK / LDWORK NBMIN = MAX( 2, ILAENV( 2, 'SORGQR', ' ', M, N, K, -1 ) ) END IF END IF END IF * IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN * * Use blocked code after the last block. * The first kk columns are handled by the block method. * KI = ( ( K-NX-1 ) / NB )*NB KK = MIN( K, KI+NB ) * * Set A(1:kk,kk+1:n) to zero. * DO 20 J = KK + 1, N DO 10 I = 1, KK A( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE KK = 0 END IF * * Use unblocked code for the last or only block. * IF( KK.LT.N ) \$ CALL SORG2R( M-KK, N-KK, K-KK, A( KK+1, KK+1 ), LDA, \$ TAU( KK+1 ), WORK, IINFO ) * IF( KK.GT.0 ) THEN * * Use blocked code * DO 50 I = KI + 1, 1, -NB IB = MIN( NB, K-I+1 ) IF( I+IB.LE.N ) THEN * * Form the triangular factor of the block reflector * H = H(i) H(i+1) . . . H(i+ib-1) * CALL SLARFT( 'Forward', 'Columnwise', M-I+1, IB, \$ A( I, I ), LDA, TAU( I ), WORK, LDWORK ) * * Apply H to A(i:m,i+ib:n) from the left * CALL SLARFB( 'Left', 'No transpose', 'Forward', \$ 'Columnwise', M-I+1, N-I-IB+1, IB, \$ A( I, I ), LDA, WORK, LDWORK, A( I, I+IB ), \$ LDA, WORK( IB+1 ), LDWORK ) END IF * * Apply H to rows i:m of current block * CALL SORG2R( M-I+1, IB, IB, A( I, I ), LDA, TAU( I ), WORK, \$ IINFO ) * * Set rows 1:i-1 of current block to zero * DO 40 J = I, I + IB - 1 DO 30 L = 1, I - 1 A( L, J ) = ZERO 30 CONTINUE 40 CONTINUE 50 CONTINUE END IF * WORK( 1 ) = IWS RETURN * * End of SORGQR * END