LAPACK  3.10.0
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 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