LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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'
32 *>
33 *> where ' 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' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp )
45 *>
46 *> If ITYPE=2, then:
47 *>
48 *> RESULT(1) = | A - V S V' | / ( |A| n ulp )
49 *>
50 *> If ITYPE=3, then:
51 *>
52 *> RESULT(1) = | I - VU' | / ( n ulp )
53 *>
54 *> For ITYPE > 1, the transformation U is expressed as a product
55 *> V = H(1)...H(n-2), where H(j) = I - tau(j) v(j) v(j)' and each
56 *> vector v(j) has its first j elements 0 and the remaining n-j elements
57 *> stored in V(j+1:n,j).
58 *> \endverbatim
59 *
60 * Arguments:
61 * ==========
62 *
63 *> \param[in] ITYPE
64 *> \verbatim
65 *> ITYPE is INTEGER
66 *> Specifies the type of tests to be performed.
67 *> 1: U expressed as a dense orthogonal matrix:
68 *> RESULT(1) = | A - U S U' | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU' | / ( n ulp )
69 *>
70 *> 2: U expressed as a product V of Housholder transformations:
71 *> RESULT(1) = | A - V S V' | / ( |A| n ulp )
72 *>
73 *> 3: U expressed both as a dense orthogonal matrix and
74 *> as a product of Housholder transformations:
75 *> RESULT(1) = | I - VU' | / ( n ulp )
76 *> \endverbatim
77 *>
78 *> \param[in] UPLO
79 *> \verbatim
80 *> UPLO is CHARACTER
81 *> If UPLO='U', the upper triangle of A and V will be used and
82 *> the (strictly) lower triangle will not be referenced.
83 *> If UPLO='L', the lower triangle of A and V will be used and
84 *> the (strictly) upper triangle will not be referenced.
85 *> \endverbatim
86 *>
87 *> \param[in] N
88 *> \verbatim
89 *> N is INTEGER
90 *> The size of the matrix. If it is zero, SSYT21 does nothing.
91 *> It must be at least zero.
92 *> \endverbatim
93 *>
94 *> \param[in] KBAND
95 *> \verbatim
96 *> KBAND is INTEGER
97 *> The bandwidth of the matrix. It may only be zero or one.
98 *> If zero, then S is diagonal, and E is not referenced. If
99 *> one, then S is symmetric tri-diagonal.
100 *> \endverbatim
101 *>
102 *> \param[in] A
103 *> \verbatim
104 *> A is REAL array, dimension (LDA, N)
105 *> The original (unfactored) matrix. It is assumed to be
106 *> symmetric, and only the upper (UPLO='U') or only the lower
107 *> (UPLO='L') will be referenced.
108 *> \endverbatim
109 *>
110 *> \param[in] LDA
111 *> \verbatim
112 *> LDA is INTEGER
113 *> The leading dimension of A. It must be at least 1
114 *> and at least N.
115 *> \endverbatim
116 *>
117 *> \param[in] D
118 *> \verbatim
119 *> D is REAL array, dimension (N)
120 *> The diagonal of the (symmetric tri-) diagonal matrix.
121 *> \endverbatim
122 *>
123 *> \param[in] E
124 *> \verbatim
125 *> E is REAL array, dimension (N-1)
126 *> The off-diagonal of the (symmetric tri-) diagonal matrix.
127 *> E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and
128 *> (3,2) element, etc.
129 *> Not referenced if KBAND=0.
130 *> \endverbatim
131 *>
132 *> \param[in] U
133 *> \verbatim
134 *> U is REAL array, dimension (LDU, N)
135 *> If ITYPE=1 or 3, this contains the orthogonal matrix in
136 *> the decomposition, expressed as a dense matrix. If ITYPE=2,
137 *> then it is not referenced.
138 *> \endverbatim
139 *>
140 *> \param[in] LDU
141 *> \verbatim
142 *> LDU is INTEGER
143 *> The leading dimension of U. LDU must be at least N and
144 *> at least 1.
145 *> \endverbatim
146 *>
147 *> \param[in] V
148 *> \verbatim
149 *> V is REAL array, dimension (LDV, N)
150 *> If ITYPE=2 or 3, the columns of this array contain the
151 *> Householder vectors used to describe the orthogonal matrix
152 *> in the decomposition. If UPLO='L', then the vectors are in
153 *> the lower triangle, if UPLO='U', then in the upper
154 *> triangle.
155 *> *NOTE* If ITYPE=2 or 3, V is modified and restored. The
156 *> subdiagonal (if UPLO='L') or the superdiagonal (if UPLO='U')
157 *> is set to one, and later reset to its original value, during
158 *> the course of the calculation.
159 *> If ITYPE=1, then it is neither referenced nor modified.
160 *> \endverbatim
161 *>
162 *> \param[in] LDV
163 *> \verbatim
164 *> LDV is INTEGER
165 *> The leading dimension of V. LDV must be at least N and
166 *> at least 1.
167 *> \endverbatim
168 *>
169 *> \param[in] TAU
170 *> \verbatim
171 *> TAU is REAL array, dimension (N)
172 *> If ITYPE >= 2, then TAU(j) is the scalar factor of
173 *> v(j) v(j)' in the Householder transformation H(j) of
174 *> the product U = H(1)...H(n-2)
175 *> If ITYPE < 2, then TAU is not referenced.
176 *> \endverbatim
177 *>
178 *> \param[out] WORK
179 *> \verbatim
180 *> WORK is REAL array, dimension (2*N**2)
181 *> \endverbatim
182 *>
183 *> \param[out] RESULT
184 *> \verbatim
185 *> RESULT is REAL array, dimension (2)
186 *> The values computed by the two tests described above. The
187 *> values are currently limited to 1/ulp, to avoid overflow.
188 *> RESULT(1) is always modified. RESULT(2) is modified only
189 *> if ITYPE=1.
190 *> \endverbatim
191 *
192 * Authors:
193 * ========
194 *
195 *> \author Univ. of Tennessee
196 *> \author Univ. of California Berkeley
197 *> \author Univ. of Colorado Denver
198 *> \author NAG Ltd.
199 *
200 *> \date November 2011
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 (version 3.4.0) --
209 * -- LAPACK is a software package provided by Univ. of Tennessee, --
210 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
211 * November 2011
212 *
213 * .. Scalar Arguments ..
214  CHARACTER uplo
215  INTEGER itype, kband, lda, ldu, ldv, n
216 * ..
217 * .. Array Arguments ..
218  REAL a( lda, * ), d( * ), e( * ), result( 2 ),
219  $ tau( * ), u( ldu, * ), v( ldv, * ), work( * )
220 * ..
221 *
222 * =====================================================================
223 *
224 * .. Parameters ..
225  REAL zero, one, ten
226  parameter( zero = 0.0e0, one = 1.0e0, ten = 10.0e0 )
227 * ..
228 * .. Local Scalars ..
229  LOGICAL lower
230  CHARACTER cuplo
231  INTEGER iinfo, j, jcol, jr, jrow
232  REAL anorm, ulp, unfl, vsave, wnorm
233 * ..
234 * .. External Functions ..
235  LOGICAL lsame
236  REAL slamch, slange, slansy
237  EXTERNAL lsame, slamch, slange, slansy
238 * ..
239 * .. External Subroutines ..
240  EXTERNAL sgemm, slacpy, slarfy, slaset, sorm2l, sorm2r,
241  $ ssyr, ssyr2
242 * ..
243 * .. Intrinsic Functions ..
244  INTRINSIC max, min, real
245 * ..
246 * .. Executable Statements ..
247 *
248  result( 1 ) = zero
249  IF( itype.EQ.1 )
250  $ result( 2 ) = zero
251  IF( n.LE.0 )
252  $ return
253 *
254  IF( lsame( uplo, 'U' ) ) THEN
255  lower = .false.
256  cuplo = 'U'
257  ELSE
258  lower = .true.
259  cuplo = 'L'
260  END IF
261 *
262  unfl = slamch( 'Safe minimum' )
263  ulp = slamch( 'Epsilon' )*slamch( 'Base' )
264 *
265 * Some Error Checks
266 *
267  IF( itype.LT.1 .OR. itype.GT.3 ) THEN
268  result( 1 ) = ten / ulp
269  return
270  END IF
271 *
272 * Do Test 1
273 *
274 * Norm of A:
275 *
276  IF( itype.EQ.3 ) THEN
277  anorm = one
278  ELSE
279  anorm = max( slansy( '1', cuplo, n, a, lda, work ), unfl )
280  END IF
281 *
282 * Compute error matrix:
283 *
284  IF( itype.EQ.1 ) THEN
285 *
286 * ITYPE=1: error = A - U S U'
287 *
288  CALL slaset( 'Full', n, n, zero, zero, work, n )
289  CALL slacpy( cuplo, n, n, a, lda, work, n )
290 *
291  DO 10 j = 1, n
292  CALL ssyr( cuplo, n, -d( j ), u( 1, j ), 1, work, n )
293  10 continue
294 *
295  IF( n.GT.1 .AND. kband.EQ.1 ) THEN
296  DO 20 j = 1, n - 1
297  CALL ssyr2( cuplo, n, -e( j ), u( 1, j ), 1, u( 1, j+1 ),
298  $ 1, work, n )
299  20 continue
300  END IF
301  wnorm = slansy( '1', cuplo, n, work, n, work( n**2+1 ) )
302 *
303  ELSE IF( itype.EQ.2 ) THEN
304 *
305 * ITYPE=2: error = V S V' - A
306 *
307  CALL slaset( 'Full', n, n, zero, zero, work, n )
308 *
309  IF( lower ) THEN
310  work( n**2 ) = d( n )
311  DO 40 j = n - 1, 1, -1
312  IF( kband.EQ.1 ) THEN
313  work( ( n+1 )*( j-1 )+2 ) = ( one-tau( j ) )*e( j )
314  DO 30 jr = j + 2, n
315  work( ( j-1 )*n+jr ) = -tau( j )*e( j )*v( jr, j )
316  30 continue
317  END IF
318 *
319  vsave = v( j+1, j )
320  v( j+1, j ) = one
321  CALL slarfy( 'L', n-j, v( j+1, j ), 1, tau( j ),
322  $ work( ( n+1 )*j+1 ), n, work( n**2+1 ) )
323  v( j+1, j ) = vsave
324  work( ( n+1 )*( j-1 )+1 ) = d( j )
325  40 continue
326  ELSE
327  work( 1 ) = d( 1 )
328  DO 60 j = 1, n - 1
329  IF( kband.EQ.1 ) THEN
330  work( ( n+1 )*j ) = ( one-tau( j ) )*e( j )
331  DO 50 jr = 1, j - 1
332  work( j*n+jr ) = -tau( j )*e( j )*v( jr, j+1 )
333  50 continue
334  END IF
335 *
336  vsave = v( j, j+1 )
337  v( j, j+1 ) = one
338  CALL slarfy( 'U', j, v( 1, j+1 ), 1, tau( j ), work, n,
339  $ work( n**2+1 ) )
340  v( j, j+1 ) = vsave
341  work( ( n+1 )*j+1 ) = d( j+1 )
342  60 continue
343  END IF
344 *
345  DO 90 jcol = 1, n
346  IF( lower ) THEN
347  DO 70 jrow = jcol, n
348  work( jrow+n*( jcol-1 ) ) = work( jrow+n*( jcol-1 ) )
349  $ - a( jrow, jcol )
350  70 continue
351  ELSE
352  DO 80 jrow = 1, jcol
353  work( jrow+n*( jcol-1 ) ) = work( jrow+n*( jcol-1 ) )
354  $ - a( jrow, jcol )
355  80 continue
356  END IF
357  90 continue
358  wnorm = slansy( '1', cuplo, n, work, n, work( n**2+1 ) )
359 *
360  ELSE IF( itype.EQ.3 ) THEN
361 *
362 * ITYPE=3: error = U V' - I
363 *
364  IF( n.LT.2 )
365  $ return
366  CALL slacpy( ' ', n, n, u, ldu, work, n )
367  IF( lower ) THEN
368  CALL sorm2r( 'R', 'T', n, n-1, n-1, v( 2, 1 ), ldv, tau,
369  $ work( n+1 ), n, work( n**2+1 ), iinfo )
370  ELSE
371  CALL sorm2l( 'R', 'T', n, n-1, n-1, v( 1, 2 ), ldv, tau,
372  $ work, n, work( n**2+1 ), iinfo )
373  END IF
374  IF( iinfo.NE.0 ) THEN
375  result( 1 ) = ten / ulp
376  return
377  END IF
378 *
379  DO 100 j = 1, n
380  work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - one
381  100 continue
382 *
383  wnorm = slange( '1', n, n, work, n, work( n**2+1 ) )
384  END IF
385 *
386  IF( anorm.GT.wnorm ) THEN
387  result( 1 ) = ( wnorm / anorm ) / ( n*ulp )
388  ELSE
389  IF( anorm.LT.one ) THEN
390  result( 1 ) = ( min( wnorm, n*anorm ) / anorm ) / ( n*ulp )
391  ELSE
392  result( 1 ) = min( wnorm / anorm, REAL( N ) ) / ( n*ulp )
393  END IF
394  END IF
395 *
396 * Do Test 2
397 *
398 * Compute UU' - I
399 *
400  IF( itype.EQ.1 ) THEN
401  CALL sgemm( 'N', 'C', n, n, n, one, u, ldu, u, ldu, zero, work,
402  $ n )
403 *
404  DO 110 j = 1, n
405  work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - one
406  110 continue
407 *
408  result( 2 ) = min( slange( '1', n, n, work, n,
409  $ work( n**2+1 ) ), REAL( N ) ) / ( n*ulp )
410  END IF
411 *
412  return
413 *
414 * End of SSYT21
415 *
416  END