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