00001 SUBROUTINE DGET54( N, A, LDA, B, LDB, S, LDS, T, LDT, U, LDU, V, 00002 $ LDV, WORK, RESULT ) 00003 * 00004 * -- LAPACK test routine (version 3.1) -- 00005 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 INTEGER LDA, LDB, LDS, LDT, LDU, LDV, N 00010 DOUBLE PRECISION RESULT 00011 * .. 00012 * .. Array Arguments .. 00013 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( LDS, * ), 00014 $ T( LDT, * ), U( LDU, * ), V( LDV, * ), 00015 $ WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DGET54 checks a generalized decomposition of the form 00022 * 00023 * A = U*S*V' and B = U*T* V' 00024 * 00025 * where ' means transpose and U and V are orthogonal. 00026 * 00027 * Specifically, 00028 * 00029 * RESULT = ||( A - U*S*V', B - U*T*V' )|| / (||( A, B )||*n*ulp ) 00030 * 00031 * Arguments 00032 * ========= 00033 * 00034 * N (input) INTEGER 00035 * The size of the matrix. If it is zero, DGET54 does nothing. 00036 * It must be at least zero. 00037 * 00038 * A (input) DOUBLE PRECISION array, dimension (LDA, N) 00039 * The original (unfactored) matrix A. 00040 * 00041 * LDA (input) INTEGER 00042 * The leading dimension of A. It must be at least 1 00043 * and at least N. 00044 * 00045 * B (input) DOUBLE PRECISION array, dimension (LDB, N) 00046 * The original (unfactored) matrix B. 00047 * 00048 * LDB (input) INTEGER 00049 * The leading dimension of B. It must be at least 1 00050 * and at least N. 00051 * 00052 * S (input) DOUBLE PRECISION array, dimension (LDS, N) 00053 * The factored matrix S. 00054 * 00055 * LDS (input) INTEGER 00056 * The leading dimension of S. It must be at least 1 00057 * and at least N. 00058 * 00059 * T (input) DOUBLE PRECISION array, dimension (LDT, N) 00060 * The factored matrix T. 00061 * 00062 * LDT (input) INTEGER 00063 * The leading dimension of T. It must be at least 1 00064 * and at least N. 00065 * 00066 * U (input) DOUBLE PRECISION array, dimension (LDU, N) 00067 * The orthogonal matrix on the left-hand side in the 00068 * decomposition. 00069 * 00070 * LDU (input) INTEGER 00071 * The leading dimension of U. LDU must be at least N and 00072 * at least 1. 00073 * 00074 * V (input) DOUBLE PRECISION array, dimension (LDV, N) 00075 * The orthogonal matrix on the left-hand side in the 00076 * decomposition. 00077 * 00078 * LDV (input) INTEGER 00079 * The leading dimension of V. LDV must be at least N and 00080 * at least 1. 00081 * 00082 * WORK (workspace) DOUBLE PRECISION array, dimension (3*N**2) 00083 * 00084 * RESULT (output) DOUBLE PRECISION 00085 * The value RESULT, It is currently limited to 1/ulp, to 00086 * avoid overflow. Errors are flagged by RESULT=10/ulp. 00087 * 00088 * ===================================================================== 00089 * 00090 * .. Parameters .. 00091 DOUBLE PRECISION ZERO, ONE 00092 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00093 * .. 00094 * .. Local Scalars .. 00095 DOUBLE PRECISION ABNORM, ULP, UNFL, WNORM 00096 * .. 00097 * .. Local Arrays .. 00098 DOUBLE PRECISION DUM( 1 ) 00099 * .. 00100 * .. External Functions .. 00101 DOUBLE PRECISION DLAMCH, DLANGE 00102 EXTERNAL DLAMCH, DLANGE 00103 * .. 00104 * .. External Subroutines .. 00105 EXTERNAL DGEMM, DLACPY 00106 * .. 00107 * .. Intrinsic Functions .. 00108 INTRINSIC DBLE, MAX, MIN 00109 * .. 00110 * .. Executable Statements .. 00111 * 00112 RESULT = ZERO 00113 IF( N.LE.0 ) 00114 $ RETURN 00115 * 00116 * Constants 00117 * 00118 UNFL = DLAMCH( 'Safe minimum' ) 00119 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 00120 * 00121 * compute the norm of (A,B) 00122 * 00123 CALL DLACPY( 'Full', N, N, A, LDA, WORK, N ) 00124 CALL DLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N ) 00125 ABNORM = MAX( DLANGE( '1', N, 2*N, WORK, N, DUM ), UNFL ) 00126 * 00127 * Compute W1 = A - U*S*V', and put in the array WORK(1:N*N) 00128 * 00129 CALL DLACPY( ' ', N, N, A, LDA, WORK, N ) 00130 CALL DGEMM( 'N', 'N', N, N, N, ONE, U, LDU, S, LDS, ZERO, 00131 $ WORK( N*N+1 ), N ) 00132 * 00133 CALL DGEMM( 'N', 'C', N, N, N, -ONE, WORK( N*N+1 ), N, V, LDV, 00134 $ ONE, WORK, N ) 00135 * 00136 * Compute W2 = B - U*T*V', and put in the workarray W(N*N+1:2*N*N) 00137 * 00138 CALL DLACPY( ' ', N, N, B, LDB, WORK( N*N+1 ), N ) 00139 CALL DGEMM( 'N', 'N', N, N, N, ONE, U, LDU, T, LDT, ZERO, 00140 $ WORK( 2*N*N+1 ), N ) 00141 * 00142 CALL DGEMM( 'N', 'C', N, N, N, -ONE, WORK( 2*N*N+1 ), N, V, LDV, 00143 $ ONE, WORK( N*N+1 ), N ) 00144 * 00145 * Compute norm(W)/ ( ulp*norm((A,B)) ) 00146 * 00147 WNORM = DLANGE( '1', N, 2*N, WORK, N, DUM ) 00148 * 00149 IF( ABNORM.GT.WNORM ) THEN 00150 RESULT = ( WNORM / ABNORM ) / ( 2*N*ULP ) 00151 ELSE 00152 IF( ABNORM.LT.ONE ) THEN 00153 RESULT = ( MIN( WNORM, 2*N*ABNORM ) / ABNORM ) / ( 2*N*ULP ) 00154 ELSE 00155 RESULT = MIN( WNORM / ABNORM, DBLE( 2*N ) ) / ( 2*N*ULP ) 00156 END IF 00157 END IF 00158 * 00159 RETURN 00160 * 00161 * End of DGET54 00162 * 00163 END
1.7.2