LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \date November 2011
136 *
137 *> \ingroup complex16_eig
138 *
139 * =====================================================================
140  SUBROUTINE zhst01( N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK,
141  $ lwork, rwork, result )
142 *
143 * -- LAPACK test routine (version 3.4.0) --
144 * -- LAPACK is a software package provided by Univ. of Tennessee, --
145 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
146 * November 2011
147 *
148 * .. Scalar Arguments ..
149  INTEGER ihi, ilo, lda, ldh, ldq, lwork, n
150 * ..
151 * .. Array Arguments ..
152  DOUBLE PRECISION result( 2 ), rwork( * )
153  COMPLEX*16 a( lda, * ), h( ldh, * ), q( ldq, * ),
154  $ work( lwork )
155 * ..
156 *
157 * =====================================================================
158 *
159 * .. Parameters ..
160  DOUBLE PRECISION one, zero
161  parameter( one = 1.0d+0, zero = 0.0d+0 )
162 * ..
163 * .. Local Scalars ..
164  INTEGER ldwork
165  DOUBLE PRECISION anorm, eps, ovfl, smlnum, unfl, wnorm
166 * ..
167 * .. External Functions ..
168  DOUBLE PRECISION dlamch, zlange
169  EXTERNAL dlamch, zlange
170 * ..
171 * .. External Subroutines ..
172  EXTERNAL dlabad, zgemm, zlacpy, zunt01
173 * ..
174 * .. Intrinsic Functions ..
175  INTRINSIC dcmplx, max, min
176 * ..
177 * .. Executable Statements ..
178 *
179 * Quick return if possible
180 *
181  IF( n.LE.0 ) THEN
182  result( 1 ) = zero
183  result( 2 ) = zero
184  return
185  END IF
186 *
187  unfl = dlamch( 'Safe minimum' )
188  eps = dlamch( 'Precision' )
189  ovfl = one / unfl
190  CALL dlabad( unfl, ovfl )
191  smlnum = unfl*n / eps
192 *
193 * Test 1: Compute norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
194 *
195 * Copy A to WORK
196 *
197  ldwork = max( 1, n )
198  CALL zlacpy( ' ', n, n, a, lda, work, ldwork )
199 *
200 * Compute Q*H
201 *
202  CALL zgemm( 'No transpose', 'No transpose', n, n, n,
203  $ dcmplx( one ), q, ldq, h, ldh, dcmplx( zero ),
204  $ work( ldwork*n+1 ), ldwork )
205 *
206 * Compute A - Q*H*Q'
207 *
208  CALL zgemm( 'No transpose', 'Conjugate transpose', n, n, n,
209  $ dcmplx( -one ), work( ldwork*n+1 ), ldwork, q, ldq,
210  $ dcmplx( one ), work, ldwork )
211 *
212  anorm = max( zlange( '1', n, n, a, lda, rwork ), unfl )
213  wnorm = zlange( '1', n, n, work, ldwork, rwork )
214 *
215 * Note that RESULT(1) cannot overflow and is bounded by 1/(N*EPS)
216 *
217  result( 1 ) = min( wnorm, anorm ) / max( smlnum, anorm*eps ) / n
218 *
219 * Test 2: Compute norm( I - Q'*Q ) / ( N * EPS )
220 *
221  CALL zunt01( 'Columns', n, n, q, ldq, work, lwork, rwork,
222  $ result( 2 ) )
223 *
224  return
225 *
226 * End of ZHST01
227 *
228  END