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
226