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