00001 SUBROUTINE SGEQRFP( M, N, A, LDA, TAU, WORK, LWORK, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.2.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * June 2010 00007 * 00008 * .. Scalar Arguments .. 00009 INTEGER INFO, LDA, LWORK, M, N 00010 * .. 00011 * .. Array Arguments .. 00012 REAL A( LDA, * ), TAU( * ), WORK( * ) 00013 * .. 00014 * 00015 * Purpose 00016 * ======= 00017 * 00018 * SGEQRFP computes a QR factorization of a real M-by-N matrix A: 00019 * A = Q * R. 00020 * 00021 * Arguments 00022 * ========= 00023 * 00024 * M (input) INTEGER 00025 * The number of rows of the matrix A. M >= 0. 00026 * 00027 * N (input) INTEGER 00028 * The number of columns of the matrix A. N >= 0. 00029 * 00030 * A (input/output) REAL array, dimension (LDA,N) 00031 * On entry, the M-by-N matrix A. 00032 * On exit, the elements on and above the diagonal of the array 00033 * contain the min(M,N)-by-N upper trapezoidal matrix R (R is 00034 * upper triangular if m >= n); the elements below the diagonal, 00035 * with the array TAU, represent the orthogonal matrix Q as a 00036 * product of min(m,n) elementary reflectors (see Further 00037 * Details). 00038 * 00039 * LDA (input) INTEGER 00040 * The leading dimension of the array A. LDA >= max(1,M). 00041 * 00042 * TAU (output) REAL array, dimension (min(M,N)) 00043 * The scalar factors of the elementary reflectors (see Further 00044 * Details). 00045 * 00046 * WORK (workspace/output) REAL array, dimension (MAX(1,LWORK)) 00047 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00048 * 00049 * LWORK (input) INTEGER 00050 * The dimension of the array WORK. LWORK >= max(1,N). 00051 * For optimum performance LWORK >= N*NB, where NB is 00052 * the optimal blocksize. 00053 * 00054 * If LWORK = -1, then a workspace query is assumed; the routine 00055 * only calculates the optimal size of the WORK array, returns 00056 * this value as the first entry of the WORK array, and no error 00057 * message related to LWORK is issued by XERBLA. 00058 * 00059 * INFO (output) INTEGER 00060 * = 0: successful exit 00061 * < 0: if INFO = -i, the i-th argument had an illegal value 00062 * 00063 * Further Details 00064 * =============== 00065 * 00066 * The matrix Q is represented as a product of elementary reflectors 00067 * 00068 * Q = H(1) H(2) . . . H(k), where k = min(m,n). 00069 * 00070 * Each H(i) has the form 00071 * 00072 * H(i) = I - tau * v * v' 00073 * 00074 * where tau is a real scalar, and v is a real vector with 00075 * v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), 00076 * and tau in TAU(i). 00077 * 00078 * ===================================================================== 00079 * 00080 * .. Local Scalars .. 00081 LOGICAL LQUERY 00082 INTEGER I, IB, IINFO, IWS, K, LDWORK, LWKOPT, NB, 00083 $ NBMIN, NX 00084 * .. 00085 * .. External Subroutines .. 00086 EXTERNAL SGEQR2P, SLARFB, SLARFT, XERBLA 00087 * .. 00088 * .. Intrinsic Functions .. 00089 INTRINSIC MAX, MIN 00090 * .. 00091 * .. External Functions .. 00092 INTEGER ILAENV 00093 EXTERNAL ILAENV 00094 * .. 00095 * .. Executable Statements .. 00096 * 00097 * Test the input arguments 00098 * 00099 INFO = 0 00100 NB = ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 ) 00101 LWKOPT = N*NB 00102 WORK( 1 ) = LWKOPT 00103 LQUERY = ( LWORK.EQ.-1 ) 00104 IF( M.LT.0 ) THEN 00105 INFO = -1 00106 ELSE IF( N.LT.0 ) THEN 00107 INFO = -2 00108 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00109 INFO = -4 00110 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN 00111 INFO = -7 00112 END IF 00113 IF( INFO.NE.0 ) THEN 00114 CALL XERBLA( 'SGEQRFP', -INFO ) 00115 RETURN 00116 ELSE IF( LQUERY ) THEN 00117 RETURN 00118 END IF 00119 * 00120 * Quick return if possible 00121 * 00122 K = MIN( M, N ) 00123 IF( K.EQ.0 ) THEN 00124 WORK( 1 ) = 1 00125 RETURN 00126 END IF 00127 * 00128 NBMIN = 2 00129 NX = 0 00130 IWS = N 00131 IF( NB.GT.1 .AND. NB.LT.K ) THEN 00132 * 00133 * Determine when to cross over from blocked to unblocked code. 00134 * 00135 NX = MAX( 0, ILAENV( 3, 'SGEQRF', ' ', M, N, -1, -1 ) ) 00136 IF( NX.LT.K ) THEN 00137 * 00138 * Determine if workspace is large enough for blocked code. 00139 * 00140 LDWORK = N 00141 IWS = LDWORK*NB 00142 IF( LWORK.LT.IWS ) THEN 00143 * 00144 * Not enough workspace to use optimal NB: reduce NB and 00145 * determine the minimum value of NB. 00146 * 00147 NB = LWORK / LDWORK 00148 NBMIN = MAX( 2, ILAENV( 2, 'SGEQRF', ' ', M, N, -1, 00149 $ -1 ) ) 00150 END IF 00151 END IF 00152 END IF 00153 * 00154 IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN 00155 * 00156 * Use blocked code initially 00157 * 00158 DO 10 I = 1, K - NX, NB 00159 IB = MIN( K-I+1, NB ) 00160 * 00161 * Compute the QR factorization of the current block 00162 * A(i:m,i:i+ib-1) 00163 * 00164 CALL SGEQR2P( M-I+1, IB, A( I, I ), LDA, TAU( I ), WORK, 00165 $ IINFO ) 00166 IF( I+IB.LE.N ) THEN 00167 * 00168 * Form the triangular factor of the block reflector 00169 * H = H(i) H(i+1) . . . H(i+ib-1) 00170 * 00171 CALL SLARFT( 'Forward', 'Columnwise', M-I+1, IB, 00172 $ A( I, I ), LDA, TAU( I ), WORK, LDWORK ) 00173 * 00174 * Apply H' to A(i:m,i+ib:n) from the left 00175 * 00176 CALL SLARFB( 'Left', 'Transpose', 'Forward', 00177 $ 'Columnwise', M-I+1, N-I-IB+1, IB, 00178 $ A( I, I ), LDA, WORK, LDWORK, A( I, I+IB ), 00179 $ LDA, WORK( IB+1 ), LDWORK ) 00180 END IF 00181 10 CONTINUE 00182 ELSE 00183 I = 1 00184 END IF 00185 * 00186 * Use unblocked code to factor the last or only block. 00187 * 00188 IF( I.LE.K ) 00189 $ CALL SGEQR2P( M-I+1, N-I+1, A( I, I ), LDA, TAU( I ), WORK, 00190 $ IINFO ) 00191 * 00192 WORK( 1 ) = IWS 00193 RETURN 00194 * 00195 * End of SGEQRFP 00196 * 00197 END
1.7.2