LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages

◆ 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*
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:101
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:112
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
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:108
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:157
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: