LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
dlatm6.f
Go to the documentation of this file.
1*> \brief \b DLATM6
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 DLATM6( TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA,
12* BETA, WX, WY, S, DIF )
13*
14* .. Scalar Arguments ..
15* INTEGER LDA, LDX, LDY, N, TYPE
16* DOUBLE PRECISION ALPHA, BETA, WX, WY
17* ..
18* .. Array Arguments ..
19* DOUBLE PRECISION A( LDA, * ), B( LDA, * ), DIF( * ), S( * ),
20* \$ X( LDX, * ), Y( LDY, * )
21* ..
22*
23*
24*> \par Purpose:
25* =============
26*>
27*> \verbatim
28*>
29*> DLATM6 generates test matrices for the generalized eigenvalue
30*> problem, their corresponding right and left eigenvector matrices,
31*> and also reciprocal condition numbers for all eigenvalues and
32*> the reciprocal condition numbers of eigenvectors corresponding to
33*> the 1th and 5th eigenvalues.
34*>
35*> Test Matrices
36*> =============
37*>
38*> Two kinds of test matrix pairs
39*>
40*> (A, B) = inverse(YH) * (Da, Db) * inverse(X)
41*>
42*> are used in the tests:
43*>
44*> Type 1:
45*> Da = 1+a 0 0 0 0 Db = 1 0 0 0 0
46*> 0 2+a 0 0 0 0 1 0 0 0
47*> 0 0 3+a 0 0 0 0 1 0 0
48*> 0 0 0 4+a 0 0 0 0 1 0
49*> 0 0 0 0 5+a , 0 0 0 0 1 , and
50*>
51*> Type 2:
52*> Da = 1 -1 0 0 0 Db = 1 0 0 0 0
53*> 1 1 0 0 0 0 1 0 0 0
54*> 0 0 1 0 0 0 0 1 0 0
55*> 0 0 0 1+a 1+b 0 0 0 1 0
56*> 0 0 0 -1-b 1+a , 0 0 0 0 1 .
57*>
58*> In both cases the same inverse(YH) and inverse(X) are used to compute
59*> (A, B), giving the exact eigenvectors to (A,B) as (YH, X):
60*>
61*> YH: = 1 0 -y y -y X = 1 0 -x -x x
62*> 0 1 -y y -y 0 1 x -x -x
63*> 0 0 1 0 0 0 0 1 0 0
64*> 0 0 0 1 0 0 0 0 1 0
65*> 0 0 0 0 1, 0 0 0 0 1 ,
66*>
67*> where a, b, x and y will have all values independently of each other.
68*> \endverbatim
69*
70* Arguments:
71* ==========
72*
73*> \param[in] TYPE
74*> \verbatim
75*> TYPE is INTEGER
76*> Specifies the problem type (see further details).
77*> \endverbatim
78*>
79*> \param[in] N
80*> \verbatim
81*> N is INTEGER
82*> Size of the matrices A and B.
83*> \endverbatim
84*>
85*> \param[out] A
86*> \verbatim
87*> A is DOUBLE PRECISION array, dimension (LDA, N).
88*> On exit A N-by-N is initialized according to TYPE.
89*> \endverbatim
90*>
91*> \param[in] LDA
92*> \verbatim
93*> LDA is INTEGER
94*> The leading dimension of A and of B.
95*> \endverbatim
96*>
97*> \param[out] B
98*> \verbatim
99*> B is DOUBLE PRECISION array, dimension (LDA, N).
100*> On exit B N-by-N is initialized according to TYPE.
101*> \endverbatim
102*>
103*> \param[out] X
104*> \verbatim
105*> X is DOUBLE PRECISION array, dimension (LDX, N).
106*> On exit X is the N-by-N matrix of right eigenvectors.
107*> \endverbatim
108*>
109*> \param[in] LDX
110*> \verbatim
111*> LDX is INTEGER
112*> The leading dimension of X.
113*> \endverbatim
114*>
115*> \param[out] Y
116*> \verbatim
117*> Y is DOUBLE PRECISION array, dimension (LDY, N).
118*> On exit Y is the N-by-N matrix of left eigenvectors.
119*> \endverbatim
120*>
121*> \param[in] LDY
122*> \verbatim
123*> LDY is INTEGER
124*> The leading dimension of Y.
125*> \endverbatim
126*>
127*> \param[in] ALPHA
128*> \verbatim
129*> ALPHA is DOUBLE PRECISION
130*> \endverbatim
131*>
132*> \param[in] BETA
133*> \verbatim
134*> BETA is DOUBLE PRECISION
135*>
136*> Weighting constants for matrix A.
137*> \endverbatim
138*>
139*> \param[in] WX
140*> \verbatim
141*> WX is DOUBLE PRECISION
142*> Constant for right eigenvector matrix.
143*> \endverbatim
144*>
145*> \param[in] WY
146*> \verbatim
147*> WY is DOUBLE PRECISION
148*> Constant for left eigenvector matrix.
149*> \endverbatim
150*>
151*> \param[out] S
152*> \verbatim
153*> S is DOUBLE PRECISION array, dimension (N)
154*> S(i) is the reciprocal condition number for eigenvalue i.
155*> \endverbatim
156*>
157*> \param[out] DIF
158*> \verbatim
159*> DIF is DOUBLE PRECISION array, dimension (N)
160*> DIF(i) is the reciprocal condition number for eigenvector i.
161*> \endverbatim
162*
163* Authors:
164* ========
165*
166*> \author Univ. of Tennessee
167*> \author Univ. of California Berkeley
168*> \author Univ. of Colorado Denver
169*> \author NAG Ltd.
170*
171*> \ingroup double_matgen
172*
173* =====================================================================
174 SUBROUTINE dlatm6( TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA,
175 \$ BETA, WX, WY, S, DIF )
176*
177* -- LAPACK computational routine --
178* -- LAPACK is a software package provided by Univ. of Tennessee, --
179* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
180*
181* .. Scalar Arguments ..
182 INTEGER LDA, LDX, LDY, N, TYPE
183 DOUBLE PRECISION ALPHA, BETA, WX, WY
184* ..
185* .. Array Arguments ..
186 DOUBLE PRECISION A( LDA, * ), B( LDA, * ), DIF( * ), S( * ),
187 \$ x( ldx, * ), y( ldy, * )
188* ..
189*
190* =====================================================================
191*
192* .. Parameters ..
193 DOUBLE PRECISION ZERO, ONE, TWO, THREE
194 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
195 \$ three = 3.0d+0 )
196* ..
197* .. Local Scalars ..
198 INTEGER I, INFO, J
199* ..
200* .. Local Arrays ..
201 DOUBLE PRECISION WORK( 100 ), Z( 12, 12 )
202* ..
203* .. Intrinsic Functions ..
204 INTRINSIC dble, sqrt
205* ..
206* .. External Subroutines ..
207 EXTERNAL dgesvd, dlacpy, dlakf2
208* ..
209* .. Executable Statements ..
210*
211* Generate test problem ...
212* (Da, Db) ...
213*
214 DO 20 i = 1, n
215 DO 10 j = 1, n
216*
217 IF( i.EQ.j ) THEN
218 a( i, i ) = dble( i ) + alpha
219 b( i, i ) = one
220 ELSE
221 a( i, j ) = zero
222 b( i, j ) = zero
223 END IF
224*
225 10 CONTINUE
226 20 CONTINUE
227*
228* Form X and Y
229*
230 CALL dlacpy( 'F', n, n, b, lda, y, ldy )
231 y( 3, 1 ) = -wy
232 y( 4, 1 ) = wy
233 y( 5, 1 ) = -wy
234 y( 3, 2 ) = -wy
235 y( 4, 2 ) = wy
236 y( 5, 2 ) = -wy
237*
238 CALL dlacpy( 'F', n, n, b, lda, x, ldx )
239 x( 1, 3 ) = -wx
240 x( 1, 4 ) = -wx
241 x( 1, 5 ) = wx
242 x( 2, 3 ) = wx
243 x( 2, 4 ) = -wx
244 x( 2, 5 ) = -wx
245*
246* Form (A, B)
247*
248 b( 1, 3 ) = wx + wy
249 b( 2, 3 ) = -wx + wy
250 b( 1, 4 ) = wx - wy
251 b( 2, 4 ) = wx - wy
252 b( 1, 5 ) = -wx + wy
253 b( 2, 5 ) = wx + wy
254 IF( type.EQ.1 ) THEN
255 a( 1, 3 ) = wx*a( 1, 1 ) + wy*a( 3, 3 )
256 a( 2, 3 ) = -wx*a( 2, 2 ) + wy*a( 3, 3 )
257 a( 1, 4 ) = wx*a( 1, 1 ) - wy*a( 4, 4 )
258 a( 2, 4 ) = wx*a( 2, 2 ) - wy*a( 4, 4 )
259 a( 1, 5 ) = -wx*a( 1, 1 ) + wy*a( 5, 5 )
260 a( 2, 5 ) = wx*a( 2, 2 ) + wy*a( 5, 5 )
261 ELSE IF( type.EQ.2 ) THEN
262 a( 1, 3 ) = two*wx + wy
263 a( 2, 3 ) = wy
264 a( 1, 4 ) = -wy*( two+alpha+beta )
265 a( 2, 4 ) = two*wx - wy*( two+alpha+beta )
266 a( 1, 5 ) = -two*wx + wy*( alpha-beta )
267 a( 2, 5 ) = wy*( alpha-beta )
268 a( 1, 1 ) = one
269 a( 1, 2 ) = -one
270 a( 2, 1 ) = one
271 a( 2, 2 ) = a( 1, 1 )
272 a( 3, 3 ) = one
273 a( 4, 4 ) = one + alpha
274 a( 4, 5 ) = one + beta
275 a( 5, 4 ) = -a( 4, 5 )
276 a( 5, 5 ) = a( 4, 4 )
277 END IF
278*
279* Compute condition numbers
280*
281 IF( type.EQ.1 ) THEN
282*
283 s( 1 ) = one / sqrt( ( one+three*wy*wy ) /
284 \$ ( one+a( 1, 1 )*a( 1, 1 ) ) )
285 s( 2 ) = one / sqrt( ( one+three*wy*wy ) /
286 \$ ( one+a( 2, 2 )*a( 2, 2 ) ) )
287 s( 3 ) = one / sqrt( ( one+two*wx*wx ) /
288 \$ ( one+a( 3, 3 )*a( 3, 3 ) ) )
289 s( 4 ) = one / sqrt( ( one+two*wx*wx ) /
290 \$ ( one+a( 4, 4 )*a( 4, 4 ) ) )
291 s( 5 ) = one / sqrt( ( one+two*wx*wx ) /
292 \$ ( one+a( 5, 5 )*a( 5, 5 ) ) )
293*
294 CALL dlakf2( 1, 4, a, lda, a( 2, 2 ), b, b( 2, 2 ), z, 12 )
295 CALL dgesvd( 'N', 'N', 8, 8, z, 12, work, work( 9 ), 1,
296 \$ work( 10 ), 1, work( 11 ), 40, info )
297 dif( 1 ) = work( 8 )
298*
299 CALL dlakf2( 4, 1, a, lda, a( 5, 5 ), b, b( 5, 5 ), z, 12 )
300 CALL dgesvd( 'N', 'N', 8, 8, z, 12, work, work( 9 ), 1,
301 \$ work( 10 ), 1, work( 11 ), 40, info )
302 dif( 5 ) = work( 8 )
303*
304 ELSE IF( type.EQ.2 ) THEN
305*
306 s( 1 ) = one / sqrt( one / three+wy*wy )
307 s( 2 ) = s( 1 )
308 s( 3 ) = one / sqrt( one / two+wx*wx )
309 s( 4 ) = one / sqrt( ( one+two*wx*wx ) /
310 \$ ( one+( one+alpha )*( one+alpha )+( one+beta )*( one+
311 \$ beta ) ) )
312 s( 5 ) = s( 4 )
313*
314 CALL dlakf2( 2, 3, a, lda, a( 3, 3 ), b, b( 3, 3 ), z, 12 )
315 CALL dgesvd( 'N', 'N', 12, 12, z, 12, work, work( 13 ), 1,
316 \$ work( 14 ), 1, work( 15 ), 60, info )
317 dif( 1 ) = work( 12 )
318*
319 CALL dlakf2( 3, 2, a, lda, a( 4, 4 ), b, b( 4, 4 ), z, 12 )
320 CALL dgesvd( 'N', 'N', 12, 12, z, 12, work, work( 13 ), 1,
321 \$ work( 14 ), 1, work( 15 ), 60, info )
322 dif( 5 ) = work( 12 )
323*
324 END IF
325*
326 RETURN
327*
328* End of DLATM6
329*
330 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 dlakf2(M, N, A, LDA, B, D, E, Z, LDZ)
DLAKF2
Definition: dlakf2.f:105
subroutine dlatm6(TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA, BETA, WX, WY, S, DIF)
DLATM6
Definition: dlatm6.f:176
subroutine dgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO)
DGESVD computes the singular value decomposition (SVD) for GE matrices
Definition: dgesvd.f:211