LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> \ingroup single_eig
130*
131* =====================================================================
132 SUBROUTINE shst01( N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK,
133 $ LWORK, RESULT )
134*
135* -- LAPACK test routine --
136* -- LAPACK is a software package provided by Univ. of Tennessee, --
137* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138*
139* .. Scalar Arguments ..
140 INTEGER IHI, ILO, LDA, LDH, LDQ, LWORK, N
141* ..
142* .. Array Arguments ..
143 REAL A( LDA, * ), H( LDH, * ), Q( LDQ, * ),
144 $ result( 2 ), work( lwork )
145* ..
146*
147* =====================================================================
148*
149* .. Parameters ..
150 REAL ONE, ZERO
151 parameter( one = 1.0e+0, zero = 0.0e+0 )
152* ..
153* .. Local Scalars ..
154 INTEGER LDWORK
155 REAL ANORM, EPS, OVFL, SMLNUM, UNFL, WNORM
156* ..
157* .. External Functions ..
158 REAL SLAMCH, SLANGE
159 EXTERNAL slamch, slange
160* ..
161* .. External Subroutines ..
162 EXTERNAL sgemm, slacpy, sort01
163* ..
164* .. Intrinsic Functions ..
165 INTRINSIC max, min
166* ..
167* .. Executable Statements ..
168*
169* Quick return if possible
170*
171 IF( n.LE.0 ) THEN
172 result( 1 ) = zero
173 result( 2 ) = zero
174 RETURN
175 END IF
176*
177 unfl = slamch( 'Safe minimum' )
178 eps = slamch( 'Precision' )
179 ovfl = one / unfl
180 smlnum = unfl*n / eps
181*
182* Test 1: Compute norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
183*
184* Copy A to WORK
185*
186 ldwork = max( 1, n )
187 CALL slacpy( ' ', n, n, a, lda, work, ldwork )
188*
189* Compute Q*H
190*
191 CALL sgemm( 'No transpose', 'No transpose', n, n, n, one, q, ldq,
192 $ h, ldh, zero, work( ldwork*n+1 ), ldwork )
193*
194* Compute A - Q*H*Q'
195*
196 CALL sgemm( 'No transpose', 'Transpose', n, n, n, -one,
197 $ work( ldwork*n+1 ), ldwork, q, ldq, one, work,
198 $ ldwork )
199*
200 anorm = max( slange( '1', n, n, a, lda, work( ldwork*n+1 ) ),
201 $ unfl )
202 wnorm = slange( '1', n, n, work, ldwork, work( ldwork*n+1 ) )
203*
204* Note that RESULT(1) cannot overflow and is bounded by 1/(N*EPS)
205*
206 result( 1 ) = min( wnorm, anorm ) / max( smlnum, anorm*eps ) / n
207*
208* Test 2: Compute norm( I - Q'*Q ) / ( N * EPS )
209*
210 CALL sort01( 'Columns', n, n, q, ldq, work, lwork, result( 2 ) )
211*
212 RETURN
213*
214* End of SHST01
215*
216 END
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
subroutine shst01(n, ilo, ihi, a, lda, h, ldh, q, ldq, work, lwork, result)
SHST01
Definition shst01.f:134
subroutine sort01(rowcol, m, n, u, ldu, work, lwork, resid)
SORT01
Definition sort01.f:116