1 /* ../netlib/claein.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib;
2 on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */
3 #include "FLA_f2c.h" /* Table of constant values */
4 static integer c__1 = 1;
5 /* > \brief \b CLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration. */
6 /* =========== DOCUMENTATION =========== */
7 /* Online html documentation available at */
8 /* http://www.netlib.org/lapack/explore-html/ */
9 /* > \htmlonly */
10 /* > Download CLAEIN + dependencies */
11 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claein. f"> */
12 /* > [TGZ]</a> */
13 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claein. f"> */
14 /* > [ZIP]</a> */
15 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claein. f"> */
16 /* > [TXT]</a> */
17 /* > \endhtmlonly */
18 /* Definition: */
19 /* =========== */
20 /* SUBROUTINE CLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK, */
21 /* EPS3, SMLNUM, INFO ) */
22 /* .. Scalar Arguments .. */
23 /* LOGICAL NOINIT, RIGHTV */
24 /* INTEGER INFO, LDB, LDH, N */
25 /* REAL EPS3, SMLNUM */
26 /* COMPLEX W */
27 /* .. */
28 /* .. Array Arguments .. */
29 /* REAL RWORK( * ) */
30 /* COMPLEX B( LDB, * ), H( LDH, * ), V( * ) */
31 /* .. */
32 /* > \par Purpose: */
33 /* ============= */
34 /* > */
35 /* > \verbatim */
36 /* > */
37 /* > CLAEIN uses inverse iteration to find a right or left eigenvector */
38 /* > corresponding to the eigenvalue W of a complex upper Hessenberg */
39 /* > matrix H. */
40 /* > \endverbatim */
41 /* Arguments: */
42 /* ========== */
43 /* > \param[in] RIGHTV */
44 /* > \verbatim */
45 /* > RIGHTV is LOGICAL */
46 /* > = .TRUE. : compute right eigenvector;
47 */
48 /* > = .FALSE.: compute left eigenvector. */
49 /* > \endverbatim */
50 /* > */
51 /* > \param[in] NOINIT */
52 /* > \verbatim */
53 /* > NOINIT is LOGICAL */
54 /* > = .TRUE. : no initial vector supplied in V */
55 /* > = .FALSE.: initial vector supplied in V. */
56 /* > \endverbatim */
57 /* > */
58 /* > \param[in] N */
59 /* > \verbatim */
60 /* > N is INTEGER */
61 /* > The order of the matrix H. N >= 0. */
62 /* > \endverbatim */
63 /* > */
64 /* > \param[in] H */
65 /* > \verbatim */
66 /* > H is COMPLEX array, dimension (LDH,N) */
67 /* > The upper Hessenberg matrix H. */
68 /* > \endverbatim */
69 /* > */
70 /* > \param[in] LDH */
71 /* > \verbatim */
72 /* > LDH is INTEGER */
73 /* > The leading dimension of the array H. LDH >= max(1,N). */
74 /* > \endverbatim */
75 /* > */
76 /* > \param[in] W */
77 /* > \verbatim */
78 /* > W is COMPLEX */
79 /* > The eigenvalue of H whose corresponding right or left */
80 /* > eigenvector is to be computed. */
81 /* > \endverbatim */
82 /* > */
83 /* > \param[in,out] V */
84 /* > \verbatim */
85 /* > V is COMPLEX array, dimension (N) */
86 /* > On entry, if NOINIT = .FALSE., V must contain a starting */
87 /* > vector for inverse iteration;
88 otherwise V need not be set. */
89 /* > On exit, V contains the computed eigenvector, normalized so */
90 /* > that the component of largest magnitude has magnitude 1;
91 here */
92 /* > the magnitude of a complex number (x,y) is taken to be */
93 /* > |x| + |y|. */
94 /* > \endverbatim */
95 /* > */
96 /* > \param[out] B */
97 /* > \verbatim */
98 /* > B is COMPLEX array, dimension (LDB,N) */
99 /* > \endverbatim */
100 /* > */
101 /* > \param[in] LDB */
102 /* > \verbatim */
103 /* > LDB is INTEGER */
104 /* > The leading dimension of the array B. LDB >= max(1,N). */
105 /* > \endverbatim */
106 /* > */
107 /* > \param[out] RWORK */
108 /* > \verbatim */
109 /* > RWORK is REAL array, dimension (N) */
110 /* > \endverbatim */
111 /* > */
112 /* > \param[in] EPS3 */
113 /* > \verbatim */
114 /* > EPS3 is REAL */
115 /* > A small machine-dependent value which is used to perturb */
116 /* > close eigenvalues, and to replace zero pivots. */
117 /* > \endverbatim */
118 /* > */
119 /* > \param[in] SMLNUM */
120 /* > \verbatim */
121 /* > SMLNUM is REAL */
122 /* > A machine-dependent value close to the underflow threshold. */
123 /* > \endverbatim */
124 /* > */
125 /* > \param[out] INFO */
126 /* > \verbatim */
127 /* > INFO is INTEGER */
128 /* > = 0: successful exit */
129 /* > = 1: inverse iteration did not converge;
130 V is set to the */
131 /* > last iterate. */
132 /* > \endverbatim */
133 /* Authors: */
134 /* ======== */
135 /* > \author Univ. of Tennessee */
136 /* > \author Univ. of California Berkeley */
137 /* > \author Univ. of Colorado Denver */
138 /* > \author NAG Ltd. */
139 /* > \date September 2012 */
140 /* > \ingroup complexOTHERauxiliary */
141 /* ===================================================================== */
142 /* Subroutine */
claein_(logical * rightv,logical * noinit,integer * n,complex * h__,integer * ldh,complex * w,complex * v,complex * b,integer * ldb,real * rwork,real * eps3,real * smlnum,integer * info)143 int claein_(logical *rightv, logical *noinit, integer *n, complex *h__, integer *ldh, complex *w, complex *v, complex *b, integer *ldb, real *rwork, real *eps3, real *smlnum, integer *info)
144 {
145 /* System generated locals */
146 integer b_dim1, b_offset, h_dim1, h_offset, i__1, i__2, i__3, i__4, i__5;
147 real r__1, r__2, r__3, r__4;
148 complex q__1, q__2;
149 /* Builtin functions */
150 double sqrt(doublereal), r_imag(complex *);
151 /* Local variables */
152 integer i__, j;
153 complex x, ei, ej;
154 integer its, ierr;
155 complex temp;
156 real scale;
157 char trans[1];
158 real rtemp, rootn, vnorm;
159 extern real scnrm2_(integer *, complex *, integer *);
160 extern integer icamax_(integer *, complex *, integer *);
161 extern /* Complex */
162 VOID cladiv_(complex *, complex *, complex *);
163 extern /* Subroutine */
164 int csscal_(integer *, real *, complex *, integer *), clatrs_(char *, char *, char *, char *, integer *, complex *, integer *, complex *, real *, real *, integer *);
165 extern real scasum_(integer *, complex *, integer *);
166 char normin[1];
167 real nrmsml, growto;
168 /* -- LAPACK auxiliary routine (version 3.4.2) -- */
169 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
170 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
171 /* September 2012 */
172 /* .. Scalar Arguments .. */
173 /* .. */
174 /* .. Array Arguments .. */
175 /* .. */
176 /* ===================================================================== */
177 /* .. Parameters .. */
178 /* .. */
179 /* .. Local Scalars .. */
180 /* .. */
181 /* .. External Functions .. */
182 /* .. */
183 /* .. External Subroutines .. */
184 /* .. */
185 /* .. Intrinsic Functions .. */
186 /* .. */
187 /* .. Statement Functions .. */
188 /* .. */
189 /* .. Statement Function definitions .. */
190 /* .. */
191 /* .. Executable Statements .. */
192 /* Parameter adjustments */
193 h_dim1 = *ldh;
194 h_offset = 1 + h_dim1;
195 h__ -= h_offset;
196 --v;
197 b_dim1 = *ldb;
198 b_offset = 1 + b_dim1;
199 b -= b_offset;
200 --rwork;
201 /* Function Body */
202 *info = 0;
203 /* GROWTO is the threshold used in the acceptance test for an */
204 /* eigenvector. */
205 rootn = sqrt((real) (*n));
206 growto = .1f / rootn;
207 /* Computing MAX */
208 r__1 = 1.f;
209 r__2 = *eps3 * rootn; // , expr subst
210 nrmsml = max(r__1,r__2) * *smlnum;
211 /* Form B = H - W*I (except that the subdiagonal elements are not */
212 /* stored). */
213 i__1 = *n;
214 for (j = 1;
215 j <= i__1;
216 ++j)
217 {
218 i__2 = j - 1;
219 for (i__ = 1;
220 i__ <= i__2;
221 ++i__)
222 {
223 i__3 = i__ + j * b_dim1;
224 i__4 = i__ + j * h_dim1;
225 b[i__3].r = h__[i__4].r;
226 b[i__3].i = h__[i__4].i; // , expr subst
227 /* L10: */
228 }
229 i__2 = j + j * b_dim1;
230 i__3 = j + j * h_dim1;
231 q__1.r = h__[i__3].r - w->r;
232 q__1.i = h__[i__3].i - w->i; // , expr subst
233 b[i__2].r = q__1.r;
234 b[i__2].i = q__1.i; // , expr subst
235 /* L20: */
236 }
237 if (*noinit)
238 {
239 /* Initialize V. */
240 i__1 = *n;
241 for (i__ = 1;
242 i__ <= i__1;
243 ++i__)
244 {
245 i__2 = i__;
246 v[i__2].r = *eps3;
247 v[i__2].i = 0.f; // , expr subst
248 /* L30: */
249 }
250 }
251 else
252 {
253 /* Scale supplied initial vector. */
254 vnorm = scnrm2_(n, &v[1], &c__1);
255 r__1 = *eps3 * rootn / max(vnorm,nrmsml);
256 csscal_(n, &r__1, &v[1], &c__1);
257 }
258 if (*rightv)
259 {
260 /* LU decomposition with partial pivoting of B, replacing zero */
261 /* pivots by EPS3. */
262 i__1 = *n - 1;
263 for (i__ = 1;
264 i__ <= i__1;
265 ++i__)
266 {
267 i__2 = i__ + 1 + i__ * h_dim1;
268 ei.r = h__[i__2].r;
269 ei.i = h__[i__2].i; // , expr subst
270 i__2 = i__ + i__ * b_dim1;
271 if ((r__1 = b[i__2].r, f2c_abs(r__1)) + (r__2 = r_imag(&b[i__ + i__ * b_dim1]), f2c_abs(r__2)) < (r__3 = ei.r, f2c_abs(r__3)) + (r__4 = r_imag(&ei), f2c_abs(r__4)))
272 {
273 /* Interchange rows and eliminate. */
274 cladiv_(&q__1, &b[i__ + i__ * b_dim1], &ei);
275 x.r = q__1.r;
276 x.i = q__1.i; // , expr subst
277 i__2 = i__ + i__ * b_dim1;
278 b[i__2].r = ei.r;
279 b[i__2].i = ei.i; // , expr subst
280 i__2 = *n;
281 for (j = i__ + 1;
282 j <= i__2;
283 ++j)
284 {
285 i__3 = i__ + 1 + j * b_dim1;
286 temp.r = b[i__3].r;
287 temp.i = b[i__3].i; // , expr subst
288 i__3 = i__ + 1 + j * b_dim1;
289 i__4 = i__ + j * b_dim1;
290 q__2.r = x.r * temp.r - x.i * temp.i;
291 q__2.i = x.r * temp.i + x.i * temp.r; // , expr subst
292 q__1.r = b[i__4].r - q__2.r;
293 q__1.i = b[i__4].i - q__2.i; // , expr subst
294 b[i__3].r = q__1.r;
295 b[i__3].i = q__1.i; // , expr subst
296 i__3 = i__ + j * b_dim1;
297 b[i__3].r = temp.r;
298 b[i__3].i = temp.i; // , expr subst
299 /* L40: */
300 }
301 }
302 else
303 {
304 /* Eliminate without interchange. */
305 i__2 = i__ + i__ * b_dim1;
306 if (b[i__2].r == 0.f && b[i__2].i == 0.f)
307 {
308 i__3 = i__ + i__ * b_dim1;
309 b[i__3].r = *eps3;
310 b[i__3].i = 0.f; // , expr subst
311 }
312 cladiv_(&q__1, &ei, &b[i__ + i__ * b_dim1]);
313 x.r = q__1.r;
314 x.i = q__1.i; // , expr subst
315 if (x.r != 0.f || x.i != 0.f)
316 {
317 i__2 = *n;
318 for (j = i__ + 1;
319 j <= i__2;
320 ++j)
321 {
322 i__3 = i__ + 1 + j * b_dim1;
323 i__4 = i__ + 1 + j * b_dim1;
324 i__5 = i__ + j * b_dim1;
325 q__2.r = x.r * b[i__5].r - x.i * b[i__5].i;
326 q__2.i = x.r * b[i__5].i + x.i * b[i__5].r; // , expr subst
327 q__1.r = b[i__4].r - q__2.r;
328 q__1.i = b[i__4].i - q__2.i; // , expr subst
329 b[i__3].r = q__1.r;
330 b[i__3].i = q__1.i; // , expr subst
331 /* L50: */
332 }
333 }
334 }
335 /* L60: */
336 }
337 i__1 = *n + *n * b_dim1;
338 if (b[i__1].r == 0.f && b[i__1].i == 0.f)
339 {
340 i__2 = *n + *n * b_dim1;
341 b[i__2].r = *eps3;
342 b[i__2].i = 0.f; // , expr subst
343 }
344 *(unsigned char *)trans = 'N';
345 }
346 else
347 {
348 /* UL decomposition with partial pivoting of B, replacing zero */
349 /* pivots by EPS3. */
350 for (j = *n;
351 j >= 2;
352 --j)
353 {
354 i__1 = j + (j - 1) * h_dim1;
355 ej.r = h__[i__1].r;
356 ej.i = h__[i__1].i; // , expr subst
357 i__1 = j + j * b_dim1;
358 if ((r__1 = b[i__1].r, f2c_abs(r__1)) + (r__2 = r_imag(&b[j + j * b_dim1]), f2c_abs(r__2)) < (r__3 = ej.r, f2c_abs(r__3)) + (r__4 = r_imag(&ej), f2c_abs(r__4)))
359 {
360 /* Interchange columns and eliminate. */
361 cladiv_(&q__1, &b[j + j * b_dim1], &ej);
362 x.r = q__1.r;
363 x.i = q__1.i; // , expr subst
364 i__1 = j + j * b_dim1;
365 b[i__1].r = ej.r;
366 b[i__1].i = ej.i; // , expr subst
367 i__1 = j - 1;
368 for (i__ = 1;
369 i__ <= i__1;
370 ++i__)
371 {
372 i__2 = i__ + (j - 1) * b_dim1;
373 temp.r = b[i__2].r;
374 temp.i = b[i__2].i; // , expr subst
375 i__2 = i__ + (j - 1) * b_dim1;
376 i__3 = i__ + j * b_dim1;
377 q__2.r = x.r * temp.r - x.i * temp.i;
378 q__2.i = x.r * temp.i + x.i * temp.r; // , expr subst
379 q__1.r = b[i__3].r - q__2.r;
380 q__1.i = b[i__3].i - q__2.i; // , expr subst
381 b[i__2].r = q__1.r;
382 b[i__2].i = q__1.i; // , expr subst
383 i__2 = i__ + j * b_dim1;
384 b[i__2].r = temp.r;
385 b[i__2].i = temp.i; // , expr subst
386 /* L70: */
387 }
388 }
389 else
390 {
391 /* Eliminate without interchange. */
392 i__1 = j + j * b_dim1;
393 if (b[i__1].r == 0.f && b[i__1].i == 0.f)
394 {
395 i__2 = j + j * b_dim1;
396 b[i__2].r = *eps3;
397 b[i__2].i = 0.f; // , expr subst
398 }
399 cladiv_(&q__1, &ej, &b[j + j * b_dim1]);
400 x.r = q__1.r;
401 x.i = q__1.i; // , expr subst
402 if (x.r != 0.f || x.i != 0.f)
403 {
404 i__1 = j - 1;
405 for (i__ = 1;
406 i__ <= i__1;
407 ++i__)
408 {
409 i__2 = i__ + (j - 1) * b_dim1;
410 i__3 = i__ + (j - 1) * b_dim1;
411 i__4 = i__ + j * b_dim1;
412 q__2.r = x.r * b[i__4].r - x.i * b[i__4].i;
413 q__2.i = x.r * b[i__4].i + x.i * b[i__4].r; // , expr subst
414 q__1.r = b[i__3].r - q__2.r;
415 q__1.i = b[i__3].i - q__2.i; // , expr subst
416 b[i__2].r = q__1.r;
417 b[i__2].i = q__1.i; // , expr subst
418 /* L80: */
419 }
420 }
421 }
422 /* L90: */
423 }
424 i__1 = b_dim1 + 1;
425 if (b[i__1].r == 0.f && b[i__1].i == 0.f)
426 {
427 i__2 = b_dim1 + 1;
428 b[i__2].r = *eps3;
429 b[i__2].i = 0.f; // , expr subst
430 }
431 *(unsigned char *)trans = 'C';
432 }
433 *(unsigned char *)normin = 'N';
434 i__1 = *n;
435 for (its = 1;
436 its <= i__1;
437 ++its)
438 {
439 /* Solve U*x = scale*v for a right eigenvector */
440 /* or U**H *x = scale*v for a left eigenvector, */
441 /* overwriting x on v. */
442 clatrs_("Upper", trans, "Nonunit", normin, n, &b[b_offset], ldb, &v[1] , &scale, &rwork[1], &ierr);
443 *(unsigned char *)normin = 'Y';
444 /* Test for sufficient growth in the norm of v. */
445 vnorm = scasum_(n, &v[1], &c__1);
446 if (vnorm >= growto * scale)
447 {
448 goto L120;
449 }
450 /* Choose new orthogonal starting vector and try again. */
451 rtemp = *eps3 / (rootn + 1.f);
452 v[1].r = *eps3;
453 v[1].i = 0.f; // , expr subst
454 i__2 = *n;
455 for (i__ = 2;
456 i__ <= i__2;
457 ++i__)
458 {
459 i__3 = i__;
460 v[i__3].r = rtemp;
461 v[i__3].i = 0.f; // , expr subst
462 /* L100: */
463 }
464 i__2 = *n - its + 1;
465 i__3 = *n - its + 1;
466 r__1 = *eps3 * rootn;
467 q__1.r = v[i__3].r - r__1;
468 q__1.i = v[i__3].i; // , expr subst
469 v[i__2].r = q__1.r;
470 v[i__2].i = q__1.i; // , expr subst
471 /* L110: */
472 }
473 /* Failure to find eigenvector in N iterations. */
474 *info = 1;
475 L120: /* Normalize eigenvector. */
476 i__ = icamax_(n, &v[1], &c__1);
477 i__1 = i__;
478 r__3 = 1.f / ((r__1 = v[i__1].r, f2c_abs(r__1)) + (r__2 = r_imag(&v[i__]), f2c_abs(r__2)));
479 csscal_(n, &r__3, &v[1], &c__1);
480 return 0;
481 /* End of CLAEIN */
482 }
483 /* claein_ */
484