 LAPACK  3.9.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.```
Date
December 2016

Definition at line 209 of file ssyt21.f.

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