00001 SUBROUTINE DTZRZF( 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 DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) 00013 * .. 00014 * 00015 * Purpose 00016 * ======= 00017 * 00018 * DTZRZF reduces the M-by-N ( M<=N ) real upper trapezoidal matrix A 00019 * to upper triangular form by means of orthogonal 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 orthogonal 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) DOUBLE PRECISION 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 * orthogonal 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) DOUBLE PRECISION array, dimension (M) 00049 * The scalar factors of the elementary reflectors. 00050 * 00051 * WORK (workspace/output) DOUBLE PRECISION 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 DOUBLE PRECISION ZERO 00104 PARAMETER ( ZERO = 0.0D+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 DLARZB, DLARZT, DLATRZ, 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 END IF 00134 * 00135 IF( INFO.EQ.0 ) THEN 00136 IF( M.EQ.0 .OR. M.EQ.N ) THEN 00137 LWKOPT = 1 00138 ELSE 00139 * 00140 * Determine the block size. 00141 * 00142 NB = ILAENV( 1, 'DGERQF', ' ', M, N, -1, -1 ) 00143 LWKOPT = M*NB 00144 END IF 00145 WORK( 1 ) = LWKOPT 00146 * 00147 IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN 00148 INFO = -7 00149 END IF 00150 END IF 00151 * 00152 IF( INFO.NE.0 ) THEN 00153 CALL XERBLA( 'DTZRZF', -INFO ) 00154 RETURN 00155 ELSE IF( LQUERY ) THEN 00156 RETURN 00157 END IF 00158 * 00159 * Quick return if possible 00160 * 00161 IF( M.EQ.0 ) THEN 00162 RETURN 00163 ELSE IF( M.EQ.N ) THEN 00164 DO 10 I = 1, N 00165 TAU( I ) = ZERO 00166 10 CONTINUE 00167 RETURN 00168 END IF 00169 * 00170 NBMIN = 2 00171 NX = 1 00172 IWS = M 00173 IF( NB.GT.1 .AND. NB.LT.M ) THEN 00174 * 00175 * Determine when to cross over from blocked to unblocked code. 00176 * 00177 NX = MAX( 0, ILAENV( 3, 'DGERQF', ' ', M, N, -1, -1 ) ) 00178 IF( NX.LT.M ) THEN 00179 * 00180 * Determine if workspace is large enough for blocked code. 00181 * 00182 LDWORK = M 00183 IWS = LDWORK*NB 00184 IF( LWORK.LT.IWS ) THEN 00185 * 00186 * Not enough workspace to use optimal NB: reduce NB and 00187 * determine the minimum value of NB. 00188 * 00189 NB = LWORK / LDWORK 00190 NBMIN = MAX( 2, ILAENV( 2, 'DGERQF', ' ', M, N, -1, 00191 $ -1 ) ) 00192 END IF 00193 END IF 00194 END IF 00195 * 00196 IF( NB.GE.NBMIN .AND. NB.LT.M .AND. NX.LT.M ) THEN 00197 * 00198 * Use blocked code initially. 00199 * The last kk rows are handled by the block method. 00200 * 00201 M1 = MIN( M+1, N ) 00202 KI = ( ( M-NX-1 ) / NB )*NB 00203 KK = MIN( M, KI+NB ) 00204 * 00205 DO 20 I = M - KK + KI + 1, M - KK + 1, -NB 00206 IB = MIN( M-I+1, NB ) 00207 * 00208 * Compute the TZ factorization of the current block 00209 * A(i:i+ib-1,i:n) 00210 * 00211 CALL DLATRZ( IB, N-I+1, N-M, A( I, I ), LDA, TAU( I ), 00212 $ WORK ) 00213 IF( I.GT.1 ) THEN 00214 * 00215 * Form the triangular factor of the block reflector 00216 * H = H(i+ib-1) . . . H(i+1) H(i) 00217 * 00218 CALL DLARZT( 'Backward', 'Rowwise', N-M, IB, A( I, M1 ), 00219 $ LDA, TAU( I ), WORK, LDWORK ) 00220 * 00221 * Apply H to A(1:i-1,i:n) from the right 00222 * 00223 CALL DLARZB( 'Right', 'No transpose', 'Backward', 00224 $ 'Rowwise', I-1, N-I+1, IB, N-M, A( I, M1 ), 00225 $ LDA, WORK, LDWORK, A( 1, I ), LDA, 00226 $ WORK( IB+1 ), LDWORK ) 00227 END IF 00228 20 CONTINUE 00229 MU = I + NB - 1 00230 ELSE 00231 MU = M 00232 END IF 00233 * 00234 * Use unblocked code to factor the last or only block 00235 * 00236 IF( MU.GT.0 ) 00237 $ CALL DLATRZ( MU, N, N-M, A, LDA, TAU, WORK ) 00238 * 00239 WORK( 1 ) = LWKOPT 00240 * 00241 RETURN 00242 * 00243 * End of DTZRZF 00244 * 00245 END
1.7.2