LAPACK 3.12.1
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  and whose scaling constants are in .
!> We shall use the letter  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*
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
ssyr2
subroutine ssyr2(uplo, n, alpha, x, incx, y, incy, a, lda)
SSYR2
Definition ssyr2.f:147
ssyr
subroutine ssyr(uplo, n, alpha, x, incx, a, lda)
SSYR
Definition ssyr.f:132
slacpy
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:101
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slange
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:112
slansy
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:120
slarfy
subroutine slarfy(uplo, n, v, incv, tau, c, ldc, work)
SLARFY
Definition slarfy.f:108
slaset
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:108
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
sorm2l
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:157
sorm2r
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:157
Here is the call graph for this function:
Here is the caller graph for this function: