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