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

◆ 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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 205 of file ssyt21.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 REAL A( LDA, * ), D( * ), E( * ), RESULT( 2 ),
218 $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
219* ..
220*
221* =====================================================================
222*
223* .. Parameters ..
224 REAL ZERO, ONE, TEN
225 parameter( zero = 0.0e0, one = 1.0e0, ten = 10.0e0 )
226* ..
227* .. Local Scalars ..
228 LOGICAL LOWER
229 CHARACTER CUPLO
230 INTEGER IINFO, J, JCOL, JR, JROW
231 REAL ANORM, ULP, UNFL, VSAVE, WNORM
232* ..
233* .. External Functions ..
234 LOGICAL LSAME
235 REAL SLAMCH, SLANGE, SLANSY
236 EXTERNAL lsame, slamch, slange, slansy
237* ..
238* .. External Subroutines ..
239 EXTERNAL sgemm, slacpy, slarfy, slaset, sorm2l, sorm2r,
240 $ ssyr, ssyr2
241* ..
242* .. Intrinsic Functions ..
243 INTRINSIC max, min, real
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 = slamch( 'Safe minimum' )
262 ulp = slamch( 'Epsilon' )*slamch( '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( slansy( '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 slaset( 'Full', n, n, zero, zero, work, n )
288 CALL slacpy( cuplo, n, n, a, lda, work, n )
289*
290 DO 10 j = 1, n
291 CALL ssyr( 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 ssyr2( cuplo, n, -e( j ), u( 1, j ), 1, u( 1, j+1 ),
297 $ 1, work, n )
298 20 CONTINUE
299 END IF
300 wnorm = slansy( '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 slaset( '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 slarfy( '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 slarfy( '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 = slansy( '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 slacpy( ' ', n, n, u, ldu, work, n )
366 IF( lower ) THEN
367 CALL sorm2r( '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 sorm2l( '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 = slange( '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, real( 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 sgemm( '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( slange( '1', n, n, work, n,
408 $ work( n**2+1 ) ), real( n ) ) / ( n*ulp )
409 END IF
410*
411 RETURN
412*
413* End of SSYT21
414*
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine ssyr2(uplo, n, alpha, x, incx, y, incy, a, lda)
SSYR2
Definition ssyr2.f:147
subroutine ssyr(uplo, n, alpha, x, incx, a, lda)
SSYR
Definition ssyr.f:132
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition slange.f:114
real function slansy(norm, uplo, n, a, lda, work)
SLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition slansy.f:122
subroutine slarfy(uplo, n, v, incv, tau, c, ldc, work)
SLARFY
Definition slarfy.f:108
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:110
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
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:159
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:159
Here is the call graph for this function:
Here is the caller graph for this function: