LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dsyt21.f
Go to the documentation of this file.
1 *> \brief \b DSYT21
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 DSYT21( 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 * DOUBLE PRECISION 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 *> DSYT21 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, DSYT21 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
122 *> The diagonal of the (symmetric tri-) diagonal matrix.
123 *> \endverbatim
124 *>
125 *> \param[in] E
126 *> \verbatim
127 *> E is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (2*N**2)
183 *> \endverbatim
184 *>
185 *> \param[out] RESULT
186 *> \verbatim
187 *> RESULT is DOUBLE PRECISION 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 double_eig
203 *
204 * =====================================================================
205  SUBROUTINE dsyt21( 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  DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), RESULT( 2 ),
218  $ tau( * ), u( ldu, * ), v( ldv, * ), work( * )
219 * ..
220 *
221 * =====================================================================
222 *
223 * .. Parameters ..
224  DOUBLE PRECISION ZERO, ONE, TEN
225  parameter( zero = 0.0d0, one = 1.0d0, ten = 10.0d0 )
226 * ..
227 * .. Local Scalars ..
228  LOGICAL LOWER
229  CHARACTER CUPLO
230  INTEGER IINFO, J, JCOL, JR, JROW
231  DOUBLE PRECISION ANORM, ULP, UNFL, VSAVE, WNORM
232 * ..
233 * .. External Functions ..
234  LOGICAL LSAME
235  DOUBLE PRECISION DLAMCH, DLANGE, DLANSY
236  EXTERNAL lsame, dlamch, dlange, dlansy
237 * ..
238 * .. External Subroutines ..
239  EXTERNAL dgemm, dlacpy, dlarfy, dlaset, dorm2l, dorm2r,
240  $ dsyr, dsyr2
241 * ..
242 * .. Intrinsic Functions ..
243  INTRINSIC dble, max, min
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 = dlamch( 'Safe minimum' )
262  ulp = dlamch( 'Epsilon' )*dlamch( '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( dlansy( '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 dlaset( 'Full', n, n, zero, zero, work, n )
288  CALL dlacpy( cuplo, n, n, a, lda, work, n )
289 *
290  DO 10 j = 1, n
291  CALL dsyr( 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 dsyr2( cuplo, n, -e( j ), u( 1, j ), 1, u( 1, j+1 ),
297  $ 1, work, n )
298  20 CONTINUE
299  END IF
300  wnorm = dlansy( '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 dlaset( '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 dlarfy( '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 dlarfy( '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 = dlansy( '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 dlacpy( ' ', n, n, u, ldu, work, n )
366  IF( lower ) THEN
367  CALL dorm2r( '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 dorm2l( '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 = dlange( '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, dble( 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 dgemm( '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( dlange( '1', n, n, work, n,
408  $ work( n**2+1 ) ), dble( n ) ) / ( n*ulp )
409  END IF
410 *
411  RETURN
412 *
413 * End of DSYT21
414 *
415  END
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine dsyr(UPLO, N, ALPHA, X, INCX, A, LDA)
DSYR
Definition: dsyr.f:132
subroutine dsyr2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DSYR2
Definition: dsyr2.f:147
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
subroutine dsyt21(ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V, LDV, TAU, WORK, RESULT)
DSYT21
Definition: dsyt21.f:207
subroutine dlarfy(UPLO, N, V, INCV, TAU, C, LDC, WORK)
DLARFY
Definition: dlarfy.f:108
subroutine dorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
DORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
Definition: dorm2r.f:159
subroutine dorm2l(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
DORM2L multiplies a general matrix by the orthogonal matrix from a QL factorization determined by sge...
Definition: dorm2l.f:159