1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 |
SUBROUTINE DHST01( N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK,
$ LWORK, RESULT ) * * -- LAPACK test routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER IHI, ILO, LDA, LDH, LDQ, LWORK, N * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), H( LDH, * ), Q( LDQ, * ), $ RESULT( 2 ), WORK( LWORK ) * .. * * Purpose * ======= * * DHST01 tests the reduction of a general matrix A to upper Hessenberg * form: A = Q*H*Q'. Two test ratios are computed; * * RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS ) * RESULT(2) = norm( I - Q'*Q ) / ( N * EPS ) * * The matrix Q is assumed to be given explicitly as it would be * following DGEHRD + DORGHR. * * In this version, ILO and IHI are not used and are assumed to be 1 and * N, respectively. * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. N >= 0. * * ILO (input) INTEGER * IHI (input) INTEGER * A is assumed to be upper triangular in rows and columns * 1:ILO-1 and IHI+1:N, so Q differs from the identity only in * rows and columns ILO+1:IHI. * * A (input) DOUBLE PRECISION array, dimension (LDA,N) * The original n by n matrix A. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * H (input) DOUBLE PRECISION array, dimension (LDH,N) * The upper Hessenberg matrix H from the reduction A = Q*H*Q' * as computed by DGEHRD. H is assumed to be zero below the * first subdiagonal. * * LDH (input) INTEGER * The leading dimension of the array H. LDH >= max(1,N). * * Q (input) DOUBLE PRECISION array, dimension (LDQ,N) * The orthogonal matrix Q from the reduction A = Q*H*Q' as * computed by DGEHRD + DORGHR. * * LDQ (input) INTEGER * The leading dimension of the array Q. LDQ >= max(1,N). * * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) * * LWORK (input) INTEGER * The length of the array WORK. LWORK >= 2*N*N. * * RESULT (output) DOUBLE PRECISION array, dimension (2) * RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS ) * RESULT(2) = norm( I - Q'*Q ) / ( N * EPS ) * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. INTEGER LDWORK DOUBLE PRECISION ANORM, EPS, OVFL, SMLNUM, UNFL, WNORM * .. * .. External Functions .. DOUBLE PRECISION DLAMCH, DLANGE EXTERNAL DLAMCH, DLANGE * .. * .. External Subroutines .. EXTERNAL DGEMM, DLABAD, DLACPY, DORT01 * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Quick return if possible * IF( N.LE.0 ) THEN RESULT( 1 ) = ZERO RESULT( 2 ) = ZERO RETURN END IF * UNFL = DLAMCH( 'Safe minimum' ) EPS = DLAMCH( 'Precision' ) OVFL = ONE / UNFL CALL DLABAD( UNFL, OVFL ) SMLNUM = UNFL*N / EPS * * Test 1: Compute norm( A - Q*H*Q' ) / ( norm(A) * N * EPS ) * * Copy A to WORK * LDWORK = MAX( 1, N ) CALL DLACPY( ' ', N, N, A, LDA, WORK, LDWORK ) * * Compute Q*H * CALL DGEMM( 'No transpose', 'No transpose', N, N, N, ONE, Q, LDQ, $ H, LDH, ZERO, WORK( LDWORK*N+1 ), LDWORK ) * * Compute A - Q*H*Q' * CALL DGEMM( 'No transpose', 'Transpose', N, N, N, -ONE, $ WORK( LDWORK*N+1 ), LDWORK, Q, LDQ, ONE, WORK, $ LDWORK ) * ANORM = MAX( DLANGE( '1', N, N, A, LDA, WORK( LDWORK*N+1 ) ), $ UNFL ) WNORM = DLANGE( '1', N, N, WORK, LDWORK, WORK( LDWORK*N+1 ) ) * * Note that RESULT(1) cannot overflow and is bounded by 1/(N*EPS) * RESULT( 1 ) = MIN( WNORM, ANORM ) / MAX( SMLNUM, ANORM*EPS ) / N * * Test 2: Compute norm( I - Q'*Q ) / ( N * EPS ) * CALL DORT01( 'Columns', N, N, Q, LDQ, WORK, LWORK, RESULT( 2 ) ) * RETURN * * End of DHST01 * END |