![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZGET51 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 ZGET51( ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, 00012 * RWORK, RESULT ) 00013 * 00014 * .. Scalar Arguments .. 00015 * INTEGER ITYPE, LDA, LDB, LDU, LDV, N 00016 * DOUBLE PRECISION RESULT 00017 * .. 00018 * .. Array Arguments .. 00019 * DOUBLE PRECISION RWORK( * ) 00020 * COMPLEX*16 A( LDA, * ), B( LDB, * ), U( LDU, * ), 00021 * $ V( LDV, * ), WORK( * ) 00022 * .. 00023 * 00024 * 00025 *> \par Purpose: 00026 * ============= 00027 *> 00028 *> \verbatim 00029 *> 00030 *> ZGET51 generally checks a decomposition of the form 00031 *> 00032 *> A = U B VC> 00033 *> where * means conjugate transpose and U and V are unitary. 00034 *> 00035 *> Specifically, if ITYPE=1 00036 *> 00037 *> RESULT = | A - U B V* | / ( |A| n ulp ) 00038 *> 00039 *> If ITYPE=2, then: 00040 *> 00041 *> RESULT = | A - B | / ( |A| n ulp ) 00042 *> 00043 *> If ITYPE=3, then: 00044 *> 00045 *> RESULT = | I - UU* | / ( n ulp ) 00046 *> \endverbatim 00047 * 00048 * Arguments: 00049 * ========== 00050 * 00051 *> \param[in] ITYPE 00052 *> \verbatim 00053 *> ITYPE is INTEGER 00054 *> Specifies the type of tests to be performed. 00055 *> =1: RESULT = | A - U B V* | / ( |A| n ulp ) 00056 *> =2: RESULT = | A - B | / ( |A| n ulp ) 00057 *> =3: RESULT = | I - UU* | / ( n ulp ) 00058 *> \endverbatim 00059 *> 00060 *> \param[in] N 00061 *> \verbatim 00062 *> N is INTEGER 00063 *> The size of the matrix. If it is zero, ZGET51 does nothing. 00064 *> It must be at least zero. 00065 *> \endverbatim 00066 *> 00067 *> \param[in] A 00068 *> \verbatim 00069 *> A is COMPLEX*16 array, dimension (LDA, N) 00070 *> The original (unfactored) matrix. 00071 *> \endverbatim 00072 *> 00073 *> \param[in] LDA 00074 *> \verbatim 00075 *> LDA is INTEGER 00076 *> The leading dimension of A. It must be at least 1 00077 *> and at least N. 00078 *> \endverbatim 00079 *> 00080 *> \param[in] B 00081 *> \verbatim 00082 *> B is COMPLEX*16 array, dimension (LDB, N) 00083 *> The factored matrix. 00084 *> \endverbatim 00085 *> 00086 *> \param[in] LDB 00087 *> \verbatim 00088 *> LDB is INTEGER 00089 *> The leading dimension of B. It must be at least 1 00090 *> and at least N. 00091 *> \endverbatim 00092 *> 00093 *> \param[in] U 00094 *> \verbatim 00095 *> U is COMPLEX*16 array, dimension (LDU, N) 00096 *> The unitary matrix on the left-hand side in the 00097 *> decomposition. 00098 *> Not referenced if ITYPE=2 00099 *> \endverbatim 00100 *> 00101 *> \param[in] LDU 00102 *> \verbatim 00103 *> LDU is INTEGER 00104 *> The leading dimension of U. LDU must be at least N and 00105 *> at least 1. 00106 *> \endverbatim 00107 *> 00108 *> \param[in] V 00109 *> \verbatim 00110 *> V is COMPLEX*16 array, dimension (LDV, N) 00111 *> The unitary matrix on the left-hand side in the 00112 *> decomposition. 00113 *> Not referenced if ITYPE=2 00114 *> \endverbatim 00115 *> 00116 *> \param[in] LDV 00117 *> \verbatim 00118 *> LDV is INTEGER 00119 *> The leading dimension of V. LDV must be at least N and 00120 *> at least 1. 00121 *> \endverbatim 00122 *> 00123 *> \param[out] WORK 00124 *> \verbatim 00125 *> WORK is COMPLEX*16 array, dimension (2*N**2) 00126 *> \endverbatim 00127 *> 00128 *> \param[out] RWORK 00129 *> \verbatim 00130 *> RWORK is DOUBLE PRECISION array, dimension (N) 00131 *> \endverbatim 00132 *> 00133 *> \param[out] RESULT 00134 *> \verbatim 00135 *> RESULT is DOUBLE PRECISION 00136 *> The values computed by the test specified by ITYPE. The 00137 *> value is currently limited to 1/ulp, to avoid overflow. 00138 *> Errors are flagged by RESULT=10/ulp. 00139 *> \endverbatim 00140 * 00141 * Authors: 00142 * ======== 00143 * 00144 *> \author Univ. of Tennessee 00145 *> \author Univ. of California Berkeley 00146 *> \author Univ. of Colorado Denver 00147 *> \author NAG Ltd. 00148 * 00149 *> \date November 2011 00150 * 00151 *> \ingroup complex16_eig 00152 * 00153 * ===================================================================== 00154 SUBROUTINE ZGET51( ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, 00155 $ RWORK, RESULT ) 00156 * 00157 * -- LAPACK test routine (version 3.4.0) -- 00158 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00159 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00160 * November 2011 00161 * 00162 * .. Scalar Arguments .. 00163 INTEGER ITYPE, LDA, LDB, LDU, LDV, N 00164 DOUBLE PRECISION RESULT 00165 * .. 00166 * .. Array Arguments .. 00167 DOUBLE PRECISION RWORK( * ) 00168 COMPLEX*16 A( LDA, * ), B( LDB, * ), U( LDU, * ), 00169 $ V( LDV, * ), WORK( * ) 00170 * .. 00171 * 00172 * ===================================================================== 00173 * 00174 * .. Parameters .. 00175 DOUBLE PRECISION ZERO, ONE, TEN 00176 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 10.0D+0 ) 00177 COMPLEX*16 CZERO, CONE 00178 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00179 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00180 * .. 00181 * .. Local Scalars .. 00182 INTEGER JCOL, JDIAG, JROW 00183 DOUBLE PRECISION ANORM, ULP, UNFL, WNORM 00184 * .. 00185 * .. External Functions .. 00186 DOUBLE PRECISION DLAMCH, ZLANGE 00187 EXTERNAL DLAMCH, ZLANGE 00188 * .. 00189 * .. External Subroutines .. 00190 EXTERNAL ZGEMM, ZLACPY 00191 * .. 00192 * .. Intrinsic Functions .. 00193 INTRINSIC DBLE, MAX, MIN 00194 * .. 00195 * .. Executable Statements .. 00196 * 00197 RESULT = ZERO 00198 IF( N.LE.0 ) 00199 $ RETURN 00200 * 00201 * Constants 00202 * 00203 UNFL = DLAMCH( 'Safe minimum' ) 00204 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 00205 * 00206 * Some Error Checks 00207 * 00208 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 00209 RESULT = TEN / ULP 00210 RETURN 00211 END IF 00212 * 00213 IF( ITYPE.LE.2 ) THEN 00214 * 00215 * Tests scaled by the norm(A) 00216 * 00217 ANORM = MAX( ZLANGE( '1', N, N, A, LDA, RWORK ), UNFL ) 00218 * 00219 IF( ITYPE.EQ.1 ) THEN 00220 * 00221 * ITYPE=1: Compute W = A - UBV' 00222 * 00223 CALL ZLACPY( ' ', N, N, A, LDA, WORK, N ) 00224 CALL ZGEMM( 'N', 'N', N, N, N, CONE, U, LDU, B, LDB, CZERO, 00225 $ WORK( N**2+1 ), N ) 00226 * 00227 CALL ZGEMM( 'N', 'C', N, N, N, -CONE, WORK( N**2+1 ), N, V, 00228 $ LDV, CONE, WORK, N ) 00229 * 00230 ELSE 00231 * 00232 * ITYPE=2: Compute W = A - B 00233 * 00234 CALL ZLACPY( ' ', N, N, B, LDB, WORK, N ) 00235 * 00236 DO 20 JCOL = 1, N 00237 DO 10 JROW = 1, N 00238 WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) ) 00239 $ - A( JROW, JCOL ) 00240 10 CONTINUE 00241 20 CONTINUE 00242 END IF 00243 * 00244 * Compute norm(W)/ ( ulp*norm(A) ) 00245 * 00246 WNORM = ZLANGE( '1', N, N, WORK, N, RWORK ) 00247 * 00248 IF( ANORM.GT.WNORM ) THEN 00249 RESULT = ( WNORM / ANORM ) / ( N*ULP ) 00250 ELSE 00251 IF( ANORM.LT.ONE ) THEN 00252 RESULT = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP ) 00253 ELSE 00254 RESULT = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP ) 00255 END IF 00256 END IF 00257 * 00258 ELSE 00259 * 00260 * Tests not scaled by norm(A) 00261 * 00262 * ITYPE=3: Compute UU' - I 00263 * 00264 CALL ZGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO, 00265 $ WORK, N ) 00266 * 00267 DO 30 JDIAG = 1, N 00268 WORK( ( N+1 )*( JDIAG-1 )+1 ) = WORK( ( N+1 )*( JDIAG-1 )+ 00269 $ 1 ) - CONE 00270 30 CONTINUE 00271 * 00272 RESULT = MIN( ZLANGE( '1', N, N, WORK, N, RWORK ), 00273 $ DBLE( N ) ) / ( N*ULP ) 00274 END IF 00275 * 00276 RETURN 00277 * 00278 * End of ZGET51 00279 * 00280 END