LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 futher 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 *> \date November 2011
172 *
173 *> \ingroup double_matgen
174 *
175 * =====================================================================
176  SUBROUTINE dlatm6( TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA,
177  $ beta, wx, wy, s, dif )
178 *
179 * -- LAPACK computational routine (version 3.4.0) --
180 * -- LAPACK is a software package provided by Univ. of Tennessee, --
181 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
182 * November 2011
183 *
184 * .. Scalar Arguments ..
185  INTEGER LDA, LDX, LDY, N, TYPE
186  DOUBLE PRECISION ALPHA, BETA, WX, WY
187 * ..
188 * .. Array Arguments ..
189  DOUBLE PRECISION A( lda, * ), B( lda, * ), DIF( * ), S( * ),
190  $ x( ldx, * ), y( ldy, * )
191 * ..
192 *
193 * =====================================================================
194 *
195 * .. Parameters ..
196  DOUBLE PRECISION ZERO, ONE, TWO, THREE
197  parameter ( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
198  $ three = 3.0d+0 )
199 * ..
200 * .. Local Scalars ..
201  INTEGER I, INFO, J
202 * ..
203 * .. Local Arrays ..
204  DOUBLE PRECISION WORK( 100 ), Z( 12, 12 )
205 * ..
206 * .. Intrinsic Functions ..
207  INTRINSIC dble, sqrt
208 * ..
209 * .. External Subroutines ..
210  EXTERNAL dgesvd, dlacpy, dlakf2
211 * ..
212 * .. Executable Statements ..
213 *
214 * Generate test problem ...
215 * (Da, Db) ...
216 *
217  DO 20 i = 1, n
218  DO 10 j = 1, n
219 *
220  IF( i.EQ.j ) THEN
221  a( i, i ) = dble( i ) + alpha
222  b( i, i ) = one
223  ELSE
224  a( i, j ) = zero
225  b( i, j ) = zero
226  END IF
227 *
228  10 CONTINUE
229  20 CONTINUE
230 *
231 * Form X and Y
232 *
233  CALL dlacpy( 'F', n, n, b, lda, y, ldy )
234  y( 3, 1 ) = -wy
235  y( 4, 1 ) = wy
236  y( 5, 1 ) = -wy
237  y( 3, 2 ) = -wy
238  y( 4, 2 ) = wy
239  y( 5, 2 ) = -wy
240 *
241  CALL dlacpy( 'F', n, n, b, lda, x, ldx )
242  x( 1, 3 ) = -wx
243  x( 1, 4 ) = -wx
244  x( 1, 5 ) = wx
245  x( 2, 3 ) = wx
246  x( 2, 4 ) = -wx
247  x( 2, 5 ) = -wx
248 *
249 * Form (A, B)
250 *
251  b( 1, 3 ) = wx + wy
252  b( 2, 3 ) = -wx + wy
253  b( 1, 4 ) = wx - wy
254  b( 2, 4 ) = wx - wy
255  b( 1, 5 ) = -wx + wy
256  b( 2, 5 ) = wx + wy
257  IF( type.EQ.1 ) THEN
258  a( 1, 3 ) = wx*a( 1, 1 ) + wy*a( 3, 3 )
259  a( 2, 3 ) = -wx*a( 2, 2 ) + wy*a( 3, 3 )
260  a( 1, 4 ) = wx*a( 1, 1 ) - wy*a( 4, 4 )
261  a( 2, 4 ) = wx*a( 2, 2 ) - wy*a( 4, 4 )
262  a( 1, 5 ) = -wx*a( 1, 1 ) + wy*a( 5, 5 )
263  a( 2, 5 ) = wx*a( 2, 2 ) + wy*a( 5, 5 )
264  ELSE IF( type.EQ.2 ) THEN
265  a( 1, 3 ) = two*wx + wy
266  a( 2, 3 ) = wy
267  a( 1, 4 ) = -wy*( two+alpha+beta )
268  a( 2, 4 ) = two*wx - wy*( two+alpha+beta )
269  a( 1, 5 ) = -two*wx + wy*( alpha-beta )
270  a( 2, 5 ) = wy*( alpha-beta )
271  a( 1, 1 ) = one
272  a( 1, 2 ) = -one
273  a( 2, 1 ) = one
274  a( 2, 2 ) = a( 1, 1 )
275  a( 3, 3 ) = one
276  a( 4, 4 ) = one + alpha
277  a( 4, 5 ) = one + beta
278  a( 5, 4 ) = -a( 4, 5 )
279  a( 5, 5 ) = a( 4, 4 )
280  END IF
281 *
282 * Compute condition numbers
283 *
284  IF( type.EQ.1 ) THEN
285 *
286  s( 1 ) = one / sqrt( ( one+three*wy*wy ) /
287  $ ( one+a( 1, 1 )*a( 1, 1 ) ) )
288  s( 2 ) = one / sqrt( ( one+three*wy*wy ) /
289  $ ( one+a( 2, 2 )*a( 2, 2 ) ) )
290  s( 3 ) = one / sqrt( ( one+two*wx*wx ) /
291  $ ( one+a( 3, 3 )*a( 3, 3 ) ) )
292  s( 4 ) = one / sqrt( ( one+two*wx*wx ) /
293  $ ( one+a( 4, 4 )*a( 4, 4 ) ) )
294  s( 5 ) = one / sqrt( ( one+two*wx*wx ) /
295  $ ( one+a( 5, 5 )*a( 5, 5 ) ) )
296 *
297  CALL dlakf2( 1, 4, a, lda, a( 2, 2 ), b, b( 2, 2 ), z, 12 )
298  CALL dgesvd( 'N', 'N', 8, 8, z, 12, work, work( 9 ), 1,
299  $ work( 10 ), 1, work( 11 ), 40, info )
300  dif( 1 ) = work( 8 )
301 *
302  CALL dlakf2( 4, 1, a, lda, a( 5, 5 ), b, b( 5, 5 ), z, 12 )
303  CALL dgesvd( 'N', 'N', 8, 8, z, 12, work, work( 9 ), 1,
304  $ work( 10 ), 1, work( 11 ), 40, info )
305  dif( 5 ) = work( 8 )
306 *
307  ELSE IF( type.EQ.2 ) THEN
308 *
309  s( 1 ) = one / sqrt( one / three+wy*wy )
310  s( 2 ) = s( 1 )
311  s( 3 ) = one / sqrt( one / two+wx*wx )
312  s( 4 ) = one / sqrt( ( one+two*wx*wx ) /
313  $ ( one+( one+alpha )*( one+alpha )+( one+beta )*( one+
314  $ beta ) ) )
315  s( 5 ) = s( 4 )
316 *
317  CALL dlakf2( 2, 3, a, lda, a( 3, 3 ), b, b( 3, 3 ), z, 12 )
318  CALL dgesvd( 'N', 'N', 12, 12, z, 12, work, work( 13 ), 1,
319  $ work( 14 ), 1, work( 15 ), 60, info )
320  dif( 1 ) = work( 12 )
321 *
322  CALL dlakf2( 3, 2, a, lda, a( 4, 4 ), b, b( 4, 4 ), z, 12 )
323  CALL dgesvd( 'N', 'N', 12, 12, z, 12, work, work( 13 ), 1,
324  $ work( 14 ), 1, work( 15 ), 60, info )
325  dif( 5 ) = work( 12 )
326 *
327  END IF
328 *
329  RETURN
330 *
331 * End of DLATM6
332 *
333  END
subroutine dlakf2(M, N, A, LDA, B, D, E, Z, LDZ)
DLAKF2
Definition: dlakf2.f:107
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
subroutine dlatm6(TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA, BETA, WX, WY, S, DIF)
DLATM6
Definition: dlatm6.f:178
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:213