00001 SUBROUTINE DHST01( N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, 00002 $ LWORK, 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 IHI, ILO, LDA, LDH, LDQ, LWORK, N 00010 * .. 00011 * .. Array Arguments .. 00012 DOUBLE PRECISION A( LDA, * ), H( LDH, * ), Q( LDQ, * ), 00013 $ RESULT( 2 ), WORK( LWORK ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * DHST01 tests the reduction of a general matrix A to upper Hessenberg 00020 * form: A = Q*H*Q'. Two test ratios are computed; 00021 * 00022 * RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS ) 00023 * RESULT(2) = norm( I - Q'*Q ) / ( N * EPS ) 00024 * 00025 * The matrix Q is assumed to be given explicitly as it would be 00026 * following DGEHRD + DORGHR. 00027 * 00028 * In this version, ILO and IHI are not used and are assumed to be 1 and 00029 * N, respectively. 00030 * 00031 * Arguments 00032 * ========= 00033 * 00034 * N (input) INTEGER 00035 * The order of the matrix A. N >= 0. 00036 * 00037 * ILO (input) INTEGER 00038 * IHI (input) INTEGER 00039 * A is assumed to be upper triangular in rows and columns 00040 * 1:ILO-1 and IHI+1:N, so Q differs from the identity only in 00041 * rows and columns ILO+1:IHI. 00042 * 00043 * A (input) DOUBLE PRECISION array, dimension (LDA,N) 00044 * The original n by n matrix A. 00045 * 00046 * LDA (input) INTEGER 00047 * The leading dimension of the array A. LDA >= max(1,N). 00048 * 00049 * H (input) DOUBLE PRECISION array, dimension (LDH,N) 00050 * The upper Hessenberg matrix H from the reduction A = Q*H*Q' 00051 * as computed by DGEHRD. H is assumed to be zero below the 00052 * first subdiagonal. 00053 * 00054 * LDH (input) INTEGER 00055 * The leading dimension of the array H. LDH >= max(1,N). 00056 * 00057 * Q (input) DOUBLE PRECISION array, dimension (LDQ,N) 00058 * The orthogonal matrix Q from the reduction A = Q*H*Q' as 00059 * computed by DGEHRD + DORGHR. 00060 * 00061 * LDQ (input) INTEGER 00062 * The leading dimension of the array Q. LDQ >= max(1,N). 00063 * 00064 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) 00065 * 00066 * LWORK (input) INTEGER 00067 * The length of the array WORK. LWORK >= 2*N*N. 00068 * 00069 * RESULT (output) DOUBLE PRECISION array, dimension (2) 00070 * RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS ) 00071 * RESULT(2) = norm( I - Q'*Q ) / ( N * EPS ) 00072 * 00073 * ===================================================================== 00074 * 00075 * .. Parameters .. 00076 DOUBLE PRECISION ONE, ZERO 00077 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00078 * .. 00079 * .. Local Scalars .. 00080 INTEGER LDWORK 00081 DOUBLE PRECISION ANORM, EPS, OVFL, SMLNUM, UNFL, WNORM 00082 * .. 00083 * .. External Functions .. 00084 DOUBLE PRECISION DLAMCH, DLANGE 00085 EXTERNAL DLAMCH, DLANGE 00086 * .. 00087 * .. External Subroutines .. 00088 EXTERNAL DGEMM, DLABAD, DLACPY, DORT01 00089 * .. 00090 * .. Intrinsic Functions .. 00091 INTRINSIC MAX, MIN 00092 * .. 00093 * .. Executable Statements .. 00094 * 00095 * Quick return if possible 00096 * 00097 IF( N.LE.0 ) THEN 00098 RESULT( 1 ) = ZERO 00099 RESULT( 2 ) = ZERO 00100 RETURN 00101 END IF 00102 * 00103 UNFL = DLAMCH( 'Safe minimum' ) 00104 EPS = DLAMCH( 'Precision' ) 00105 OVFL = ONE / UNFL 00106 CALL DLABAD( UNFL, OVFL ) 00107 SMLNUM = UNFL*N / EPS 00108 * 00109 * Test 1: Compute norm( A - Q*H*Q' ) / ( norm(A) * N * EPS ) 00110 * 00111 * Copy A to WORK 00112 * 00113 LDWORK = MAX( 1, N ) 00114 CALL DLACPY( ' ', N, N, A, LDA, WORK, LDWORK ) 00115 * 00116 * Compute Q*H 00117 * 00118 CALL DGEMM( 'No transpose', 'No transpose', N, N, N, ONE, Q, LDQ, 00119 $ H, LDH, ZERO, WORK( LDWORK*N+1 ), LDWORK ) 00120 * 00121 * Compute A - Q*H*Q' 00122 * 00123 CALL DGEMM( 'No transpose', 'Transpose', N, N, N, -ONE, 00124 $ WORK( LDWORK*N+1 ), LDWORK, Q, LDQ, ONE, WORK, 00125 $ LDWORK ) 00126 * 00127 ANORM = MAX( DLANGE( '1', N, N, A, LDA, WORK( LDWORK*N+1 ) ), 00128 $ UNFL ) 00129 WNORM = DLANGE( '1', N, N, WORK, LDWORK, WORK( LDWORK*N+1 ) ) 00130 * 00131 * Note that RESULT(1) cannot overflow and is bounded by 1/(N*EPS) 00132 * 00133 RESULT( 1 ) = MIN( WNORM, ANORM ) / MAX( SMLNUM, ANORM*EPS ) / N 00134 * 00135 * Test 2: Compute norm( I - Q'*Q ) / ( N * EPS ) 00136 * 00137 CALL DORT01( 'Columns', N, N, Q, LDQ, WORK, LWORK, RESULT( 2 ) ) 00138 * 00139 RETURN 00140 * 00141 * End of DHST01 00142 * 00143 END
1.7.2