![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CGET54 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 CGET54( 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 * REAL RESULT 00017 * .. 00018 * .. Array Arguments .. 00019 * COMPLEX 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 *> CGET54 checks a generalized decomposition of the form 00031 *> 00032 *> A = U*S*V' and B = U*T* V' 00033 *> 00034 *> where ' means conjugate transpose and U and V are unitary. 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, SGET54 does nothing. 00048 *> It must be at least zero. 00049 *> \endverbatim 00050 *> 00051 *> \param[in] A 00052 *> \verbatim 00053 *> A is COMPLEX 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 COMPLEX array, dimension (3*N**2) 00134 *> \endverbatim 00135 *> 00136 *> \param[out] RESULT 00137 *> \verbatim 00138 *> RESULT is REAL 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 complex_eig 00154 * 00155 * ===================================================================== 00156 SUBROUTINE CGET54( 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 REAL RESULT 00167 * .. 00168 * .. Array Arguments .. 00169 COMPLEX A( LDA, * ), B( LDB, * ), S( LDS, * ), 00170 $ T( LDT, * ), U( LDU, * ), V( LDV, * ), 00171 $ WORK( * ) 00172 * .. 00173 * 00174 * ===================================================================== 00175 * 00176 * .. Parameters .. 00177 REAL ZERO, ONE 00178 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00179 COMPLEX CZERO, CONE 00180 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00181 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00182 * .. 00183 * .. Local Scalars .. 00184 REAL ABNORM, ULP, UNFL, WNORM 00185 * .. 00186 * .. Local Arrays .. 00187 REAL DUM( 1 ) 00188 * .. 00189 * .. External Functions .. 00190 REAL CLANGE, SLAMCH 00191 EXTERNAL CLANGE, SLAMCH 00192 * .. 00193 * .. External Subroutines .. 00194 EXTERNAL CGEMM, CLACPY 00195 * .. 00196 * .. Intrinsic Functions .. 00197 INTRINSIC MAX, MIN, REAL 00198 * .. 00199 * .. Executable Statements .. 00200 * 00201 RESULT = ZERO 00202 IF( N.LE.0 ) 00203 $ RETURN 00204 * 00205 * Constants 00206 * 00207 UNFL = SLAMCH( 'Safe minimum' ) 00208 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 00209 * 00210 * compute the norm of (A,B) 00211 * 00212 CALL CLACPY( 'Full', N, N, A, LDA, WORK, N ) 00213 CALL CLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N ) 00214 ABNORM = MAX( CLANGE( '1', N, 2*N, WORK, N, DUM ), UNFL ) 00215 * 00216 * Compute W1 = A - U*S*V', and put in the array WORK(1:N*N) 00217 * 00218 CALL CLACPY( ' ', N, N, A, LDA, WORK, N ) 00219 CALL CGEMM( 'N', 'N', N, N, N, CONE, U, LDU, S, LDS, CZERO, 00220 $ WORK( N*N+1 ), N ) 00221 * 00222 CALL CGEMM( 'N', 'C', N, N, N, -CONE, WORK( N*N+1 ), N, V, LDV, 00223 $ CONE, WORK, N ) 00224 * 00225 * Compute W2 = B - U*T*V', and put in the workarray W(N*N+1:2*N*N) 00226 * 00227 CALL CLACPY( ' ', N, N, B, LDB, WORK( N*N+1 ), N ) 00228 CALL CGEMM( 'N', 'N', N, N, N, CONE, U, LDU, T, LDT, CZERO, 00229 $ WORK( 2*N*N+1 ), N ) 00230 * 00231 CALL CGEMM( 'N', 'C', N, N, N, -CONE, WORK( 2*N*N+1 ), N, V, LDV, 00232 $ CONE, WORK( N*N+1 ), N ) 00233 * 00234 * Compute norm(W)/ ( ulp*norm((A,B)) ) 00235 * 00236 WNORM = CLANGE( '1', N, 2*N, WORK, N, DUM ) 00237 * 00238 IF( ABNORM.GT.WNORM ) THEN 00239 RESULT = ( WNORM / ABNORM ) / ( 2*N*ULP ) 00240 ELSE 00241 IF( ABNORM.LT.ONE ) THEN 00242 RESULT = ( MIN( WNORM, 2*N*ABNORM ) / ABNORM ) / ( 2*N*ULP ) 00243 ELSE 00244 RESULT = MIN( WNORM / ABNORM, REAL( 2*N ) ) / ( 2*N*ULP ) 00245 END IF 00246 END IF 00247 * 00248 RETURN 00249 * 00250 * End of CGET54 00251 * 00252 END