LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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.

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*
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
subroutine dsyr2(uplo, n, alpha, x, incx, y, incy, a, lda)
DSYR2
Definition dsyr2.f:147
subroutine dsyr(uplo, n, alpha, x, incx, a, lda)
DSYR
Definition dsyr.f:132
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
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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
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
subroutine dlarfy(uplo, n, v, incv, tau, c, ldc, work)
DLARFY
Definition dlarfy.f:108
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:48
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
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
Here is the call graph for this function:
Here is the caller graph for this function: