![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DGET51 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 DGET51( ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, 00012 * RESULT ) 00013 * 00014 * .. Scalar Arguments .. 00015 * INTEGER ITYPE, LDA, LDB, LDU, LDV, N 00016 * DOUBLE PRECISION RESULT 00017 * .. 00018 * .. Array Arguments .. 00019 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), U( LDU, * ), 00020 * $ V( LDV, * ), WORK( * ) 00021 * .. 00022 * 00023 * 00024 *> \par Purpose: 00025 * ============= 00026 *> 00027 *> \verbatim 00028 *> 00029 *> DGET51 generally checks a decomposition of the form 00030 *> 00031 *> A = U B V' 00032 *> 00033 *> where ' means transpose and U and V are orthogonal. 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, DGET51 does nothing. 00064 *> It must be at least zero. 00065 *> \endverbatim 00066 *> 00067 *> \param[in] A 00068 *> \verbatim 00069 *> A is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDU, N) 00096 *> The orthogonal 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 DOUBLE PRECISION array, dimension (LDV, N) 00111 *> The orthogonal 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 DOUBLE PRECISION array, dimension (2*N**2) 00126 *> \endverbatim 00127 *> 00128 *> \param[out] RESULT 00129 *> \verbatim 00130 *> RESULT is DOUBLE PRECISION 00131 *> The values computed by the test specified by ITYPE. The 00132 *> value is currently limited to 1/ulp, to avoid overflow. 00133 *> Errors are flagged by RESULT=10/ulp. 00134 *> \endverbatim 00135 * 00136 * Authors: 00137 * ======== 00138 * 00139 *> \author Univ. of Tennessee 00140 *> \author Univ. of California Berkeley 00141 *> \author Univ. of Colorado Denver 00142 *> \author NAG Ltd. 00143 * 00144 *> \date November 2011 00145 * 00146 *> \ingroup double_eig 00147 * 00148 * ===================================================================== 00149 SUBROUTINE DGET51( ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, 00150 $ RESULT ) 00151 * 00152 * -- LAPACK test routine (version 3.4.0) -- 00153 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00154 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00155 * November 2011 00156 * 00157 * .. Scalar Arguments .. 00158 INTEGER ITYPE, LDA, LDB, LDU, LDV, N 00159 DOUBLE PRECISION RESULT 00160 * .. 00161 * .. Array Arguments .. 00162 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), U( LDU, * ), 00163 $ V( LDV, * ), WORK( * ) 00164 * .. 00165 * 00166 * ===================================================================== 00167 * 00168 * .. Parameters .. 00169 DOUBLE PRECISION ZERO, ONE, TEN 00170 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0 ) 00171 * .. 00172 * .. Local Scalars .. 00173 INTEGER JCOL, JDIAG, JROW 00174 DOUBLE PRECISION ANORM, ULP, UNFL, WNORM 00175 * .. 00176 * .. External Functions .. 00177 DOUBLE PRECISION DLAMCH, DLANGE 00178 EXTERNAL DLAMCH, DLANGE 00179 * .. 00180 * .. External Subroutines .. 00181 EXTERNAL DGEMM, DLACPY 00182 * .. 00183 * .. Intrinsic Functions .. 00184 INTRINSIC DBLE, MAX, MIN 00185 * .. 00186 * .. Executable Statements .. 00187 * 00188 RESULT = ZERO 00189 IF( N.LE.0 ) 00190 $ RETURN 00191 * 00192 * Constants 00193 * 00194 UNFL = DLAMCH( 'Safe minimum' ) 00195 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 00196 * 00197 * Some Error Checks 00198 * 00199 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 00200 RESULT = TEN / ULP 00201 RETURN 00202 END IF 00203 * 00204 IF( ITYPE.LE.2 ) THEN 00205 * 00206 * Tests scaled by the norm(A) 00207 * 00208 ANORM = MAX( DLANGE( '1', N, N, A, LDA, WORK ), UNFL ) 00209 * 00210 IF( ITYPE.EQ.1 ) THEN 00211 * 00212 * ITYPE=1: Compute W = A - UBV' 00213 * 00214 CALL DLACPY( ' ', N, N, A, LDA, WORK, N ) 00215 CALL DGEMM( 'N', 'N', N, N, N, ONE, U, LDU, B, LDB, ZERO, 00216 $ WORK( N**2+1 ), N ) 00217 * 00218 CALL DGEMM( 'N', 'C', N, N, N, -ONE, WORK( N**2+1 ), N, V, 00219 $ LDV, ONE, WORK, N ) 00220 * 00221 ELSE 00222 * 00223 * ITYPE=2: Compute W = A - B 00224 * 00225 CALL DLACPY( ' ', N, N, B, LDB, WORK, N ) 00226 * 00227 DO 20 JCOL = 1, N 00228 DO 10 JROW = 1, N 00229 WORK( JROW+N*( JCOL-1 ) ) = WORK( JROW+N*( JCOL-1 ) ) 00230 $ - A( JROW, JCOL ) 00231 10 CONTINUE 00232 20 CONTINUE 00233 END IF 00234 * 00235 * Compute norm(W)/ ( ulp*norm(A) ) 00236 * 00237 WNORM = DLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) ) 00238 * 00239 IF( ANORM.GT.WNORM ) THEN 00240 RESULT = ( WNORM / ANORM ) / ( N*ULP ) 00241 ELSE 00242 IF( ANORM.LT.ONE ) THEN 00243 RESULT = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP ) 00244 ELSE 00245 RESULT = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP ) 00246 END IF 00247 END IF 00248 * 00249 ELSE 00250 * 00251 * Tests not scaled by norm(A) 00252 * 00253 * ITYPE=3: Compute UU' - I 00254 * 00255 CALL DGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK, 00256 $ N ) 00257 * 00258 DO 30 JDIAG = 1, N 00259 WORK( ( N+1 )*( JDIAG-1 )+1 ) = WORK( ( N+1 )*( JDIAG-1 )+ 00260 $ 1 ) - ONE 00261 30 CONTINUE 00262 * 00263 RESULT = MIN( DLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) ), 00264 $ DBLE( N ) ) / ( N*ULP ) 00265 END IF 00266 * 00267 RETURN 00268 * 00269 * End of DGET51 00270 * 00271 END