LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine ssyt21 ( integer  ITYPE,
character  UPLO,
integer  N,
integer  KBAND,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( * )  D,
real, dimension( * )  E,
real, dimension( ldu, * )  U,
integer  LDU,
real, dimension( ldv, * )  V,
integer  LDV,
real, dimension( * )  TAU,
real, dimension( * )  WORK,
real, dimension( 2 )  RESULT 
)

SSYT21

Purpose: 
 SSYT21 generally checks a decomposition of the form

    A = U S U'

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

 If ITYPE=1, then U is represented as a dense matrix; otherwise U is
 expressed as a product of Householder transformations, whose vectors
 are stored in the array "V" and whose scaling constants are in "TAU".
 We shall use the letter "V" to refer to the product of Householder
 transformations (which should be equal to U).

 Specifically, if ITYPE=1, then:

    RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC>    RESULT(2) = | I - UU' | / ( n ulp )

 If ITYPE=2, then:

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

 If ITYPE=3, then:

    RESULT(1) = | I - VU' | / ( n ulp )

 For ITYPE > 1, the transformation U is expressed as a product
 V = H(1)...H(n-2),  where H(j) = I  -  tau(j) v(j) v(j)' and each
 vector v(j) has its first j elements 0 and the remaining n-j elements
 stored in V(j+1:n,j).
Parameters
[in]ITYPE
          ITYPE is INTEGER
          Specifies the type of tests to be performed.
          1: U expressed as a dense orthogonal matrix:
             RESULT(1) = | A - U S U' | / ( |A| n ulp )   *andC>             RESULT(2) = | I - UU' | / ( n ulp )

          2: U expressed as a product V of Housholder transformations:
             RESULT(1) = | A - V S V' | / ( |A| n ulp )

          3: U expressed both as a dense orthogonal matrix and
             as a product of Housholder transformations:
             RESULT(1) = | I - VU' | / ( n ulp )
[in]UPLO
          UPLO is CHARACTER
          If UPLO='U', the upper triangle of A and V will be used and
          the (strictly) lower triangle will not be referenced.
          If UPLO='L', the lower triangle of A and V will be used and
          the (strictly) upper triangle will not be referenced.
[in]N
          N is INTEGER
          The size of the matrix.  If it is zero, SSYT21 does nothing.
          It must be at least zero.
[in]KBAND
          KBAND is INTEGER
          The bandwidth of the matrix.  It may only be zero or one.
          If zero, then S is diagonal, and E is not referenced.  If
          one, then S is symmetric tri-diagonal.
[in]A
          A is REAL array, dimension (LDA, N)
          The original (unfactored) matrix.  It is assumed to be
          symmetric, and only the upper (UPLO='U') or only the lower
          (UPLO='L') will be referenced.
[in]LDA
          LDA is INTEGER
          The leading dimension of A.  It must be at least 1
          and at least N.
[in]D
          D is REAL array, dimension (N)
          The diagonal of the (symmetric tri-) diagonal matrix.
[in]E
          E is REAL array, dimension (N-1)
          The off-diagonal of the (symmetric tri-) diagonal matrix.
          E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and
          (3,2) element, etc.
          Not referenced if KBAND=0.
[in]U
          U is REAL array, dimension (LDU, N)
          If ITYPE=1 or 3, this contains the orthogonal matrix in
          the decomposition, expressed as a dense matrix.  If ITYPE=2,
          then it is not referenced.
[in]LDU
          LDU is INTEGER
          The leading dimension of U.  LDU must be at least N and
          at least 1.
[in]V
          V is REAL array, dimension (LDV, N)
          If ITYPE=2 or 3, the columns of this array contain the
          Householder vectors used to describe the orthogonal matrix
          in the decomposition.  If UPLO='L', then the vectors are in
          the lower triangle, if UPLO='U', then in the upper
          triangle.
          *NOTE* If ITYPE=2 or 3, V is modified and restored.  The
          subdiagonal (if UPLO='L') or the superdiagonal (if UPLO='U')
          is set to one, and later reset to its original value, during
          the course of the calculation.
          If ITYPE=1, then it is neither referenced nor modified.
[in]LDV
          LDV is INTEGER
          The leading dimension of V.  LDV must be at least N and
          at least 1.
[in]TAU
          TAU is REAL array, dimension (N)
          If ITYPE >= 2, then TAU(j) is the scalar factor of
          v(j) v(j)' in the Householder transformation H(j) of
          the product  U = H(1)...H(n-2)
          If ITYPE < 2, then TAU is not referenced.
