00001 SUBROUTINE CTZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.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 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 INTEGER INFO, LDA, LWORK, M, N 00010 * .. 00011 * .. Array Arguments .. 00012 COMPLEX A( LDA, * ), TAU( * ), WORK( * ) 00013 * .. 00014 * 00015 * Purpose 00016 * ======= 00017 * 00018 * CTZRZF reduces the M-by-N ( M<=N ) complex upper trapezoidal matrix A 00019 * to upper triangular form by means of unitary transformations. 00020 * 00021 * The upper trapezoidal matrix A is factored as 00022 * 00023 * A = ( R 0 ) * Z, 00024 * 00025 * where Z is an N-by-N unitary matrix and R is an M-by-M upper 00026 * triangular matrix. 00027 * 00028 * Arguments 00029 * ========= 00030 * 00031 * M (input) INTEGER 00032 * The number of rows of the matrix A. M >= 0. 00033 * 00034 * N (input) INTEGER 00035 * The number of columns of the matrix A. N >= M. 00036 * 00037 * A (input/output) COMPLEX array, dimension (LDA,N) 00038 * On entry, the leading M-by-N upper trapezoidal part of the 00039 * array A must contain the matrix to be factorized. 00040 * On exit, the leading M-by-M upper triangular part of A 00041 * contains the upper triangular matrix R, and elements M+1 to 00042 * N of the first M rows of A, with the array TAU, represent the 00043 * unitary matrix Z as a product of M elementary reflectors. 00044 * 00045 * LDA (input) INTEGER 00046 * The leading dimension of the array A. LDA >= max(1,M). 00047 * 00048 * TAU (output) COMPLEX array, dimension (M) 00049 * The scalar factors of the elementary reflectors. 00050 * 00051 * WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) 00052 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00053 * 00054 * LWORK (input) INTEGER 00055 * The dimension of the array WORK. LWORK >= max(1,M). 00056 * For optimum performance LWORK >= M*NB, where NB is 00057 * the optimal blocksize. 00058 * 00059 * If LWORK = -1, then a workspace query is assumed; the routine 00060 * only calculates the optimal size of the WORK array, returns 00061 * this value as the first entry of the WORK array, and no error 00062 * message related to LWORK is issued by XERBLA. 00063 * 00064 * INFO (output) INTEGER 00065 * = 0: successful exit 00066 * < 0: if INFO = -i, the i-th argument had an illegal value 00067 * 00068 * Further Details 00069 * =============== 00070 * 00071 * Based on contributions by 00072 * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA 00073 * 00074 * The factorization is obtained by Householder's method. The kth 00075 * transformation matrix, Z( k ), which is used to introduce zeros into 00076 * the ( m - k + 1 )th row of A, is given in the form 00077 * 00078 * Z( k ) = ( I 0 ), 00079 * ( 0 T( k ) ) 00080 * 00081 * where 00082 * 00083 * T( k ) = I - tau*u( k )*u( k )', u( k ) = ( 1 ), 00084 * ( 0 ) 00085 * ( z( k ) ) 00086 * 00087 * tau is a scalar and z( k ) is an ( n - m ) element vector. 00088 * tau and z( k ) are chosen to annihilate the elements of the kth row 00089 * of X. 00090 * 00091 * The scalar tau is returned in the kth element of TAU and the vector 00092 * u( k ) in the kth row of A, such that the elements of z( k ) are 00093 * in a( k, m + 1 ), ..., a( k, n ). The elements of R are returned in 00094 * the upper triangular part of A. 00095 * 00096 * Z is given by 00097 * 00098 * Z = Z( 1 ) * Z( 2 ) * ... * Z( m ). 00099 * 00100 * ===================================================================== 00101 * 00102 * .. Parameters .. 00103 COMPLEX ZERO 00104 PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ) ) 00105 * .. 00106 * .. Local Scalars .. 00107 LOGICAL LQUERY 00108 INTEGER I, IB, IWS, KI, KK, LDWORK, LWKOPT, M1, MU, NB, 00109 $ NBMIN, NX 00110 * .. 00111 * .. External Subroutines .. 00112 EXTERNAL CLARZB, CLARZT, CLATRZ, XERBLA 00113 * .. 00114 * .. Intrinsic Functions .. 00115 INTRINSIC MAX, MIN 00116 * .. 00117 * .. External Functions .. 00118 INTEGER ILAENV 00119 EXTERNAL ILAENV 00120 * .. 00121 * .. Executable Statements .. 00122 * 00123 * Test the input arguments 00124 * 00125 INFO = 0 00126 LQUERY = ( LWORK.EQ.-1 ) 00127 IF( M.LT.0 ) THEN 00128 INFO = -1 00129 ELSE IF( N.LT.M ) THEN 00130 INFO = -2 00131 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00132 INFO = -4 00133 ELSE IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN 00134 INFO = -7 00135 END IF 00136 * 00137 IF( INFO.EQ.0 ) THEN 00138 IF( M.EQ.0 .OR. M.EQ.N ) THEN 00139 LWKOPT = 1 00140 ELSE 00141 * 00142 * Determine the block size. 00143 * 00144 NB = ILAENV( 1, 'CGERQF', ' ', M, N, -1, -1 ) 00145 LWKOPT = M*NB 00146 END IF 00147 WORK( 1 ) = LWKOPT 00148 * 00149 IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN 00150 INFO = -7 00151 END IF 00152 END IF 00153 * 00154 IF( INFO.NE.0 ) THEN 00155 CALL XERBLA( 'CTZRZF', -INFO ) 00156 RETURN 00157 ELSE IF( LQUERY ) THEN 00158 RETURN 00159 END IF 00160 * 00161 * Quick return if possible 00162 * 00163 IF( M.EQ.0 ) THEN 00164 RETURN 00165 ELSE IF( M.EQ.N ) THEN 00166 DO 10 I = 1, N 00167 TAU( I ) = ZERO 00168 10 CONTINUE 00169 RETURN 00170 END IF 00171 * 00172 NBMIN = 2 00173 NX = 1 00174 IWS = M 00175 IF( NB.GT.1 .AND. NB.LT.M ) THEN 00176 * 00177 * Determine when to cross over from blocked to unblocked code. 00178 * 00179 NX = MAX( 0, ILAENV( 3, 'CGERQF', ' ', M, N, -1, -1 ) ) 00180 IF( NX.LT.M ) THEN 00181 * 00182 * Determine if workspace is large enough for blocked code. 00183 * 00184 LDWORK = M 00185 IWS = LDWORK*NB 00186 IF( LWORK.LT.IWS ) THEN 00187 * 00188 * Not enough workspace to use optimal NB: reduce NB and 00189 * determine the minimum value of NB. 00190 * 00191 NB = LWORK / LDWORK 00192 NBMIN = MAX( 2, ILAENV( 2, 'CGERQF', ' ', M, N, -1, 00193 $ -1 ) ) 00194 END IF 00195 END IF 00196 END IF 00197 * 00198 IF( NB.GE.NBMIN .AND. NB.LT.M .AND. NX.LT.M ) THEN 00199 * 00200 * Use blocked code initially. 00201 * The last kk rows are handled by the block method. 00202 * 00203 M1 = MIN( M+1, N ) 00204 KI = ( ( M-NX-1 ) / NB )*NB 00205 KK = MIN( M, KI+NB ) 00206 * 00207 DO 20 I = M - KK + KI + 1, M - KK + 1, -NB 00208 IB = MIN( M-I+1, NB ) 00209 * 00210 * Compute the TZ factorization of the current block 00211 * A(i:i+ib-1,i:n) 00212 * 00213 CALL CLATRZ( IB, N-I+1, N-M, A( I, I ), LDA, TAU( I ), 00214 $ WORK ) 00215 IF( I.GT.1 ) THEN 00216 * 00217 * Form the triangular factor of the block reflector 00218 * H = H(i+ib-1) . . . H(i+1) H(i) 00219 * 00220 CALL CLARZT( 'Backward', 'Rowwise', N-M, IB, A( I, M1 ), 00221 $ LDA, TAU( I ), WORK, LDWORK ) 00222 * 00223 * Apply H to A(1:i-1,i:n) from the right 00224 * 00225 CALL CLARZB( 'Right', 'No transpose', 'Backward', 00226 $ 'Rowwise', I-1, N-I+1, IB, N-M, A( I, M1 ), 00227 $ LDA, WORK, LDWORK, A( 1, I ), LDA, 00228 $ WORK( IB+1 ), LDWORK ) 00229 END IF 00230 20 CONTINUE 00231 MU = I + NB - 1 00232 ELSE 00233 MU = M 00234 END IF 00235 * 00236 * Use unblocked code to factor the last or only block 00237 * 00238 IF( MU.GT.0 ) 00239 $ CALL CLATRZ( MU, N, N-M, A, LDA, TAU, WORK ) 00240 * 00241 WORK( 1 ) = LWKOPT 00242 * 00243 RETURN 00244 * 00245 * End of CTZRZF 00246 * 00247 END
1.7.2