LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ dsyt21()

 subroutine dsyt21 ( integer ITYPE, character UPLO, integer N, integer KBAND, double precision, dimension( lda, * ) A, integer LDA, double precision, dimension( * ) D, double precision, dimension( * ) E, double precision, dimension( ldu, * ) U, integer LDU, double precision, dimension( ldv, * ) V, integer LDV, double precision, dimension( * ) TAU, double precision, dimension( * ) WORK, double precision, dimension( 2 ) RESULT )

DSYT21

Purpose:
``` DSYT21 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, DSYT21 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N) The diagonal of the (symmetric tri-) diagonal matrix.``` [in] E ``` E is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (2*N**2)` [out] RESULT ``` RESULT is DOUBLE PRECISION 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.```

Definition at line 205 of file dsyt21.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  DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), RESULT( 2 ),
218  \$ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
219 * ..
220 *
221 * =====================================================================
222 *
223 * .. Parameters ..
224  DOUBLE PRECISION ZERO, ONE, TEN
225  parameter( zero = 0.0d0, one = 1.0d0, ten = 10.0d0 )
226 * ..
227 * .. Local Scalars ..
228  LOGICAL LOWER
229  CHARACTER CUPLO
230  INTEGER IINFO, J, JCOL, JR, JROW
231  DOUBLE PRECISION ANORM, ULP, UNFL, VSAVE, WNORM
232 * ..
233 * .. External Functions ..
234  LOGICAL LSAME
235  DOUBLE PRECISION DLAMCH, DLANGE, DLANSY
236  EXTERNAL lsame, dlamch, dlange, dlansy
237 * ..
238 * .. External Subroutines ..
239  EXTERNAL dgemm, dlacpy, dlarfy, dlaset, dorm2l, dorm2r,
240  \$ dsyr, dsyr2
241 * ..
242 * .. Intrinsic Functions ..
243  INTRINSIC dble, max, min
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 = dlamch( 'Safe minimum' )
262  ulp = dlamch( 'Epsilon' )*dlamch( '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( dlansy( '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 dlaset( 'Full', n, n, zero, zero, work, n )
288  CALL dlacpy( cuplo, n, n, a, lda, work, n )
289 *
290  DO 10 j = 1, n
291  CALL dsyr( 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 dsyr2( cuplo, n, -e( j ), u( 1, j ), 1, u( 1, j+1 ),
297  \$ 1, work, n )
298  20 CONTINUE
299  END IF
300  wnorm = dlansy( '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 dlaset( '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 dlarfy( '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 dlarfy( '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 = dlansy( '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 dlacpy( ' ', n, n, u, ldu, work, n )
366  IF( lower ) THEN
367  CALL dorm2r( '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 dorm2l( '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 = dlange( '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, dble( 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 dgemm( '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( dlange( '1', n, n, work, n,
408  \$ work( n**2+1 ) ), dble( n ) ) / ( n*ulp )
409  END IF
410 *
411  RETURN
412 *
413 * End of DSYT21
414 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dsyr(UPLO, N, ALPHA, X, INCX, A, LDA)
DSYR
Definition: dsyr.f:132
subroutine dsyr2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DSYR2
Definition: dsyr2.f:147
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: dlange.f:114
subroutine dlarfy(UPLO, N, V, INCV, TAU, C, LDC, WORK)
DLARFY
Definition: dlarfy.f:108
subroutine dorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
DORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
Definition: dorm2r.f:159
subroutine dorm2l(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
DORM2L multiplies a general matrix by the orthogonal matrix from a QL factorization determined by sge...
Definition: dorm2l.f:159
double precision function dlansy(NORM, UPLO, N, A, LDA, WORK)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: dlansy.f:122
Here is the call graph for this function:
Here is the caller graph for this function: