LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ sstt21()

subroutine sstt21 ( integer  N,
integer  KBAND,
real, dimension( * )  AD,
real, dimension( * )  AE,
real, dimension( * )  SD,
real, dimension( * )  SE,
real, dimension( ldu, * )  U,
integer  LDU,
real, dimension( * )  WORK,
real, dimension( 2 )  RESULT 
)

SSTT21

Purpose:
 SSTT21 checks a decomposition of the form

    A = U S U'

 where ' means transpose, A is symmetric tridiagonal, U is orthogonal,
 and S is diagonal (if KBAND=0) or symmetric tridiagonal (if KBAND=1).
 Two tests are performed:

    RESULT(1) = | A - U S U' | / ( |A| n ulp )

    RESULT(2) = | I - UU' | / ( n ulp )
Parameters
[in]N
          N is INTEGER
          The size of the matrix.  If it is zero, SSTT21 does nothing.
          It must be at least zero.
[in]KBAND
          KBAND is INTEGER
          The bandwidth of the matrix S.  It may only be zero or one.
          If zero, then S is diagonal, and SE is not referenced.  If
          one, then S is symmetric tri-diagonal.
[in]AD
          AD is REAL array, dimension (N)
          The diagonal of the original (unfactored) matrix A.  A is
          assumed to be symmetric tridiagonal.
[in]AE
          AE is REAL array, dimension (N-1)
          The off-diagonal of the original (unfactored) matrix A.  A
          is assumed to be symmetric tridiagonal.  AE(1) is the (1,2)
          and (2,1) element, AE(2) is the (2,3) and (3,2) element, etc.
[in]SD
          SD is REAL array, dimension (N)
          The diagonal of the (symmetric tri-) diagonal matrix S.
[in]SE
          SE is REAL array, dimension (N-1)
          The off-diagonal of the (symmetric tri-) diagonal matrix S.
          Not referenced if KBSND=0.  If KBAND=1, then AE(1) is the
          (1,2) and (2,1) element, SE(2) is the (2,3) and (3,2)
          element, etc.
[in]U
          U is REAL array, dimension (LDU, N)
          The orthogonal matrix in the decomposition.
[in]LDU
          LDU is INTEGER
          The leading dimension of U.  LDU must be at least N.
[out]WORK
          WORK is REAL array, dimension (N*(N+1))
[out]RESULT
          RESULT is REAL array, dimension (2)
          The values computed by the two tests described above.  The
          values are currently limited to 1/ulp, to avoid overflow.
          RESULT(1) is always modified.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 129 of file sstt21.f.

129 *
130 * -- LAPACK test routine (version 3.7.0) --
131 * -- LAPACK is a software package provided by Univ. of Tennessee, --
132 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
133 * December 2016
134 *
135 * .. Scalar Arguments ..
136  INTEGER kband, ldu, n
137 * ..
138 * .. Array Arguments ..
139  REAL ad( * ), ae( * ), result( 2 ), sd( * ),
140  $ se( * ), u( ldu, * ), work( * )
141 * ..
142 *
143 * =====================================================================
144 *
145 * .. Parameters ..
146  REAL zero, one
147  parameter( zero = 0.0e0, one = 1.0e0 )
148 * ..
149 * .. Local Scalars ..
150  INTEGER j
151  REAL anorm, temp1, temp2, ulp, unfl, wnorm
152 * ..
153 * .. External Functions ..
154  REAL slamch, slange, slansy
155  EXTERNAL slamch, slange, slansy
156 * ..
157 * .. External Subroutines ..
158  EXTERNAL sgemm, slaset, ssyr, ssyr2
159 * ..
160 * .. Intrinsic Functions ..
161  INTRINSIC abs, max, min, real
162 * ..
163 * .. Executable Statements ..
164 *
165 * 1) Constants
166 *
167  result( 1 ) = zero
168  result( 2 ) = zero
169  IF( n.LE.0 )
170  $ RETURN
171 *
172  unfl = slamch( 'Safe minimum' )
173  ulp = slamch( 'Precision' )
174 *
175 * Do Test 1
176 *
177 * Copy A & Compute its 1-Norm:
178 *
179  CALL slaset( 'Full', n, n, zero, zero, work, n )
180 *
181  anorm = zero
182  temp1 = zero
183 *
184  DO 10 j = 1, n - 1
185  work( ( n+1 )*( j-1 )+1 ) = ad( j )
186  work( ( n+1 )*( j-1 )+2 ) = ae( j )
187  temp2 = abs( ae( j ) )
188  anorm = max( anorm, abs( ad( j ) )+temp1+temp2 )
189  temp1 = temp2
190  10 CONTINUE
191 *
192  work( n**2 ) = ad( n )
193  anorm = max( anorm, abs( ad( n ) )+temp1, unfl )
194 *
195 * Norm of A - USU'
196 *
197  DO 20 j = 1, n
198  CALL ssyr( 'L', n, -sd( j ), u( 1, j ), 1, work, n )
199  20 CONTINUE
200 *
201  IF( n.GT.1 .AND. kband.EQ.1 ) THEN
202  DO 30 j = 1, n - 1
203  CALL ssyr2( 'L', n, -se( j ), u( 1, j ), 1, u( 1, j+1 ), 1,
204  $ work, n )
205  30 CONTINUE
206  END IF
207 *
208  wnorm = slansy( '1', 'L', n, work, n, work( n**2+1 ) )
209 *
210  IF( anorm.GT.wnorm ) THEN
211  result( 1 ) = ( wnorm / anorm ) / ( n*ulp )
212  ELSE
213  IF( anorm.LT.one ) THEN
214  result( 1 ) = ( min( wnorm, n*anorm ) / anorm ) / ( n*ulp )
215  ELSE
216  result( 1 ) = min( wnorm / anorm, REAL( N ) ) / ( n*ulp )
217  END IF
218  END IF
219 *
220 * Do Test 2
221 *
222 * Compute UU' - I
223 *
224  CALL sgemm( 'N', 'C', n, n, n, one, u, ldu, u, ldu, zero, work,
225  $ n )
226 *
227  DO 40 j = 1, n
228  work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - one
229  40 CONTINUE
230 *
231  result( 2 ) = min( REAL( N ), slange( '1', n, n, work, n,
232  $ work( n**2+1 ) ) ) / ( n*ulp )
233 *
234  RETURN
235 *
236 * End of SSTT21
237 *
subroutine ssyr2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
SSYR2
Definition: ssyr2.f:149
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
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
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:112
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine ssyr(UPLO, N, ALPHA, X, INCX, A, LDA)
SSYR
Definition: ssyr.f:134
real function slansy(NORM, UPLO, N, A, LDA, WORK)
SLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric matrix.
Definition: slansy.f:124
Here is the call graph for this function:
Here is the caller graph for this function: