LAPACK  3.9.0
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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 209 of file dsyt21.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  DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), RESULT( 2 ),
221  $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
222 * ..
223 *
224 * =====================================================================
225 *
226 * .. Parameters ..
227  DOUBLE PRECISION ZERO, ONE, TEN
228  parameter( zero = 0.0d0, one = 1.0d0, ten = 10.0d0 )
229 * ..
230 * .. Local Scalars ..
231  LOGICAL LOWER
232  CHARACTER CUPLO
233  INTEGER IINFO, J, JCOL, JR, JROW
234  DOUBLE PRECISION ANORM, ULP, UNFL, VSAVE, WNORM
235 * ..
236 * .. External Functions ..
237  LOGICAL LSAME
238  DOUBLE PRECISION DLAMCH, DLANGE, DLANSY
239  EXTERNAL lsame, dlamch, dlange, dlansy
240 * ..
241 * .. External Subroutines ..
242  EXTERNAL dgemm, dlacpy, dlarfy, dlaset, dorm2l, dorm2r,
243  $ dsyr, dsyr2
244 * ..
245 * .. Intrinsic Functions ..
246  INTRINSIC dble, max, min
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 = dlamch( 'Safe minimum' )
265  ulp = dlamch( 'Epsilon' )*dlamch( '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( dlansy( '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 dlaset( 'Full', n, n, zero, zero, work, n )
291  CALL dlacpy( cuplo, n, n, a, lda, work, n )
292 *
293  DO 10 j = 1, n
294  CALL dsyr( 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 dsyr2( cuplo, n, -e( j ), u( 1, j ), 1, u( 1, j+1 ),
300  $ 1, work, n )
301  20 CONTINUE
302  END IF
303  wnorm = dlansy( '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 dlaset( '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 dlarfy( '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 dlarfy( '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 = dlansy( '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 dlacpy( ' ', n, n, u, ldu, work, n )
369  IF( lower ) THEN
370  CALL dorm2r( '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 dorm2l( '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 = dlange( '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, dble( 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 dgemm( '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( dlange( '1', n, n, work, n,
411  $ work( n**2+1 ) ), dble( n ) ) / ( n*ulp )
412  END IF
413 *
414  RETURN
415 *
416 * End of DSYT21
417 *
Here is the call graph for this function:
Here is the caller graph for this function:
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
dsyr
subroutine dsyr(UPLO, N, ALPHA, X, INCX, A, LDA)
DSYR
Definition: dsyr.f:134
dsyr2
subroutine dsyr2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DSYR2
Definition: dsyr2.f:149
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
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
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
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