![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZTZRZF 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZTZRZF + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztzrzf.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztzrzf.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztzrzf.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZTZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * INTEGER INFO, LDA, LWORK, M, N 00025 * .. 00026 * .. Array Arguments .. 00027 * COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * ) 00028 * .. 00029 * 00030 * 00031 *> \par Purpose: 00032 * ============= 00033 *> 00034 *> \verbatim 00035 *> 00036 *> ZTZRZF reduces the M-by-N ( M<=N ) complex upper trapezoidal matrix A 00037 *> to upper triangular form by means of unitary transformations. 00038 *> 00039 *> The upper trapezoidal matrix A is factored as 00040 *> 00041 *> A = ( R 0 ) * Z, 00042 *> 00043 *> where Z is an N-by-N unitary matrix and R is an M-by-M upper 00044 *> triangular matrix. 00045 *> \endverbatim 00046 * 00047 * Arguments: 00048 * ========== 00049 * 00050 *> \param[in] M 00051 *> \verbatim 00052 *> M is INTEGER 00053 *> The number of rows of the matrix A. M >= 0. 00054 *> \endverbatim 00055 *> 00056 *> \param[in] N 00057 *> \verbatim 00058 *> N is INTEGER 00059 *> The number of columns of the matrix A. N >= M. 00060 *> \endverbatim 00061 *> 00062 *> \param[in,out] A 00063 *> \verbatim 00064 *> A is COMPLEX*16 array, dimension (LDA,N) 00065 *> On entry, the leading M-by-N upper trapezoidal part of the 00066 *> array A must contain the matrix to be factorized. 00067 *> On exit, the leading M-by-M upper triangular part of A 00068 *> contains the upper triangular matrix R, and elements M+1 to 00069 *> N of the first M rows of A, with the array TAU, represent the 00070 *> unitary matrix Z as a product of M elementary reflectors. 00071 *> \endverbatim 00072 *> 00073 *> \param[in] LDA 00074 *> \verbatim 00075 *> LDA is INTEGER 00076 *> The leading dimension of the array A. LDA >= max(1,M). 00077 *> \endverbatim 00078 *> 00079 *> \param[out] TAU 00080 *> \verbatim 00081 *> TAU is COMPLEX*16 array, dimension (M) 00082 *> The scalar factors of the elementary reflectors. 00083 *> \endverbatim 00084 *> 00085 *> \param[out] WORK 00086 *> \verbatim 00087 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) 00088 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00089 *> \endverbatim 00090 *> 00091 *> \param[in] LWORK 00092 *> \verbatim 00093 *> LWORK is INTEGER 00094 *> The dimension of the array WORK. LWORK >= max(1,M). 00095 *> For optimum performance LWORK >= M*NB, where NB is 00096 *> the optimal blocksize. 00097 *> 00098 *> If LWORK = -1, then a workspace query is assumed; the routine 00099 *> only calculates the optimal size of the WORK array, returns 00100 *> this value as the first entry of the WORK array, and no error 00101 *> message related to LWORK is issued by XERBLA. 00102 *> \endverbatim 00103 *> 00104 *> \param[out] INFO 00105 *> \verbatim 00106 *> INFO is INTEGER 00107 *> = 0: successful exit 00108 *> < 0: if INFO = -i, the i-th argument had an illegal value 00109 *> \endverbatim 00110 * 00111 * Authors: 00112 * ======== 00113 * 00114 *> \author Univ. of Tennessee 00115 *> \author Univ. of California Berkeley 00116 *> \author Univ. of Colorado Denver 00117 *> \author NAG Ltd. 00118 * 00119 *> \date April 2012 00120 * 00121 *> \ingroup complex16OTHERcomputational 00122 * 00123 *> \par Contributors: 00124 * ================== 00125 *> 00126 *> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA 00127 * 00128 *> \par Further Details: 00129 * ===================== 00130 *> 00131 *> \verbatim 00132 *> 00133 *> The N-by-N matrix Z can be computed by 00134 *> 00135 *> Z = Z(1)*Z(2)* ... *Z(M) 00136 *> 00137 *> where each N-by-N Z(k) is given by 00138 *> 00139 *> Z(k) = I - tau(k)*v(k)*v(k)**H 00140 *> 00141 *> with v(k) is the kth row vector of the M-by-N matrix 00142 *> 00143 *> V = ( I A(:,M+1:N) ) 00144 *> 00145 *> I is the M-by-M identity matrix, A(:,M+1:N) 00146 *> is the output stored in A on exit from DTZRZF, 00147 *> and tau(k) is the kth element of the array TAU. 00148 *> 00149 *> \endverbatim 00150 *> 00151 * ===================================================================== 00152 SUBROUTINE ZTZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO ) 00153 * 00154 * -- LAPACK computational routine (version 3.4.1) -- 00155 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00156 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00157 * April 2012 00158 * 00159 * .. Scalar Arguments .. 00160 INTEGER INFO, LDA, LWORK, M, N 00161 * .. 00162 * .. Array Arguments .. 00163 COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * ) 00164 * .. 00165 * 00166 * ===================================================================== 00167 * 00168 * .. Parameters .. 00169 COMPLEX*16 ZERO 00170 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) 00171 * .. 00172 * .. Local Scalars .. 00173 LOGICAL LQUERY 00174 INTEGER I, IB, IWS, KI, KK, LDWORK, LWKMIN, LWKOPT, 00175 $ M1, MU, NB, NBMIN, NX 00176 * .. 00177 * .. External Subroutines .. 00178 EXTERNAL XERBLA, ZLARZB, ZLARZT, ZLATRZ 00179 * .. 00180 * .. Intrinsic Functions .. 00181 INTRINSIC MAX, MIN 00182 * .. 00183 * .. External Functions .. 00184 INTEGER ILAENV 00185 EXTERNAL ILAENV 00186 * .. 00187 * .. Executable Statements .. 00188 * 00189 * Test the input arguments 00190 * 00191 INFO = 0 00192 LQUERY = ( LWORK.EQ.-1 ) 00193 IF( M.LT.0 ) THEN 00194 INFO = -1 00195 ELSE IF( N.LT.M ) THEN 00196 INFO = -2 00197 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00198 INFO = -4 00199 END IF 00200 * 00201 IF( INFO.EQ.0 ) THEN 00202 IF( M.EQ.0 .OR. M.EQ.N ) THEN 00203 LWKOPT = 1 00204 LWKMIN = 1 00205 ELSE 00206 * 00207 * Determine the block size. 00208 * 00209 NB = ILAENV( 1, 'ZGERQF', ' ', M, N, -1, -1 ) 00210 LWKOPT = M*NB 00211 LWKMIN = MAX( 1, M ) 00212 END IF 00213 WORK( 1 ) = LWKOPT 00214 * 00215 IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN 00216 INFO = -7 00217 END IF 00218 END IF 00219 * 00220 IF( INFO.NE.0 ) THEN 00221 CALL XERBLA( 'ZTZRZF', -INFO ) 00222 RETURN 00223 ELSE IF( LQUERY ) THEN 00224 RETURN 00225 END IF 00226 * 00227 * Quick return if possible 00228 * 00229 IF( M.EQ.0 ) THEN 00230 RETURN 00231 ELSE IF( M.EQ.N ) THEN 00232 DO 10 I = 1, N 00233 TAU( I ) = ZERO 00234 10 CONTINUE 00235 RETURN 00236 END IF 00237 * 00238 NBMIN = 2 00239 NX = 1 00240 IWS = M 00241 IF( NB.GT.1 .AND. NB.LT.M ) THEN 00242 * 00243 * Determine when to cross over from blocked to unblocked code. 00244 * 00245 NX = MAX( 0, ILAENV( 3, 'ZGERQF', ' ', M, N, -1, -1 ) ) 00246 IF( NX.LT.M ) THEN 00247 * 00248 * Determine if workspace is large enough for blocked code. 00249 * 00250 LDWORK = M 00251 IWS = LDWORK*NB 00252 IF( LWORK.LT.IWS ) THEN 00253 * 00254 * Not enough workspace to use optimal NB: reduce NB and 00255 * determine the minimum value of NB. 00256 * 00257 NB = LWORK / LDWORK 00258 NBMIN = MAX( 2, ILAENV( 2, 'ZGERQF', ' ', M, N, -1, 00259 $ -1 ) ) 00260 END IF 00261 END IF 00262 END IF 00263 * 00264 IF( NB.GE.NBMIN .AND. NB.LT.M .AND. NX.LT.M ) THEN 00265 * 00266 * Use blocked code initially. 00267 * The last kk rows are handled by the block method. 00268 * 00269 M1 = MIN( M+1, N ) 00270 KI = ( ( M-NX-1 ) / NB )*NB 00271 KK = MIN( M, KI+NB ) 00272 * 00273 DO 20 I = M - KK + KI + 1, M - KK + 1, -NB 00274 IB = MIN( M-I+1, NB ) 00275 * 00276 * Compute the TZ factorization of the current block 00277 * A(i:i+ib-1,i:n) 00278 * 00279 CALL ZLATRZ( IB, N-I+1, N-M, A( I, I ), LDA, TAU( I ), 00280 $ WORK ) 00281 IF( I.GT.1 ) THEN 00282 * 00283 * Form the triangular factor of the block reflector 00284 * H = H(i+ib-1) . . . H(i+1) H(i) 00285 * 00286 CALL ZLARZT( 'Backward', 'Rowwise', N-M, IB, A( I, M1 ), 00287 $ LDA, TAU( I ), WORK, LDWORK ) 00288 * 00289 * Apply H to A(1:i-1,i:n) from the right 00290 * 00291 CALL ZLARZB( 'Right', 'No transpose', 'Backward', 00292 $ 'Rowwise', I-1, N-I+1, IB, N-M, A( I, M1 ), 00293 $ LDA, WORK, LDWORK, A( 1, I ), LDA, 00294 $ WORK( IB+1 ), LDWORK ) 00295 END IF 00296 20 CONTINUE 00297 MU = I + NB - 1 00298 ELSE 00299 MU = M 00300 END IF 00301 * 00302 * Use unblocked code to factor the last or only block 00303 * 00304 IF( MU.GT.0 ) 00305 $ CALL ZLATRZ( MU, N, N-M, A, LDA, TAU, WORK ) 00306 * 00307 WORK( 1 ) = LWKOPT 00308 * 00309 RETURN 00310 * 00311 * End of ZTZRZF 00312 * 00313 END