LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ shst01()

subroutine shst01 ( integer  N,
integer  ILO,
integer  IHI,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( ldh, * )  H,
integer  LDH,
real, dimension( ldq, * )  Q,
integer  LDQ,
real, dimension( lwork )  WORK,
integer  LWORK,
real, dimension( 2 )  RESULT 
)

SHST01

Purpose:
 SHST01 tests the reduction of a general matrix A to upper Hessenberg
 form:  A = Q*H*Q'.  Two test ratios are computed;

 RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
 RESULT(2) = norm( I - Q'*Q ) / ( N * EPS )

 The matrix Q is assumed to be given explicitly as it would be
 following SGEHRD + SORGHR.

 In this version, ILO and IHI are not used and are assumed to be 1 and
 N, respectively.
Parameters
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]ILO
          ILO is INTEGER
[in]IHI
          IHI is INTEGER

          A is assumed to be upper triangular in rows and columns
          1:ILO-1 and IHI+1:N, so Q differs from the identity only in
          rows and columns ILO+1:IHI.
[in]A
          A is REAL array, dimension (LDA,N)
          The original n by n matrix A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]H
          H is REAL array, dimension (LDH,N)
          The upper Hessenberg matrix H from the reduction A = Q*H*Q'
          as computed by SGEHRD.  H is assumed to be zero below the
          first subdiagonal.
[in]LDH
          LDH is INTEGER
          The leading dimension of the array H.  LDH >= max(1,N).
[in]Q
          Q is REAL array, dimension (LDQ,N)
          The orthogonal matrix Q from the reduction A = Q*H*Q' as
          computed by SGEHRD + SORGHR.
[in]LDQ
          LDQ is INTEGER
          The leading dimension of the array Q.  LDQ >= max(1,N).
[out]WORK
          WORK is REAL array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The length of the array WORK.  LWORK >= 2*N*N.
[out]RESULT
          RESULT is REAL array, dimension (2)
          RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
          RESULT(2) = norm( I - Q'*Q ) / ( N * EPS )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 136 of file shst01.f.

136 *
137 * -- LAPACK test routine (version 3.7.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 * December 2016
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 *
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
subroutine sort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
SORT01
Definition: sort01.f:118
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: slange.f:116
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
Here is the call graph for this function:
Here is the caller graph for this function: