![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DGEQRFP 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DGEQRFP + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeqrfp.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeqrfp.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeqrfp.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DGEQRFP( M, N, A, LDA, TAU, WORK, LWORK, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * INTEGER INFO, LDA, LWORK, M, N 00025 * .. 00026 * .. Array Arguments .. 00027 * DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) 00028 * .. 00029 * 00030 * 00031 *> \par Purpose: 00032 * ============= 00033 *> 00034 *> \verbatim 00035 *> 00036 *> DGEQRFP computes a QR factorization of a real M-by-N matrix A: 00037 *> A = Q * R. 00038 *> \endverbatim 00039 * 00040 * Arguments: 00041 * ========== 00042 * 00043 *> \param[in] M 00044 *> \verbatim 00045 *> M is INTEGER 00046 *> The number of rows of the matrix A. M >= 0. 00047 *> \endverbatim 00048 *> 00049 *> \param[in] N 00050 *> \verbatim 00051 *> N is INTEGER 00052 *> The number of columns of the matrix A. N >= 0. 00053 *> \endverbatim 00054 *> 00055 *> \param[in,out] A 00056 *> \verbatim 00057 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00058 *> On entry, the M-by-N matrix A. 00059 *> On exit, the elements on and above the diagonal of the array 00060 *> contain the min(M,N)-by-N upper trapezoidal matrix R (R is 00061 *> upper triangular if m >= n); the elements below the diagonal, 00062 *> with the array TAU, represent the orthogonal matrix Q as a 00063 *> product of min(m,n) elementary reflectors (see Further 00064 *> Details). 00065 *> \endverbatim 00066 *> 00067 *> \param[in] LDA 00068 *> \verbatim 00069 *> LDA is INTEGER 00070 *> The leading dimension of the array A. LDA >= max(1,M). 00071 *> \endverbatim 00072 *> 00073 *> \param[out] TAU 00074 *> \verbatim 00075 *> TAU is DOUBLE PRECISION array, dimension (min(M,N)) 00076 *> The scalar factors of the elementary reflectors (see Further 00077 *> Details). 00078 *> \endverbatim 00079 *> 00080 *> \param[out] WORK 00081 *> \verbatim 00082 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00083 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00084 *> \endverbatim 00085 *> 00086 *> \param[in] LWORK 00087 *> \verbatim 00088 *> LWORK is INTEGER 00089 *> The dimension of the array WORK. LWORK >= max(1,N). 00090 *> For optimum performance LWORK >= N*NB, where NB is 00091 *> the optimal blocksize. 00092 *> 00093 *> If LWORK = -1, then a workspace query is assumed; the routine 00094 *> only calculates the optimal size of the WORK array, returns 00095 *> this value as the first entry of the WORK array, and no error 00096 *> message related to LWORK is issued by XERBLA. 00097 *> \endverbatim 00098 *> 00099 *> \param[out] INFO 00100 *> \verbatim 00101 *> INFO is INTEGER 00102 *> = 0: successful exit 00103 *> < 0: if INFO = -i, the i-th argument had an illegal value 00104 *> \endverbatim 00105 * 00106 * Authors: 00107 * ======== 00108 * 00109 *> \author Univ. of Tennessee 00110 *> \author Univ. of California Berkeley 00111 *> \author Univ. of Colorado Denver 00112 *> \author NAG Ltd. 00113 * 00114 *> \date November 2011 00115 * 00116 *> \ingroup doubleGEcomputational 00117 * 00118 *> \par Further Details: 00119 * ===================== 00120 *> 00121 *> \verbatim 00122 *> 00123 *> The matrix Q is represented as a product of elementary reflectors 00124 *> 00125 *> Q = H(1) H(2) . . . H(k), where k = min(m,n). 00126 *> 00127 *> Each H(i) has the form 00128 *> 00129 *> H(i) = I - tau * v * v**T 00130 *> 00131 *> where tau is a real scalar, and v is a real vector with 00132 *> v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), 00133 *> and tau in TAU(i). 00134 *> \endverbatim 00135 *> 00136 * ===================================================================== 00137 SUBROUTINE DGEQRFP( M, N, A, LDA, TAU, WORK, LWORK, INFO ) 00138 * 00139 * -- LAPACK computational routine (version 3.4.0) -- 00140 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00141 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00142 * November 2011 00143 * 00144 * .. Scalar Arguments .. 00145 INTEGER INFO, LDA, LWORK, M, N 00146 * .. 00147 * .. Array Arguments .. 00148 DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) 00149 * .. 00150 * 00151 * ===================================================================== 00152 * 00153 * .. Local Scalars .. 00154 LOGICAL LQUERY 00155 INTEGER I, IB, IINFO, IWS, K, LDWORK, LWKOPT, NB, 00156 $ NBMIN, NX 00157 * .. 00158 * .. External Subroutines .. 00159 EXTERNAL DGEQR2P, DLARFB, DLARFT, XERBLA 00160 * .. 00161 * .. Intrinsic Functions .. 00162 INTRINSIC MAX, MIN 00163 * .. 00164 * .. External Functions .. 00165 INTEGER ILAENV 00166 EXTERNAL ILAENV 00167 * .. 00168 * .. Executable Statements .. 00169 * 00170 * Test the input arguments 00171 * 00172 INFO = 0 00173 NB = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) 00174 LWKOPT = N*NB 00175 WORK( 1 ) = LWKOPT 00176 LQUERY = ( LWORK.EQ.-1 ) 00177 IF( M.LT.0 ) THEN 00178 INFO = -1 00179 ELSE IF( N.LT.0 ) THEN 00180 INFO = -2 00181 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00182 INFO = -4 00183 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN 00184 INFO = -7 00185 END IF 00186 IF( INFO.NE.0 ) THEN 00187 CALL XERBLA( 'DGEQRFP', -INFO ) 00188 RETURN 00189 ELSE IF( LQUERY ) THEN 00190 RETURN 00191 END IF 00192 * 00193 * Quick return if possible 00194 * 00195 K = MIN( M, N ) 00196 IF( K.EQ.0 ) THEN 00197 WORK( 1 ) = 1 00198 RETURN 00199 END IF 00200 * 00201 NBMIN = 2 00202 NX = 0 00203 IWS = N 00204 IF( NB.GT.1 .AND. NB.LT.K ) THEN 00205 * 00206 * Determine when to cross over from blocked to unblocked code. 00207 * 00208 NX = MAX( 0, ILAENV( 3, 'DGEQRF', ' ', M, N, -1, -1 ) ) 00209 IF( NX.LT.K ) THEN 00210 * 00211 * Determine if workspace is large enough for blocked code. 00212 * 00213 LDWORK = N 00214 IWS = LDWORK*NB 00215 IF( LWORK.LT.IWS ) THEN 00216 * 00217 * Not enough workspace to use optimal NB: reduce NB and 00218 * determine the minimum value of NB. 00219 * 00220 NB = LWORK / LDWORK 00221 NBMIN = MAX( 2, ILAENV( 2, 'DGEQRF', ' ', M, N, -1, 00222 $ -1 ) ) 00223 END IF 00224 END IF 00225 END IF 00226 * 00227 IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN 00228 * 00229 * Use blocked code initially 00230 * 00231 DO 10 I = 1, K - NX, NB 00232 IB = MIN( K-I+1, NB ) 00233 * 00234 * Compute the QR factorization of the current block 00235 * A(i:m,i:i+ib-1) 00236 * 00237 CALL DGEQR2P( M-I+1, IB, A( I, I ), LDA, TAU( I ), WORK, 00238 $ IINFO ) 00239 IF( I+IB.LE.N ) THEN 00240 * 00241 * Form the triangular factor of the block reflector 00242 * H = H(i) H(i+1) . . . H(i+ib-1) 00243 * 00244 CALL DLARFT( 'Forward', 'Columnwise', M-I+1, IB, 00245 $ A( I, I ), LDA, TAU( I ), WORK, LDWORK ) 00246 * 00247 * Apply H**T to A(i:m,i+ib:n) from the left 00248 * 00249 CALL DLARFB( 'Left', 'Transpose', 'Forward', 00250 $ 'Columnwise', M-I+1, N-I-IB+1, IB, 00251 $ A( I, I ), LDA, WORK, LDWORK, A( I, I+IB ), 00252 $ LDA, WORK( IB+1 ), LDWORK ) 00253 END IF 00254 10 CONTINUE 00255 ELSE 00256 I = 1 00257 END IF 00258 * 00259 * Use unblocked code to factor the last or only block. 00260 * 00261 IF( I.LE.K ) 00262 $ CALL DGEQR2P( M-I+1, N-I+1, A( I, I ), LDA, TAU( I ), WORK, 00263 $ IINFO ) 00264 * 00265 WORK( 1 ) = IWS 00266 RETURN 00267 * 00268 * End of DGEQRFP 00269 * 00270 END