LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
zhst01.f
Go to the documentation of this file.
1 *> \brief \b ZHST01
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE ZHST01( N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK,
12 * LWORK, RWORK, RESULT )
13 *
14 * .. Scalar Arguments ..
15 * INTEGER IHI, ILO, LDA, LDH, LDQ, LWORK, N
16 * ..
17 * .. Array Arguments ..
18 * DOUBLE PRECISION RESULT( 2 ), RWORK( * )
19 * COMPLEX*16 A( LDA, * ), H( LDH, * ), Q( LDQ, * ),
20 * $ WORK( LWORK )
21 * ..
22 *
23 *
24 *> \par Purpose:
25 * =============
26 *>
27 *> \verbatim
28 *>
29 *> ZHST01 tests the reduction of a general matrix A to upper Hessenberg
30 *> form: A = Q*H*Q'. Two test ratios are computed;
31 *>
32 *> RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
33 *> RESULT(2) = norm( I - Q'*Q ) / ( N * EPS )
34 *>
35 *> The matrix Q is assumed to be given explicitly as it would be
36 *> following ZGEHRD + ZUNGHR.
37 *>
38 *> In this version, ILO and IHI are not used, but they could be used
39 *> to save some work if this is desired.
40 *> \endverbatim
41 *
42 * Arguments:
43 * ==========
44 *
45 *> \param[in] N
46 *> \verbatim
47 *> N is INTEGER
48 *> The order of the matrix A. N >= 0.
49 *> \endverbatim
50 *>
51 *> \param[in] ILO
52 *> \verbatim
53 *> ILO is INTEGER
54 *> \endverbatim
55 *>
56 *> \param[in] IHI
57 *> \verbatim
58 *> IHI is INTEGER
59 *>
60 *> A is assumed to be upper triangular in rows and columns
61 *> 1:ILO-1 and IHI+1:N, so Q differs from the identity only in
62 *> rows and columns ILO+1:IHI.
63 *> \endverbatim
64 *>
65 *> \param[in] A
66 *> \verbatim
67 *> A is COMPLEX*16 array, dimension (LDA,N)
68 *> The original n by n matrix A.
69 *> \endverbatim
70 *>
71 *> \param[in] LDA
72 *> \verbatim
73 *> LDA is INTEGER
74 *> The leading dimension of the array A. LDA >= max(1,N).
75 *> \endverbatim
76 *>
77 *> \param[in] H
78 *> \verbatim
79 *> H is COMPLEX*16 array, dimension (LDH,N)
80 *> The upper Hessenberg matrix H from the reduction A = Q*H*Q'
81 *> as computed by ZGEHRD. H is assumed to be zero below the
82 *> first subdiagonal.
83 *> \endverbatim
84 *>
85 *> \param[in] LDH
86 *> \verbatim
87 *> LDH is INTEGER
88 *> The leading dimension of the array H. LDH >= max(1,N).
89 *> \endverbatim
90 *>
91 *> \param[in] Q
92 *> \verbatim
93 *> Q is COMPLEX*16 array, dimension (LDQ,N)
94 *> The orthogonal matrix Q from the reduction A = Q*H*Q' as
95 *> computed by ZGEHRD + ZUNGHR.
96 *> \endverbatim
97 *>
98 *> \param[in] LDQ
99 *> \verbatim
100 *> LDQ is INTEGER
101 *> The leading dimension of the array Q. LDQ >= max(1,N).
102 *> \endverbatim
103 *>
104 *> \param[out] WORK
105 *> \verbatim
106 *> WORK is COMPLEX*16 array, dimension (LWORK)
107 *> \endverbatim
108 *>
109 *> \param[in] LWORK
110 *> \verbatim
111 *> LWORK is INTEGER
112 *> The length of the array WORK. LWORK >= 2*N*N.
113 *> \endverbatim
114 *>
115 *> \param[out] RWORK
116 *> \verbatim
117 *> RWORK is DOUBLE PRECISION array, dimension (N)
118 *> \endverbatim
119 *>
120 *> \param[out] RESULT
121 *> \verbatim
122 *> RESULT is DOUBLE PRECISION array, dimension (2)
123 *> RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
124 *> RESULT(2) = norm( I - Q'*Q ) / ( N * EPS )
125 *> \endverbatim
126 *
127 * Authors:
128 * ========
129 *
130 *> \author Univ. of Tennessee
131 *> \author Univ. of California Berkeley
132 *> \author Univ. of Colorado Denver
133 *> \author NAG Ltd.
134 *
135 *> \ingroup complex16_eig
136 *
137 * =====================================================================
138  SUBROUTINE zhst01( N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK,
139  $ LWORK, RWORK, RESULT )
140 *
141 * -- LAPACK test routine --
142 * -- LAPACK is a software package provided by Univ. of Tennessee, --
143 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
144 *
145 * .. Scalar Arguments ..
146  INTEGER IHI, ILO, LDA, LDH, LDQ, LWORK, N
147 * ..
148 * .. Array Arguments ..
149  DOUBLE PRECISION RESULT( 2 ), RWORK( * )
150  COMPLEX*16 A( LDA, * ), H( LDH, * ), Q( LDQ, * ),
151  $ work( lwork )
152 * ..
153 *
154 * =====================================================================
155 *
156 * .. Parameters ..
157  DOUBLE PRECISION ONE, ZERO
158  parameter( one = 1.0d+0, zero = 0.0d+0 )
159 * ..
160 * .. Local Scalars ..
161  INTEGER LDWORK
162  DOUBLE PRECISION ANORM, EPS, OVFL, SMLNUM, UNFL, WNORM
163 * ..
164 * .. External Functions ..
165  DOUBLE PRECISION DLAMCH, ZLANGE
166  EXTERNAL dlamch, zlange
167 * ..
168 * .. External Subroutines ..
169  EXTERNAL dlabad, zgemm, zlacpy, zunt01
170 * ..
171 * .. Intrinsic Functions ..
172  INTRINSIC dcmplx, max, min
173 * ..
174 * .. Executable Statements ..
175 *
176 * Quick return if possible
177 *
178  IF( n.LE.0 ) THEN
179  result( 1 ) = zero
180  result( 2 ) = zero
181  RETURN
182  END IF
183 *
184  unfl = dlamch( 'Safe minimum' )
185  eps = dlamch( 'Precision' )
186  ovfl = one / unfl
187  CALL dlabad( unfl, ovfl )
188  smlnum = unfl*n / eps
189 *
190 * Test 1: Compute norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
191 *
192 * Copy A to WORK
193 *
194  ldwork = max( 1, n )
195  CALL zlacpy( ' ', n, n, a, lda, work, ldwork )
196 *
197 * Compute Q*H
198 *
199  CALL zgemm( 'No transpose', 'No transpose', n, n, n,
200  $ dcmplx( one ), q, ldq, h, ldh, dcmplx( zero ),
201  $ work( ldwork*n+1 ), ldwork )
202 *
203 * Compute A - Q*H*Q'
204 *
205  CALL zgemm( 'No transpose', 'Conjugate transpose', n, n, n,
206  $ dcmplx( -one ), work( ldwork*n+1 ), ldwork, q, ldq,
207  $ dcmplx( one ), work, ldwork )
208 *
209  anorm = max( zlange( '1', n, n, a, lda, rwork ), unfl )
210  wnorm = zlange( '1', n, n, work, ldwork, rwork )
211 *
212 * Note that RESULT(1) cannot overflow and is bounded by 1/(N*EPS)
213 *
214  result( 1 ) = min( wnorm, anorm ) / max( smlnum, anorm*eps ) / n
215 *
216 * Test 2: Compute norm( I - Q'*Q ) / ( N * EPS )
217 *
218  CALL zunt01( 'Columns', n, n, q, ldq, work, lwork, rwork,
219  $ result( 2 ) )
220 *
221  RETURN
222 *
223 * End of ZHST01
224 *
225  END
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine zhst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RWORK, RESULT)
ZHST01
Definition: zhst01.f:140
subroutine zunt01(ROWCOL, M, N, U, LDU, WORK, LWORK, RWORK, RESID)
ZUNT01
Definition: zunt01.f:126
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103