LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
clatm6.f
Go to the documentation of this file.
1 *> \brief \b CLATM6
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 CLATM6( 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 ALPHA, BETA, WX, WY
17 * ..
18 * .. Array Arguments ..
19 * REAL DIF( * ), S( * )
20 * COMPLEX A( LDA, * ), B( LDA, * ), X( LDX, * ),
21 * $ Y( LDY, * )
22 * ..
23 *
24 *
25 *> \par Purpose:
26 * =============
27 *>
28 *> \verbatim
29 *>
30 *> CLATM6 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 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 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 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 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
128 *> \endverbatim
129 *>
130 *> \param[in] BETA
131 *> \verbatim
132 *> BETA is COMPLEX
133 *>
134 *> Weighting constants for matrix A.
135 *> \endverbatim
136 *>
137 *> \param[in] WX
138 *> \verbatim
139 *> WX is COMPLEX
140 *> Constant for right eigenvector matrix.
141 *> \endverbatim
142 *>
143 *> \param[in] WY
144 *> \verbatim
145 *> WY is COMPLEX
146 *> Constant for left eigenvector matrix.
147 *> \endverbatim
148 *>
149 *> \param[out] S
150 *> \verbatim
151 *> S is REAL 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 REAL 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 *> \date December 2016
170 *
171 *> \ingroup complex_matgen
172 *
173 * =====================================================================
174  SUBROUTINE clatm6( TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA,
175  $ BETA, WX, WY, S, DIF )
176 *
177 * -- LAPACK computational routine (version 3.7.0) --
178 * -- LAPACK is a software package provided by Univ. of Tennessee, --
179 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
180 * December 2016
181 *
182 * .. Scalar Arguments ..
183  INTEGER LDA, LDX, LDY, N, TYPE
184  COMPLEX ALPHA, BETA, WX, WY
185 * ..
186 * .. Array Arguments ..
187  REAL DIF( * ), S( * )
188  COMPLEX A( LDA, * ), B( LDA, * ), X( LDX, * ),
189  $ y( ldy, * )
190 * ..
191 *
192 * =====================================================================
193 *
194 * .. Parameters ..
195  REAL RONE, TWO, THREE
196  parameter( rone = 1.0e+0, two = 2.0e+0, three = 3.0e+0 )
197  COMPLEX ZERO, ONE
198  parameter( zero = ( 0.0e+0, 0.0e+0 ),
199  $ one = ( 1.0e+0, 0.0e+0 ) )
200 * ..
201 * .. Local Scalars ..
202  INTEGER I, INFO, J
203 * ..
204 * .. Local Arrays ..
205  REAL RWORK( 50 )
206  COMPLEX WORK( 26 ), Z( 8, 8 )
207 * ..
208 * .. Intrinsic Functions ..
209  INTRINSIC cabs, cmplx, conjg, real, sqrt
210 * ..
211 * .. External Subroutines ..
212  EXTERNAL cgesvd, clacpy, clakf2
213 * ..
214 * .. Executable Statements ..
215 *
216 * Generate test problem ...
217 * (Da, Db) ...
218 *
219  DO 20 i = 1, n
220  DO 10 j = 1, n
221 *
222  IF( i.EQ.j ) THEN
223  a( i, i ) = cmplx( i ) + alpha
224  b( i, i ) = one
225  ELSE
226  a( i, j ) = zero
227  b( i, j ) = zero
228  END IF
229 *
230  10 CONTINUE
231  20 CONTINUE
232  IF( type.EQ.2 ) THEN
233  a( 1, 1 ) = cmplx( rone, rone )
234  a( 2, 2 ) = conjg( a( 1, 1 ) )
235  a( 3, 3 ) = one
236  a( 4, 4 ) = cmplx( real( one+alpha ), real( one+beta ) )
237  a( 5, 5 ) = conjg( a( 4, 4 ) )
238  END IF
239 *
240 * Form X and Y
241 *
242  CALL clacpy( 'F', n, n, b, lda, y, ldy )
243  y( 3, 1 ) = -conjg( wy )
244  y( 4, 1 ) = conjg( wy )
245  y( 5, 1 ) = -conjg( wy )
246  y( 3, 2 ) = -conjg( wy )
247  y( 4, 2 ) = conjg( wy )
248  y( 5, 2 ) = -conjg( wy )
249 *
250  CALL clacpy( 'F', n, n, b, lda, x, ldx )
251  x( 1, 3 ) = -wx
252  x( 1, 4 ) = -wx
253  x( 1, 5 ) = wx
254  x( 2, 3 ) = wx
255  x( 2, 4 ) = -wx
256  x( 2, 5 ) = -wx
257 *
258 * Form (A, B)
259 *
260  b( 1, 3 ) = wx + wy
261  b( 2, 3 ) = -wx + wy
262  b( 1, 4 ) = wx - wy
263  b( 2, 4 ) = wx - wy
264  b( 1, 5 ) = -wx + wy
265  b( 2, 5 ) = wx + wy
266  a( 1, 3 ) = wx*a( 1, 1 ) + wy*a( 3, 3 )
267  a( 2, 3 ) = -wx*a( 2, 2 ) + wy*a( 3, 3 )
268  a( 1, 4 ) = wx*a( 1, 1 ) - wy*a( 4, 4 )
269  a( 2, 4 ) = wx*a( 2, 2 ) - wy*a( 4, 4 )
270  a( 1, 5 ) = -wx*a( 1, 1 ) + wy*a( 5, 5 )
271  a( 2, 5 ) = wx*a( 2, 2 ) + wy*a( 5, 5 )
272 *
273 * Compute condition numbers
274 *
275  s( 1 ) = rone / sqrt( ( rone+three*cabs( wy )*cabs( wy ) ) /
276  $ ( rone+cabs( a( 1, 1 ) )*cabs( a( 1, 1 ) ) ) )
277  s( 2 ) = rone / sqrt( ( rone+three*cabs( wy )*cabs( wy ) ) /
278  $ ( rone+cabs( a( 2, 2 ) )*cabs( a( 2, 2 ) ) ) )
279  s( 3 ) = rone / sqrt( ( rone+two*cabs( wx )*cabs( wx ) ) /
280  $ ( rone+cabs( a( 3, 3 ) )*cabs( a( 3, 3 ) ) ) )
281  s( 4 ) = rone / sqrt( ( rone+two*cabs( wx )*cabs( wx ) ) /
282  $ ( rone+cabs( a( 4, 4 ) )*cabs( a( 4, 4 ) ) ) )
283  s( 5 ) = rone / sqrt( ( rone+two*cabs( wx )*cabs( wx ) ) /
284  $ ( rone+cabs( a( 5, 5 ) )*cabs( a( 5, 5 ) ) ) )
285 *
286  CALL clakf2( 1, 4, a, lda, a( 2, 2 ), b, b( 2, 2 ), z, 8 )
287  CALL cgesvd( 'N', 'N', 8, 8, z, 8, rwork, work, 1, work( 2 ), 1,
288  $ work( 3 ), 24, rwork( 9 ), info )
289  dif( 1 ) = rwork( 8 )
290 *
291  CALL clakf2( 4, 1, a, lda, a( 5, 5 ), b, b( 5, 5 ), z, 8 )
292  CALL cgesvd( 'N', 'N', 8, 8, z, 8, rwork, work, 1, work( 2 ), 1,
293  $ work( 3 ), 24, rwork( 9 ), info )
294  dif( 5 ) = rwork( 8 )
295 *
296  RETURN
297 *
298 * End of CLATM6
299 *
300  END
clacpy
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
clakf2
subroutine clakf2(M, N, A, LDA, B, D, E, Z, LDZ)
CLAKF2
Definition: clakf2.f:107
cgesvd
subroutine cgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO)
CGESVD computes the singular value decomposition (SVD) for GE matrices
Definition: cgesvd.f:216
clatm6
subroutine clatm6(TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA, BETA, WX, WY, S, DIF)
CLATM6
Definition: clatm6.f:176