LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
shst01.f
Go to the documentation of this file.
1 *> \brief \b SHST01
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 SHST01( N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK,
12 * LWORK, RESULT )
13 *
14 * .. Scalar Arguments ..
15 * INTEGER IHI, ILO, LDA, LDH, LDQ, LWORK, N
16 * ..
17 * .. Array Arguments ..
18 * REAL A( LDA, * ), H( LDH, * ), Q( LDQ, * ),
19 * $ RESULT( 2 ), WORK( LWORK )
20 * ..
21 *
22 *
23 *> \par Purpose:
24 * =============
25 *>
26 *> \verbatim
27 *>
28 *> SHST01 tests the reduction of a general matrix A to upper Hessenberg
29 *> form: A = Q*H*Q'. Two test ratios are computed;
30 *>
31 *> RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
32 *> RESULT(2) = norm( I - Q'*Q ) / ( N * EPS )
33 *>
34 *> The matrix Q is assumed to be given explicitly as it would be
35 *> following SGEHRD + SORGHR.
36 *>
37 *> In this version, ILO and IHI are not used and are assumed to be 1 and
38 *> N, respectively.
39 *> \endverbatim
40 *
41 * Arguments:
42 * ==========
43 *
44 *> \param[in] N
45 *> \verbatim
46 *> N is INTEGER
47 *> The order of the matrix A. N >= 0.
48 *> \endverbatim
49 *>
50 *> \param[in] ILO
51 *> \verbatim
52 *> ILO is INTEGER
53 *> \endverbatim
54 *>
55 *> \param[in] IHI
56 *> \verbatim
57 *> IHI is INTEGER
58 *>
59 *> A is assumed to be upper triangular in rows and columns
60 *> 1:ILO-1 and IHI+1:N, so Q differs from the identity only in
61 *> rows and columns ILO+1:IHI.
62 *> \endverbatim
63 *>
64 *> \param[in] A
65 *> \verbatim
66 *> A is REAL array, dimension (LDA,N)
67 *> The original n by n matrix A.
68 *> \endverbatim
69 *>
70 *> \param[in] LDA
71 *> \verbatim
72 *> LDA is INTEGER
73 *> The leading dimension of the array A. LDA >= max(1,N).
74 *> \endverbatim
75 *>
76 *> \param[in] H
77 *> \verbatim
78 *> H is REAL array, dimension (LDH,N)
79 *> The upper Hessenberg matrix H from the reduction A = Q*H*Q'
80 *> as computed by SGEHRD. H is assumed to be zero below the
81 *> first subdiagonal.
82 *> \endverbatim
83 *>
84 *> \param[in] LDH
85 *> \verbatim
86 *> LDH is INTEGER
87 *> The leading dimension of the array H. LDH >= max(1,N).
88 *> \endverbatim
89 *>
90 *> \param[in] Q
91 *> \verbatim
92 *> Q is REAL array, dimension (LDQ,N)
93 *> The orthogonal matrix Q from the reduction A = Q*H*Q' as
94 *> computed by SGEHRD + SORGHR.
95 *> \endverbatim
96 *>
97 *> \param[in] LDQ
98 *> \verbatim
99 *> LDQ is INTEGER
100 *> The leading dimension of the array Q. LDQ >= max(1,N).
101 *> \endverbatim
102 *>
103 *> \param[out] WORK
104 *> \verbatim
105 *> WORK is REAL array, dimension (LWORK)
106 *> \endverbatim
107 *>
108 *> \param[in] LWORK
109 *> \verbatim
110 *> LWORK is INTEGER
111 *> The length of the array WORK. LWORK >= 2*N*N.
112 *> \endverbatim
113 *>
114 *> \param[out] RESULT
115 *> \verbatim
116 *> RESULT is REAL array, dimension (2)
117 *> RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
118 *> RESULT(2) = norm( I - Q'*Q ) / ( N * EPS )
119 *> \endverbatim
120 *
121 * Authors:
122 * ========
123 *
124 *> \author Univ. of Tennessee
125 *> \author Univ. of California Berkeley
126 *> \author Univ. of Colorado Denver
127 *> \author NAG Ltd.
128 *
129 *> \date November 2011
130 *
131 *> \ingroup single_eig
132 *
133 * =====================================================================
134  SUBROUTINE shst01( N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK,
135  $ lwork, result )
136 *
137 * -- LAPACK test routine (version 3.4.0) --
138 * -- LAPACK is a software package provided by Univ. of Tennessee, --
139 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
140 * November 2011
141 *
142 * .. Scalar Arguments ..
143  INTEGER ihi, ilo, lda, ldh, ldq, lwork, n
144 * ..
145 * .. Array Arguments ..
146  REAL a( lda, * ), h( ldh, * ), q( ldq, * ),
147  $ result( 2 ), work( lwork )
148 * ..
149 *
150 * =====================================================================
151 *
152 * .. Parameters ..
153  REAL one, zero
154  parameter( one = 1.0e+0, zero = 0.0e+0 )
155 * ..
156 * .. Local Scalars ..
157  INTEGER ldwork
158  REAL anorm, eps, ovfl, smlnum, unfl, wnorm
159 * ..
160 * .. External Functions ..
161  REAL slamch, slange
162  EXTERNAL slamch, slange
163 * ..
164 * .. External Subroutines ..
165  EXTERNAL sgemm, slabad, slacpy, sort01
166 * ..
167 * .. Intrinsic Functions ..
168  INTRINSIC max, min
169 * ..
170 * .. Executable Statements ..
171 *
172 * Quick return if possible
173 *
174  IF( n.LE.0 ) THEN
175  result( 1 ) = zero
176  result( 2 ) = zero
177  return
178  END IF
179 *
180  unfl = slamch( 'Safe minimum' )
181  eps = slamch( 'Precision' )
182  ovfl = one / unfl
183  CALL slabad( unfl, ovfl )
184  smlnum = unfl*n / eps
185 *
186 * Test 1: Compute norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
187 *
188 * Copy A to WORK
189 *
190  ldwork = max( 1, n )
191  CALL slacpy( ' ', n, n, a, lda, work, ldwork )
192 *
193 * Compute Q*H
194 *
195  CALL sgemm( 'No transpose', 'No transpose', n, n, n, one, q, ldq,
196  $ h, ldh, zero, work( ldwork*n+1 ), ldwork )
197 *
198 * Compute A - Q*H*Q'
199 *
200  CALL sgemm( 'No transpose', 'Transpose', n, n, n, -one,
201  $ work( ldwork*n+1 ), ldwork, q, ldq, one, work,
202  $ ldwork )
203 *
204  anorm = max( slange( '1', n, n, a, lda, work( ldwork*n+1 ) ),
205  $ unfl )
206  wnorm = slange( '1', n, n, work, ldwork, work( ldwork*n+1 ) )
207 *
208 * Note that RESULT(1) cannot overflow and is bounded by 1/(N*EPS)
209 *
210  result( 1 ) = min( wnorm, anorm ) / max( smlnum, anorm*eps ) / n
211 *
212 * Test 2: Compute norm( I - Q'*Q ) / ( N * EPS )
213 *
214  CALL sort01( 'Columns', n, n, q, ldq, work, lwork, result( 2 ) )
215 *
216  return
217 *
218 * End of SHST01
219 *
220  END