LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ssyt21.f
Go to the documentation of this file.
1*> \brief \b SSYT21
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE SSYT21( ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V,
12* LDV, TAU, WORK, RESULT )
13*
14* .. Scalar Arguments ..
15* CHARACTER UPLO
16* INTEGER ITYPE, KBAND, LDA, LDU, LDV, N
17* ..
18* .. Array Arguments ..
19* REAL A( LDA, * ), D( * ), E( * ), RESULT( 2 ),
20* $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
21* ..
22*
23*
24*> \par Purpose:
25* =============
26*>
27*> \verbatim
28*>
29*> SSYT21 generally checks a decomposition of the form
30*>
31*> A = U S U**T
32*>
33*> where **T means transpose, A is symmetric, U is orthogonal, and S is
34*> diagonal (if KBAND=0) or symmetric tridiagonal (if KBAND=1).
35*>
36*> If ITYPE=1, then U is represented as a dense matrix; otherwise U is
37*> expressed as a product of Householder transformations, whose vectors
38*> are stored in the array "V" and whose scaling constants are in "TAU".
39*> We shall use the letter "V" to refer to the product of Householder
40*> transformations (which should be equal to U).
41*>
42*> Specifically, if ITYPE=1, then:
43*>
44*> RESULT(1) = | A - U S U**T | / ( |A| n ulp ) and
45*> RESULT(2) = | I - U U**T | / ( n ulp )
46*>
47*> If ITYPE=2, then:
48*>
49*> RESULT(1) = | A - V S V**T | / ( |A| n ulp )
50*>
51*> If ITYPE=3, then:
52*>
53*> RESULT(1) = | I - V U**T | / ( n ulp )
54*>
55*> For ITYPE > 1, the transformation U is expressed as a product
56*> V = H(1)...H(n-2), where H(j) = I - tau(j) v(j) v(j)**T and each
57*> vector v(j) has its first j elements 0 and the remaining n-j elements
58*> stored in V(j+1:n,j).
59*> \endverbatim
60*
61* Arguments:
62* ==========
63*
64*> \param[in] ITYPE
65*> \verbatim
66*> ITYPE is INTEGER
67*> Specifies the type of tests to be performed.
68*> 1: U expressed as a dense orthogonal matrix:
69*> RESULT(1) = | A - U S U**T | / ( |A| n ulp ) and
70*> RESULT(2) = | I - U U**T | / ( n ulp )
71*>
72*> 2: U expressed as a product V of Housholder transformations:
73*> RESULT(1) = | A - V S V**T | / ( |A| n ulp )
74*>
75*> 3: U expressed both as a dense orthogonal matrix and
76*> as a product of Housholder transformations:
77*> RESULT(1) = | I - V U**T | / ( n ulp )
78*> \endverbatim
79*>
80*> \param[in] UPLO
81*> \verbatim
82*> UPLO is CHARACTER
83*> If UPLO='U', the upper triangle of A and V will be used and
84*> the (strictly) lower triangle will not be referenced.
85*> If UPLO='L', the lower triangle of A and V will be used and
86*> the (strictly) upper triangle will not be referenced.
87*> \endverbatim
88*>
89*> \param[in] N
90*> \verbatim
91*> N is INTEGER
92*> The size of the matrix. If it is zero, SSYT21 does nothing.
93*> It must be at least zero.
94*> \endverbatim
95*>
96*> \param[in] KBAND
97*> \verbatim
98*> KBAND is INTEGER
99*> The bandwidth of the matrix. It may only be zero or one.
100*> If zero, then S is diagonal, and E is not referenced. If
101*> one, then S is symmetric tri-diagonal.
102*> \endverbatim
103*>
104*> \param[in] A
105*> \verbatim
106*> A is REAL array, dimension (LDA, N)
107*> The original (unfactored) matrix. It is assumed to be
108*> symmetric, and only the upper (UPLO='U') or only the lower
109*> (UPLO='L') will be referenced.
110*> \endverbatim
111*>
112*> \param[in] LDA
113*> \verbatim
114*> LDA is INTEGER
115*> The leading dimension of A. It must be at least 1
116*> and at least N.
117*> \endverbatim
118*>
119*> \param[in] D
120*> \verbatim
121*> D is REAL array, dimension (N)
122*> The diagonal of the (symmetric tri-) diagonal matrix.
123*> \endverbatim
124*>
125*> \param[in] E
126*> \verbatim
127*> E is REAL array, dimension (N-1)
128*> The off-diagonal of the (symmetric tri-) diagonal matrix.
129*> E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and
130*> (3,2) element, etc.
131*> Not referenced if KBAND=0.
132*> \endverbatim
133*>
134*> \param[in] U
135*> \verbatim
136*> U is REAL array, dimension (LDU, N)
137*> If ITYPE=1 or 3, this contains the orthogonal matrix in
138*> the decomposition, expressed as a dense matrix. If ITYPE=2,
139*> then it is not referenced.
140*> \endverbatim
141*>
142*> \param[in] LDU
143*> \verbatim
144*> LDU is INTEGER
145*> The leading dimension of U. LDU must be at least N and
146*> at least 1.
147*> \endverbatim
148*>
149*> \param[in] V
150*> \verbatim
151*> V is REAL array, dimension (LDV, N)
152*> If ITYPE=2 or 3, the columns of this array contain the
153*> Householder vectors used to describe the orthogonal matrix
154*> in the decomposition. If UPLO='L', then the vectors are in
155*> the lower triangle, if UPLO='U', then in the upper
156*> triangle.
157*> *NOTE* If ITYPE=2 or 3, V is modified and restored. The
158*> subdiagonal (if UPLO='L') or the superdiagonal (if UPLO='U')
159*> is set to one, and later reset to its original value, during
160*> the course of the calculation.
161*> If ITYPE=1, then it is neither referenced nor modified.
162*> \endverbatim
163*>
164*> \param[in] LDV
165*> \verbatim
166*> LDV is INTEGER
167*> The leading dimension of V. LDV must be at least N and
168*> at least 1.
169*> \endverbatim
170*>
171*> \param[in] TAU
172*> \verbatim
173*> TAU is REAL array, dimension (N)
174*> If ITYPE >= 2, then TAU(j) is the scalar factor of
175*> v(j) v(j)**T in the Householder transformation H(j) of
176*> the product U = H(1)...H(n-2)
177*> If ITYPE < 2, then TAU is not referenced.
178*> \endverbatim
179*>
180*> \param[out] WORK
181*> \verbatim
182*> WORK is REAL array, dimension (2*N**2)
183*> \endverbatim
184*>
185*> \param[out] RESULT
186*> \verbatim
187*> RESULT is REAL array, dimension (2)
188*> The values computed by the two tests described above. The
189*> values are currently limited to 1/ulp, to avoid overflow.
190*> RESULT(1) is always modified. RESULT(2) is modified only
191*> if ITYPE=1.
192*> \endverbatim
193*
194* Authors:
195* ========
196*
197*> \author Univ. of Tennessee
198*> \author Univ. of California Berkeley
199*> \author Univ. of Colorado Denver
200*> \author NAG Ltd.
201*
202*> \ingroup single_eig
203*
204* =====================================================================
205 SUBROUTINE ssyt21( ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V,
206 $ LDV, TAU, WORK, RESULT )
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*
415 END
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
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
subroutine slarfy(UPLO, N, V, INCV, TAU, C, LDC, WORK)
SLARFY
Definition: slarfy.f:108
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
subroutine ssyr(UPLO, N, ALPHA, X, INCX, A, LDA)
SSYR
Definition: ssyr.f:132
subroutine ssyr2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
SSYR2
Definition: ssyr2.f:147
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
subroutine ssyt21(ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V, LDV, TAU, WORK, RESULT)
SSYT21
Definition: ssyt21.f:207