LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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'

 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, 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)' 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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 207 of file dsyt21.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  DOUBLE PRECISION a( lda, * ), d( * ), e( * ), result( 2 ),
219  $ tau( * ), u( ldu, * ), v( ldv, * ), work( * )
220 * ..
221 *
222 * =====================================================================
223 *
224 * .. Parameters ..
225  DOUBLE PRECISION zero, one, ten
226  parameter ( zero = 0.0d0, one = 1.0d0, ten = 10.0d0 )
227 * ..
228 * .. Local Scalars ..
229  LOGICAL lower
230  CHARACTER cuplo
231  INTEGER iinfo, j, jcol, jr, jrow
232  DOUBLE PRECISION anorm, ulp, unfl, vsave, wnorm
233 * ..
234 * .. External Functions ..
235  LOGICAL lsame
236  DOUBLE PRECISION dlamch, dlange, dlansy
237  EXTERNAL lsame, dlamch, dlange, dlansy
238 * ..
239 * .. External Subroutines ..
240  EXTERNAL dgemm, dlacpy, dlarfy, dlaset, dorm2l, dorm2r,
241  $ dsyr, dsyr2
242 * ..
243 * .. Intrinsic Functions ..
244  INTRINSIC dble, max, min
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 = dlamch( 'Safe minimum' )
263  ulp = dlamch( 'Epsilon' )*dlamch( '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( dlansy( '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 dlaset( 'Full', n, n, zero, zero, work, n )
289  CALL dlacpy( cuplo, n, n, a, lda, work, n )
290 *
291  DO 10 j = 1, n
292  CALL dsyr( 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 dsyr2( cuplo, n, -e( j ), u( 1, j ), 1, u( 1, j+1 ),
298  $ 1, work, n )
299  20 CONTINUE
300  END IF
301  wnorm = dlansy( '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 dlaset( '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 dlarfy( '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 dlarfy( '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 = dlansy( '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 dlacpy( ' ', n, n, u, ldu, work, n )
367  IF( lower ) THEN
368  CALL dorm2r( '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 dorm2l( '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 = dlange( '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, dble( 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 dgemm( '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( dlange( '1', n, n, work, n,
409  $ work( n**2+1 ) ), dble( n ) ) / ( n*ulp )
410  END IF
411 *
412  RETURN
413 *
414 * End of DSYT21
415 *
dlaset
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:112
dlansy
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, or the element of largest absolute value of a real symmetric matrix.
Definition: dlansy.f:124
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
dsyr
subroutine dsyr(UPLO, N, ALPHA, X, INCX, A, LDA)
DSYR
Definition: dsyr.f:134
dlacpy
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
dlange
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:116
dorm2r
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:161
dsyr2
subroutine dsyr2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DSYR2
Definition: dsyr2.f:149
dorm2l
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:161
dlarfy
subroutine dlarfy(UPLO, N, V, INCV, TAU, C, LDC, WORK)
DLARFY
Definition: dlarfy.f:110
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: