![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZQLT03 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE ZQLT03( M, N, K, AF, C, CC, Q, LDA, TAU, WORK, LWORK, 00012 * RWORK, RESULT ) 00013 * 00014 * .. Scalar Arguments .. 00015 * INTEGER K, LDA, LWORK, M, N 00016 * .. 00017 * .. Array Arguments .. 00018 * DOUBLE PRECISION RESULT( * ), RWORK( * ) 00019 * COMPLEX*16 AF( LDA, * ), C( LDA, * ), CC( LDA, * ), 00020 * $ Q( LDA, * ), TAU( * ), WORK( LWORK ) 00021 * .. 00022 * 00023 * 00024 *> \par Purpose: 00025 * ============= 00026 *> 00027 *> \verbatim 00028 *> 00029 *> ZQLT03 tests ZUNMQL, which computes Q*C, Q'*C, C*Q or C*Q'. 00030 *> 00031 *> ZQLT03 compares the results of a call to ZUNMQL with the results of 00032 *> forming Q explicitly by a call to ZUNGQL and then performing matrix 00033 *> multiplication by a call to ZGEMM. 00034 *> \endverbatim 00035 * 00036 * Arguments: 00037 * ========== 00038 * 00039 *> \param[in] M 00040 *> \verbatim 00041 *> M is INTEGER 00042 *> The order of the orthogonal matrix Q. M >= 0. 00043 *> \endverbatim 00044 *> 00045 *> \param[in] N 00046 *> \verbatim 00047 *> N is INTEGER 00048 *> The number of rows or columns of the matrix C; C is m-by-n if 00049 *> Q is applied from the left, or n-by-m if Q is applied from 00050 *> the right. N >= 0. 00051 *> \endverbatim 00052 *> 00053 *> \param[in] K 00054 *> \verbatim 00055 *> K is INTEGER 00056 *> The number of elementary reflectors whose product defines the 00057 *> orthogonal matrix Q. M >= K >= 0. 00058 *> \endverbatim 00059 *> 00060 *> \param[in] AF 00061 *> \verbatim 00062 *> AF is COMPLEX*16 array, dimension (LDA,N) 00063 *> Details of the QL factorization of an m-by-n matrix, as 00064 *> returned by ZGEQLF. See CGEQLF for further details. 00065 *> \endverbatim 00066 *> 00067 *> \param[out] C 00068 *> \verbatim 00069 *> C is COMPLEX*16 array, dimension (LDA,N) 00070 *> \endverbatim 00071 *> 00072 *> \param[out] CC 00073 *> \verbatim 00074 *> CC is COMPLEX*16 array, dimension (LDA,N) 00075 *> \endverbatim 00076 *> 00077 *> \param[out] Q 00078 *> \verbatim 00079 *> Q is COMPLEX*16 array, dimension (LDA,M) 00080 *> \endverbatim 00081 *> 00082 *> \param[in] LDA 00083 *> \verbatim 00084 *> LDA is INTEGER 00085 *> The leading dimension of the arrays AF, C, CC, and Q. 00086 *> \endverbatim 00087 *> 00088 *> \param[in] TAU 00089 *> \verbatim 00090 *> TAU is COMPLEX*16 array, dimension (min(M,N)) 00091 *> The scalar factors of the elementary reflectors corresponding 00092 *> to the QL factorization in AF. 00093 *> \endverbatim 00094 *> 00095 *> \param[out] WORK 00096 *> \verbatim 00097 *> WORK is COMPLEX*16 array, dimension (LWORK) 00098 *> \endverbatim 00099 *> 00100 *> \param[in] LWORK 00101 *> \verbatim 00102 *> LWORK is INTEGER 00103 *> The length of WORK. LWORK must be at least M, and should be 00104 *> M*NB, where NB is the blocksize for this environment. 00105 *> \endverbatim 00106 *> 00107 *> \param[out] RWORK 00108 *> \verbatim 00109 *> RWORK is DOUBLE PRECISION array, dimension (M) 00110 *> \endverbatim 00111 *> 00112 *> \param[out] RESULT 00113 *> \verbatim 00114 *> RESULT is DOUBLE PRECISION array, dimension (4) 00115 *> The test ratios compare two techniques for multiplying a 00116 *> random matrix C by an m-by-m orthogonal matrix Q. 00117 *> RESULT(1) = norm( Q*C - Q*C ) / ( M * norm(C) * EPS ) 00118 *> RESULT(2) = norm( C*Q - C*Q ) / ( M * norm(C) * EPS ) 00119 *> RESULT(3) = norm( Q'*C - Q'*C )/ ( M * norm(C) * EPS ) 00120 *> RESULT(4) = norm( C*Q' - C*Q' )/ ( M * norm(C) * EPS ) 00121 *> \endverbatim 00122 * 00123 * Authors: 00124 * ======== 00125 * 00126 *> \author Univ. of Tennessee 00127 *> \author Univ. of California Berkeley 00128 *> \author Univ. of Colorado Denver 00129 *> \author NAG Ltd. 00130 * 00131 *> \date November 2011 00132 * 00133 *> \ingroup complex16_lin 00134 * 00135 * ===================================================================== 00136 SUBROUTINE ZQLT03( M, N, K, AF, C, CC, Q, LDA, TAU, WORK, LWORK, 00137 $ RWORK, RESULT ) 00138 * 00139 * -- LAPACK test 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 K, LDA, LWORK, M, N 00146 * .. 00147 * .. Array Arguments .. 00148 DOUBLE PRECISION RESULT( * ), RWORK( * ) 00149 COMPLEX*16 AF( LDA, * ), C( LDA, * ), CC( LDA, * ), 00150 $ Q( LDA, * ), TAU( * ), WORK( LWORK ) 00151 * .. 00152 * 00153 * ===================================================================== 00154 * 00155 * .. Parameters .. 00156 DOUBLE PRECISION ZERO, ONE 00157 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00158 COMPLEX*16 ROGUE 00159 PARAMETER ( ROGUE = ( -1.0D+10, -1.0D+10 ) ) 00160 * .. 00161 * .. Local Scalars .. 00162 CHARACTER SIDE, TRANS 00163 INTEGER INFO, ISIDE, ITRANS, J, MC, MINMN, NC 00164 DOUBLE PRECISION CNORM, EPS, RESID 00165 * .. 00166 * .. External Functions .. 00167 LOGICAL LSAME 00168 DOUBLE PRECISION DLAMCH, ZLANGE 00169 EXTERNAL LSAME, DLAMCH, ZLANGE 00170 * .. 00171 * .. External Subroutines .. 00172 EXTERNAL ZGEMM, ZLACPY, ZLARNV, ZLASET, ZUNGQL, ZUNMQL 00173 * .. 00174 * .. Local Arrays .. 00175 INTEGER ISEED( 4 ) 00176 * .. 00177 * .. Intrinsic Functions .. 00178 INTRINSIC DBLE, DCMPLX, MAX, MIN 00179 * .. 00180 * .. Scalars in Common .. 00181 CHARACTER*32 SRNAMT 00182 * .. 00183 * .. Common blocks .. 00184 COMMON / SRNAMC / SRNAMT 00185 * .. 00186 * .. Data statements .. 00187 DATA ISEED / 1988, 1989, 1990, 1991 / 00188 * .. 00189 * .. Executable Statements .. 00190 * 00191 EPS = DLAMCH( 'Epsilon' ) 00192 MINMN = MIN( M, N ) 00193 * 00194 * Quick return if possible 00195 * 00196 IF( MINMN.EQ.0 ) THEN 00197 RESULT( 1 ) = ZERO 00198 RESULT( 2 ) = ZERO 00199 RESULT( 3 ) = ZERO 00200 RESULT( 4 ) = ZERO 00201 RETURN 00202 END IF 00203 * 00204 * Copy the last k columns of the factorization to the array Q 00205 * 00206 CALL ZLASET( 'Full', M, M, ROGUE, ROGUE, Q, LDA ) 00207 IF( K.GT.0 .AND. M.GT.K ) 00208 $ CALL ZLACPY( 'Full', M-K, K, AF( 1, N-K+1 ), LDA, 00209 $ Q( 1, M-K+1 ), LDA ) 00210 IF( K.GT.1 ) 00211 $ CALL ZLACPY( 'Upper', K-1, K-1, AF( M-K+1, N-K+2 ), LDA, 00212 $ Q( M-K+1, M-K+2 ), LDA ) 00213 * 00214 * Generate the m-by-m matrix Q 00215 * 00216 SRNAMT = 'ZUNGQL' 00217 CALL ZUNGQL( M, M, K, Q, LDA, TAU( MINMN-K+1 ), WORK, LWORK, 00218 $ INFO ) 00219 * 00220 DO 30 ISIDE = 1, 2 00221 IF( ISIDE.EQ.1 ) THEN 00222 SIDE = 'L' 00223 MC = M 00224 NC = N 00225 ELSE 00226 SIDE = 'R' 00227 MC = N 00228 NC = M 00229 END IF 00230 * 00231 * Generate MC by NC matrix C 00232 * 00233 DO 10 J = 1, NC 00234 CALL ZLARNV( 2, ISEED, MC, C( 1, J ) ) 00235 10 CONTINUE 00236 CNORM = ZLANGE( '1', MC, NC, C, LDA, RWORK ) 00237 IF( CNORM.EQ.ZERO ) 00238 $ CNORM = ONE 00239 * 00240 DO 20 ITRANS = 1, 2 00241 IF( ITRANS.EQ.1 ) THEN 00242 TRANS = 'N' 00243 ELSE 00244 TRANS = 'C' 00245 END IF 00246 * 00247 * Copy C 00248 * 00249 CALL ZLACPY( 'Full', MC, NC, C, LDA, CC, LDA ) 00250 * 00251 * Apply Q or Q' to C 00252 * 00253 SRNAMT = 'ZUNMQL' 00254 IF( K.GT.0 ) 00255 $ CALL ZUNMQL( SIDE, TRANS, MC, NC, K, AF( 1, N-K+1 ), LDA, 00256 $ TAU( MINMN-K+1 ), CC, LDA, WORK, LWORK, 00257 $ INFO ) 00258 * 00259 * Form explicit product and subtract 00260 * 00261 IF( LSAME( SIDE, 'L' ) ) THEN 00262 CALL ZGEMM( TRANS, 'No transpose', MC, NC, MC, 00263 $ DCMPLX( -ONE ), Q, LDA, C, LDA, 00264 $ DCMPLX( ONE ), CC, LDA ) 00265 ELSE 00266 CALL ZGEMM( 'No transpose', TRANS, MC, NC, NC, 00267 $ DCMPLX( -ONE ), C, LDA, Q, LDA, 00268 $ DCMPLX( ONE ), CC, LDA ) 00269 END IF 00270 * 00271 * Compute error in the difference 00272 * 00273 RESID = ZLANGE( '1', MC, NC, CC, LDA, RWORK ) 00274 RESULT( ( ISIDE-1 )*2+ITRANS ) = RESID / 00275 $ ( DBLE( MAX( 1, M ) )*CNORM*EPS ) 00276 * 00277 20 CONTINUE 00278 30 CONTINUE 00279 * 00280 RETURN 00281 * 00282 * End of ZQLT03 00283 * 00284 END