LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dstt22.f
Go to the documentation of this file.
1*> \brief \b DSTT22
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 DSTT22( N, M, KBAND, AD, AE, SD, SE, U, LDU, WORK,
12* LDWORK, RESULT )
13*
14* .. Scalar Arguments ..
15* INTEGER KBAND, LDU, LDWORK, M, N
16* ..
17* .. Array Arguments ..
18* DOUBLE PRECISION AD( * ), AE( * ), RESULT( 2 ), SD( * ),
19* $ SE( * ), U( LDU, * ), WORK( LDWORK, * )
20* ..
21*
22*
23*> \par Purpose:
24* =============
25*>
26*> \verbatim
27*>
28*> DSTT22 checks a set of M eigenvalues and eigenvectors,
29*>
30*> A U = U S
31*>
32*> where A is symmetric tridiagonal, the columns of U are orthogonal,
33*> and S is diagonal (if KBAND=0) or symmetric tridiagonal (if KBAND=1).
34*> Two tests are performed:
35*>
36*> RESULT(1) = | U' A U - S | / ( |A| m ulp )
37*>
38*> RESULT(2) = | I - U'U | / ( m ulp )
39*> \endverbatim
40*
41* Arguments:
42* ==========
43*
44*> \param[in] N
45*> \verbatim
46*> N is INTEGER
47*> The size of the matrix. If it is zero, DSTT22 does nothing.
48*> It must be at least zero.
49*> \endverbatim
50*>
51*> \param[in] M
52*> \verbatim
53*> M is INTEGER
54*> The number of eigenpairs to check. If it is zero, DSTT22
55*> does nothing. It must be at least zero.
56*> \endverbatim
57*>
58*> \param[in] KBAND
59*> \verbatim
60*> KBAND is INTEGER
61*> The bandwidth of the matrix S. It may only be zero or one.
62*> If zero, then S is diagonal, and SE is not referenced. If
63*> one, then S is symmetric tri-diagonal.
64*> \endverbatim
65*>
66*> \param[in] AD
67*> \verbatim
68*> AD is DOUBLE PRECISION array, dimension (N)
69*> The diagonal of the original (unfactored) matrix A. A is
70*> assumed to be symmetric tridiagonal.
71*> \endverbatim
72*>
73*> \param[in] AE
74*> \verbatim
75*> AE is DOUBLE PRECISION array, dimension (N)
76*> The off-diagonal of the original (unfactored) matrix A. A
77*> is assumed to be symmetric tridiagonal. AE(1) is ignored,
78*> AE(2) is the (1,2) and (2,1) element, etc.
79*> \endverbatim
80*>
81*> \param[in] SD
82*> \verbatim
83*> SD is DOUBLE PRECISION array, dimension (N)
84*> The diagonal of the (symmetric tri-) diagonal matrix S.
85*> \endverbatim
86*>
87*> \param[in] SE
88*> \verbatim
89*> SE is DOUBLE PRECISION array, dimension (N)
90*> The off-diagonal of the (symmetric tri-) diagonal matrix S.
91*> Not referenced if KBSND=0. If KBAND=1, then AE(1) is
92*> ignored, SE(2) is the (1,2) and (2,1) element, etc.
93*> \endverbatim
94*>
95*> \param[in] U
96*> \verbatim
97*> U is DOUBLE PRECISION array, dimension (LDU, N)
98*> The orthogonal matrix in the decomposition.
99*> \endverbatim
100*>
101*> \param[in] LDU
102*> \verbatim
103*> LDU is INTEGER
104*> The leading dimension of U. LDU must be at least N.
105*> \endverbatim
106*>
107*> \param[out] WORK
108*> \verbatim
109*> WORK is DOUBLE PRECISION array, dimension (LDWORK, M+1)
110*> \endverbatim
111*>
112*> \param[in] LDWORK
113*> \verbatim
114*> LDWORK is INTEGER
115*> The leading dimension of WORK. LDWORK must be at least
116*> max(1,M).
117*> \endverbatim
118*>
119*> \param[out] RESULT
120*> \verbatim
121*> RESULT is DOUBLE PRECISION array, dimension (2)
122*> The values computed by the two tests described above. The
123*> values are currently limited to 1/ulp, to avoid overflow.
124*> \endverbatim
125*
126* Authors:
127* ========
128*
129*> \author Univ. of Tennessee
130*> \author Univ. of California Berkeley
131*> \author Univ. of Colorado Denver
132*> \author NAG Ltd.
133*
134*> \ingroup double_eig
135*
136* =====================================================================
137 SUBROUTINE dstt22( N, M, KBAND, AD, AE, SD, SE, U, LDU, WORK,
138 $ LDWORK, RESULT )
139*
140* -- LAPACK test routine --
141* -- LAPACK is a software package provided by Univ. of Tennessee, --
142* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
143*
144* .. Scalar Arguments ..
145 INTEGER KBAND, LDU, LDWORK, M, N
146* ..
147* .. Array Arguments ..
148 DOUBLE PRECISION AD( * ), AE( * ), RESULT( 2 ), SD( * ),
149 $ se( * ), u( ldu, * ), work( ldwork, * )
150* ..
151*
152* =====================================================================
153*
154* .. Parameters ..
155 DOUBLE PRECISION ZERO, ONE
156 parameter( zero = 0.0d0, one = 1.0d0 )
157* ..
158* .. Local Scalars ..
159 INTEGER I, J, K
160 DOUBLE PRECISION ANORM, AUKJ, ULP, UNFL, WNORM
161* ..
162* .. External Functions ..
163 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY
164 EXTERNAL dlamch, dlange, dlansy
165* ..
166* .. External Subroutines ..
167 EXTERNAL dgemm
168* ..
169* .. Intrinsic Functions ..
170 INTRINSIC abs, dble, max, min
171* ..
172* .. Executable Statements ..
173*
174 result( 1 ) = zero
175 result( 2 ) = zero
176 IF( n.LE.0 .OR. m.LE.0 )
177 $ RETURN
178*
179 unfl = dlamch( 'Safe minimum' )
180 ulp = dlamch( 'Epsilon' )
181*
182* Do Test 1
183*
184* Compute the 1-norm of A.
185*
186 IF( n.GT.1 ) THEN
187 anorm = abs( ad( 1 ) ) + abs( ae( 1 ) )
188 DO 10 j = 2, n - 1
189 anorm = max( anorm, abs( ad( j ) )+abs( ae( j ) )+
190 $ abs( ae( j-1 ) ) )
191 10 CONTINUE
192 anorm = max( anorm, abs( ad( n ) )+abs( ae( n-1 ) ) )
193 ELSE
194 anorm = abs( ad( 1 ) )
195 END IF
196 anorm = max( anorm, unfl )
197*
198* Norm of U'AU - S
199*
200 DO 40 i = 1, m
201 DO 30 j = 1, m
202 work( i, j ) = zero
203 DO 20 k = 1, n
204 aukj = ad( k )*u( k, j )
205 IF( k.NE.n )
206 $ aukj = aukj + ae( k )*u( k+1, j )
207 IF( k.NE.1 )
208 $ aukj = aukj + ae( k-1 )*u( k-1, j )
209 work( i, j ) = work( i, j ) + u( k, i )*aukj
210 20 CONTINUE
211 30 CONTINUE
212 work( i, i ) = work( i, i ) - sd( i )
213 IF( kband.EQ.1 ) THEN
214 IF( i.NE.1 )
215 $ work( i, i-1 ) = work( i, i-1 ) - se( i-1 )
216 IF( i.NE.n )
217 $ work( i, i+1 ) = work( i, i+1 ) - se( i )
218 END IF
219 40 CONTINUE
220*
221 wnorm = dlansy( '1', 'L', m, work, m, work( 1, m+1 ) )
222*
223 IF( anorm.GT.wnorm ) THEN
224 result( 1 ) = ( wnorm / anorm ) / ( m*ulp )
225 ELSE
226 IF( anorm.LT.one ) THEN
227 result( 1 ) = ( min( wnorm, m*anorm ) / anorm ) / ( m*ulp )
228 ELSE
229 result( 1 ) = min( wnorm / anorm, dble( m ) ) / ( m*ulp )
230 END IF
231 END IF
232*
233* Do Test 2
234*
235* Compute U'U - I
236*
237 CALL dgemm( 'T', 'N', m, m, n, one, u, ldu, u, ldu, zero, work,
238 $ m )
239*
240 DO 50 j = 1, m
241 work( j, j ) = work( j, j ) - one
242 50 CONTINUE
243*
244 result( 2 ) = min( dble( m ), dlange( '1', m, m, work, m, work( 1,
245 $ m+1 ) ) ) / ( m*ulp )
246*
247 RETURN
248*
249* End of DSTT22
250*
251 END
subroutine dstt22(n, m, kband, ad, ae, sd, se, u, ldu, work, ldwork, result)
DSTT22
Definition dstt22.f:139
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188