LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ ssyt21()

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**T

 where **T 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**T | / ( |A| n ulp ) and
    RESULT(2) = | I - U U**T | / ( n ulp )

 If ITYPE=2, then:

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

 If ITYPE=3, then:

    RESULT(1) = | I - V U**T | / ( 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)**T 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**T | / ( |A| n ulp ) and
             RESULT(2) = | I - U U**T | / ( n ulp )

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

          3: U expressed both as a dense orthogonal matrix and
             as a product of Housholder transformations:
             RESULT(1) = | I - V U**T | / ( 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)**T 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.

Definition at line 205 of file ssyt21.f.

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