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