[out]WORK
          WORK is REAL array, dimension (2*N**2)
[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.  RESULT(2) is modified only
          if ITYPE=1.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 207 of file ssyt21.f.

207 *
208 * -- LAPACK test routine (version 3.4.0) --
209 * -- LAPACK is a software package provided by Univ. of Tennessee, --
210 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
211 * November 2011
212 *
213 * .. Scalar Arguments ..
214  CHARACTER uplo
215  INTEGER itype, kband, lda, ldu, ldv, n
216 * ..
217 * .. Array Arguments ..
218  REAL a( lda, * ), d( * ), e( * ), result( 2 ),
219  $ tau( * ), u( ldu, * ), v( ldv, * ), work( * )
220 * ..
221 *
222 * =====================================================================
223 *
224 * .. Parameters ..
225  REAL zero, one, ten
226  parameter ( zero = 0.0e0, one = 1.0e0, ten = 10.0e0 )
227 * ..
228 * .. Local Scalars ..
229  LOGICAL lower
230  CHARACTER cuplo
231  INTEGER iinfo, j, jcol, jr, jrow
232  REAL anorm, ulp, unfl, vsave, wnorm
233 * ..
234 * .. External Functions ..
235  LOGICAL lsame
236  REAL slamch, slange, slansy
237  EXTERNAL lsame, slamch, slange, slansy
238 * ..
239 * .. External Subroutines ..
240  EXTERNAL sgemm, slacpy, slarfy, slaset, sorm2l, sorm2r,
241  $ ssyr, ssyr2
242 * ..
243 * .. Intrinsic Functions ..
244  INTRINSIC max, min, real
245 * ..
246 * .. Executable Statements ..
247 *
248  result( 1 ) = zero
249  IF( itype.EQ.1 )
250  $ result( 2 ) = zero
251  IF( n.LE.0 )
252  $ RETURN
253 *
254  IF( lsame( uplo, 'U' ) ) THEN
255  lower = .false.
256  cuplo = 'U'
257  ELSE
258  lower = .true.
259  cuplo = 'L'
260  END IF
261 *
262  unfl = slamch( 'Safe minimum' )
263  ulp = slamch( 'Epsilon' )*slamch( 'Base' )
264 *
265 * Some Error Checks
266 *
267  IF( itype.LT.1 .OR. itype.GT.3 ) THEN
268  result( 1 ) = ten / ulp
269  RETURN
270  END IF
271 *
272 * Do Test 1
273 *
274 * Norm of A:
275 *
276  IF( itype.EQ.3 ) THEN
277  anorm = one
278  ELSE
279  anorm = max( slansy( '1', cuplo, n, a, lda, work ), unfl )
280  END IF
281 *
282 * Compute error matrix:
283 *
284  IF( itype.EQ.1 ) THEN
285 *
286 * ITYPE=1: error = A - U S U'
287 *
288  CALL slaset( 'Full', n, n, zero, zero, work, n )
289  CALL slacpy( cuplo, n, n, a, lda, work, n )
290 *
291  DO 10 j = 1, n
292  CALL ssyr( cuplo, n, -d( j ), u( 1, j ), 1, work, n )
293  10 CONTINUE
294 *
295  IF( n.GT.1 .AND. kband.EQ.1 ) THEN
296  DO 20 j = 1, n - 1
297  CALL ssyr2( cuplo, n, -e( j ), u( 1, j ), 1, u( 1, j+1 ),
298  $ 1, work, n )
299  20 CONTINUE
300  END IF
301  wnorm = slansy( '1', cuplo, n, work, n, work( n**2+1 ) )
302 *
303  ELSE IF( itype.EQ.2 ) THEN
304 *
305 * ITYPE=2: error = V S V' - A
306 *
307  CALL slaset( 'Full', n, n, zero, zero, work, n )
308 *
309  IF( lower ) THEN
310  work( n**2 ) = d( n )
311  DO 40 j = n - 1, 1, -1
312  IF( kband.EQ.1 ) THEN
313  work( ( n+1 )*( j-1 )+2 ) = ( one-tau( j ) )*e( j )
314  DO 30 jr = j + 2, n
315  work( ( j-1 )*n+jr ) = -tau( j )*e( j )*v( jr, j )
316  30 CONTINUE
317  END IF
318 *
319  vsave = v( j+1, j )
320  v( j+1, j ) = one
321  CALL slarfy( 'L', n-j, v( j+1, j ), 1, tau( j ),
322  $ work( ( n+1 )*j+1 ), n, work( n**2+1 ) )
323  v( j+1, j ) = vsave
324  work( ( n+1 )*( j-1 )+1 ) = d( j )
325  40 CONTINUE
326  ELSE
327  work( 1 ) = d( 1 )
328  DO 60 j = 1, n - 1
329  IF( kband.EQ.1 ) THEN
330  work( ( n+1 )*j ) = ( one-tau( j ) )*e( j )
331  DO 50 jr = 1, j - 1
332  work( j*n+jr ) = -tau( j )*e( j )*v( jr, j+1 )
333  50 CONTINUE
334  END IF
335 *
336  vsave = v( j, j+1 )
337  v( j, j+1 ) = one
338  CALL slarfy( 'U', j, v( 1, j+1 ), 1, tau( j ), work, n,
339  $ work( n**2+1 ) )
340  v( j, j+1 ) = vsave
341  work( ( n+1 )*j+1 ) = d( j+1 )
342  60 CONTINUE
343  END IF
344 *
345  DO 90 jcol = 1, n
346  IF( lower ) THEN
347  DO 70 jrow = jcol, n
348  work( jrow+n*( jcol-1 ) ) = work( jrow+n*( jcol-1 ) )
349  $ - a( jrow, jcol )
350  70 CONTINUE
351  ELSE
352  DO 80 jrow = 1, jcol
353  work( jrow+n*( jcol-1 ) ) = work( jrow+n*( jcol-1 ) )
354  $ - a( jrow, jcol )
355  80 CONTINUE
356  END IF
357  90 CONTINUE
358  wnorm = slansy( '1', cuplo, n, work, n, work( n**2+1 ) )
359 *
360  ELSE IF( itype.EQ.3 ) THEN
361 *
362 * ITYPE=3: error = U V' - I
363 *
364  IF( n.LT.2 )
365  $ RETURN
366  CALL slacpy( ' ', n, n, u, ldu, work, n )
367  IF( lower ) THEN
368  CALL sorm2r( 'R', 'T', n, n-1, n-1, v( 2, 1 ), ldv, tau,
369  $ work( n+1 ), n, work( n**2+1 ), iinfo )
370  ELSE
371  CALL sorm2l( 'R', 'T', n, n-1, n-1, v( 1, 2 ), ldv, tau,
372  $ work, n, work( n**2+1 ), iinfo )
373  END IF
374  IF( iinfo.NE.0 ) THEN
375  result( 1 ) = ten / ulp
376  RETURN
377  END IF
378 *
379  DO 100 j = 1, n
380  work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - one
381  100 CONTINUE
382 *
383  wnorm = slange( '1', n, n, work, n, work( n**2+1 ) )
384  END IF
385 *
386  IF( anorm.GT.wnorm ) THEN
387  result( 1 ) = ( wnorm / anorm ) / ( n*ulp )
388  ELSE
389  IF( anorm.LT.one ) THEN
390  result( 1 ) = ( min( wnorm, n*anorm ) / anorm ) / ( n*ulp )
391  ELSE
392  result( 1 ) = min( wnorm / anorm, REAL( N ) ) / ( n*ulp )
393  END IF
394  END IF
395 *
396 * Do Test 2
397 *
398 * Compute UU' - I
399 *
400  IF( itype.EQ.1 ) THEN
401  CALL sgemm( 'N', 'C', n, n, n, one, u, ldu, u, ldu, zero, work,
402  $ n )
403 *
404  DO 110 j = 1, n
405  work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - one
406  110 CONTINUE
407 *
408  result( 2 ) = min( slange( '1', n, n, work, n,
409  $ work( n**2+1 ) ), REAL( N ) ) / ( n*ulp )
410  END IF
411 *
412  RETURN
413 *
414 * End of SSYT21
415 *
ssyr2
subroutine ssyr2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
SSYR2
Definition: ssyr2.f:149
sgemm
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
slacpy
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
slaset
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
slange
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
sorm2r
subroutine sorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
SORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
Definition: sorm2r.f:161
sorm2l
subroutine sorm2l(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
SORM2L multiplies a general matrix by the orthogonal matrix from a QL factorization determined by sge...
Definition: sorm2l.f:161
slarfy
subroutine slarfy(UPLO, N, V, INCV, TAU, C, LDC, WORK)
SLARFY
Definition: slarfy.f:110
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
ssyr
subroutine ssyr(UPLO, N, ALPHA, X, INCX, A, LDA)
SSYR
Definition: ssyr.f:134
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
slansy
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: