1 /* ../netlib/zlaein.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 ZLAEIN 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 ZLAEIN + dependencies */
11 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaein. f"> */
12 /* > [TGZ]</a> */
13 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaein. f"> */
14 /* > [ZIP]</a> */
15 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaein. f"> */
16 /* > [TXT]</a> */
17 /* > \endhtmlonly */
18 /* Definition: */
19 /* =========== */
20 /* SUBROUTINE ZLAEIN( 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 /* DOUBLE PRECISION EPS3, SMLNUM */
26 /* COMPLEX*16 W */
27 /* .. */
28 /* .. Array Arguments .. */
29 /* DOUBLE PRECISION RWORK( * ) */
30 /* COMPLEX*16 B( LDB, * ), H( LDH, * ), V( * ) */
31 /* .. */
32 /* > \par Purpose: */
33 /* ============= */
34 /* > */
35 /* > \verbatim */
36 /* > */
37 /* > ZLAEIN 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*16 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*16 */
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*16 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*16 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 DOUBLE PRECISION array, dimension (N) */
110 /* > \endverbatim */
111 /* > */
112 /* > \param[in] EPS3 */
113 /* > \verbatim */
114 /* > EPS3 is DOUBLE PRECISION */
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 DOUBLE PRECISION */
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 complex16OTHERauxiliary */
141 /* ===================================================================== */
142 /* Subroutine */
zlaein_(logical * rightv,logical * noinit,integer * n,doublecomplex * h__,integer * ldh,doublecomplex * w,doublecomplex * v,doublecomplex * b,integer * ldb,doublereal * rwork,doublereal * eps3,doublereal * smlnum,integer * info)143 int zlaein_(logical *rightv, logical *noinit, integer *n, doublecomplex *h__, integer *ldh, doublecomplex *w, doublecomplex *v, doublecomplex *b, integer *ldb, doublereal *rwork, doublereal *eps3, doublereal *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 doublereal d__1, d__2, d__3, d__4;
148 doublecomplex z__1, z__2;
149 /* Builtin functions */
150 double sqrt(doublereal), d_imag(doublecomplex *);
151 /* Local variables */
152 integer i__, j;
153 doublecomplex x, ei, ej;
154 integer its, ierr;
155 doublecomplex temp;
156 doublereal scale;
157 char trans[1];
158 doublereal rtemp, rootn, vnorm;
159 extern doublereal dznrm2_(integer *, doublecomplex *, integer *);
160 extern /* Subroutine */
161 int zdscal_(integer *, doublereal *, doublecomplex *, integer *);
162 extern integer izamax_(integer *, doublecomplex *, integer *);
163 extern /* Double Complex */
164 VOID zladiv_(doublecomplex *, doublecomplex *, doublecomplex *);
165 char normin[1];
166 extern doublereal dzasum_(integer *, doublecomplex *, integer *);
167 doublereal nrmsml;
168 extern /* Subroutine */
169 int zlatrs_(char *, char *, char *, char *, integer *, doublecomplex *, integer *, doublecomplex *, doublereal *, doublereal *, integer *);
170 doublereal growto;
171 /* -- LAPACK auxiliary routine (version 3.4.2) -- */
172 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
173 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
174 /* September 2012 */
175 /* .. Scalar Arguments .. */
176 /* .. */
177 /* .. Array Arguments .. */
178 /* .. */
179 /* ===================================================================== */
180 /* .. Parameters .. */
181 /* .. */
182 /* .. Local Scalars .. */
183 /* .. */
184 /* .. External Functions .. */
185 /* .. */
186 /* .. External Subroutines .. */
187 /* .. */
188 /* .. Intrinsic Functions .. */
189 /* .. */
190 /* .. Statement Functions .. */
191 /* .. */
192 /* .. Statement Function definitions .. */
193 /* .. */
194 /* .. Executable Statements .. */
195 /* Parameter adjustments */
196 h_dim1 = *ldh;
197 h_offset = 1 + h_dim1;
198 h__ -= h_offset;
199 --v;
200 b_dim1 = *ldb;
201 b_offset = 1 + b_dim1;
202 b -= b_offset;
203 --rwork;
204 /* Function Body */
205 *info = 0;
206 /* GROWTO is the threshold used in the acceptance test for an */
207 /* eigenvector. */
208 rootn = sqrt((doublereal) (*n));
209 growto = .1 / rootn;
210 /* Computing MAX */
211 d__1 = 1.;
212 d__2 = *eps3 * rootn; // , expr subst
213 nrmsml = max(d__1,d__2) * *smlnum;
214 /* Form B = H - W*I (except that the subdiagonal elements are not */
215 /* stored). */
216 i__1 = *n;
217 for (j = 1;
218 j <= i__1;
219 ++j)
220 {
221 i__2 = j - 1;
222 for (i__ = 1;
223 i__ <= i__2;
224 ++i__)
225 {
226 i__3 = i__ + j * b_dim1;
227 i__4 = i__ + j * h_dim1;
228 b[i__3].r = h__[i__4].r;
229 b[i__3].i = h__[i__4].i; // , expr subst
230 /* L10: */
231 }
232 i__2 = j + j * b_dim1;
233 i__3 = j + j * h_dim1;
234 z__1.r = h__[i__3].r - w->r;
235 z__1.i = h__[i__3].i - w->i; // , expr subst
236 b[i__2].r = z__1.r;
237 b[i__2].i = z__1.i; // , expr subst
238 /* L20: */
239 }
240 if (*noinit)
241 {
242 /* Initialize V. */
243 i__1 = *n;
244 for (i__ = 1;
245 i__ <= i__1;
246 ++i__)
247 {
248 i__2 = i__;
249 v[i__2].r = *eps3;
250 v[i__2].i = 0.; // , expr subst
251 /* L30: */
252 }
253 }
254 else
255 {
256 /* Scale supplied initial vector. */
257 vnorm = dznrm2_(n, &v[1], &c__1);
258 d__1 = *eps3 * rootn / max(vnorm,nrmsml);
259 zdscal_(n, &d__1, &v[1], &c__1);
260 }
261 if (*rightv)
262 {
263 /* LU decomposition with partial pivoting of B, replacing zero */
264 /* pivots by EPS3. */
265 i__1 = *n - 1;
266 for (i__ = 1;
267 i__ <= i__1;
268 ++i__)
269 {
270 i__2 = i__ + 1 + i__ * h_dim1;
271 ei.r = h__[i__2].r;
272 ei.i = h__[i__2].i; // , expr subst
273 i__2 = i__ + i__ * b_dim1;
274 if ((d__1 = b[i__2].r, f2c_abs(d__1)) + (d__2 = d_imag(&b[i__ + i__ * b_dim1]), f2c_abs(d__2)) < (d__3 = ei.r, f2c_abs(d__3)) + (d__4 = d_imag(&ei), f2c_abs(d__4)))
275 {
276 /* Interchange rows and eliminate. */
277 zladiv_(&z__1, &b[i__ + i__ * b_dim1], &ei);
278 x.r = z__1.r;
279 x.i = z__1.i; // , expr subst
280 i__2 = i__ + i__ * b_dim1;
281 b[i__2].r = ei.r;
282 b[i__2].i = ei.i; // , expr subst
283 i__2 = *n;
284 for (j = i__ + 1;
285 j <= i__2;
286 ++j)
287 {
288 i__3 = i__ + 1 + j * b_dim1;
289 temp.r = b[i__3].r;
290 temp.i = b[i__3].i; // , expr subst
291 i__3 = i__ + 1 + j * b_dim1;
292 i__4 = i__ + j * b_dim1;
293 z__2.r = x.r * temp.r - x.i * temp.i;
294 z__2.i = x.r * temp.i + x.i * temp.r; // , expr subst
295 z__1.r = b[i__4].r - z__2.r;
296 z__1.i = b[i__4].i - z__2.i; // , expr subst
297 b[i__3].r = z__1.r;
298 b[i__3].i = z__1.i; // , expr subst
299 i__3 = i__ + j * b_dim1;
300 b[i__3].r = temp.r;
301 b[i__3].i = temp.i; // , expr subst
302 /* L40: */
303 }
304 }
305 else
306 {
307 /* Eliminate without interchange. */
308 i__2 = i__ + i__ * b_dim1;
309 if (b[i__2].r == 0. && b[i__2].i == 0.)
310 {
311 i__3 = i__ + i__ * b_dim1;
312 b[i__3].r = *eps3;
313 b[i__3].i = 0.; // , expr subst
314 }
315 zladiv_(&z__1, &ei, &b[i__ + i__ * b_dim1]);
316 x.r = z__1.r;
317 x.i = z__1.i; // , expr subst
318 if (x.r != 0. || x.i != 0.)
319 {
320 i__2 = *n;
321 for (j = i__ + 1;
322 j <= i__2;
323 ++j)
324 {
325 i__3 = i__ + 1 + j * b_dim1;
326 i__4 = i__ + 1 + j * b_dim1;
327 i__5 = i__ + j * b_dim1;
328 z__2.r = x.r * b[i__5].r - x.i * b[i__5].i;
329 z__2.i = x.r * b[i__5].i + x.i * b[i__5].r; // , expr subst
330 z__1.r = b[i__4].r - z__2.r;
331 z__1.i = b[i__4].i - z__2.i; // , expr subst
332 b[i__3].r = z__1.r;
333 b[i__3].i = z__1.i; // , expr subst
334 /* L50: */
335 }
336 }
337 }
338 /* L60: */
339 }
340 i__1 = *n + *n * b_dim1;
341 if (b[i__1].r == 0. && b[i__1].i == 0.)
342 {
343 i__2 = *n + *n * b_dim1;
344 b[i__2].r = *eps3;
345 b[i__2].i = 0.; // , expr subst
346 }
347 *(unsigned char *)trans = 'N';
348 }
349 else
350 {
351 /* UL decomposition with partial pivoting of B, replacing zero */
352 /* pivots by EPS3. */
353 for (j = *n;
354 j >= 2;
355 --j)
356 {
357 i__1 = j + (j - 1) * h_dim1;
358 ej.r = h__[i__1].r;
359 ej.i = h__[i__1].i; // , expr subst
360 i__1 = j + j * b_dim1;
361 if ((d__1 = b[i__1].r, f2c_abs(d__1)) + (d__2 = d_imag(&b[j + j * b_dim1]), f2c_abs(d__2)) < (d__3 = ej.r, f2c_abs(d__3)) + (d__4 = d_imag(&ej), f2c_abs(d__4)))
362 {
363 /* Interchange columns and eliminate. */
364 zladiv_(&z__1, &b[j + j * b_dim1], &ej);
365 x.r = z__1.r;
366 x.i = z__1.i; // , expr subst
367 i__1 = j + j * b_dim1;
368 b[i__1].r = ej.r;
369 b[i__1].i = ej.i; // , expr subst
370 i__1 = j - 1;
371 for (i__ = 1;
372 i__ <= i__1;
373 ++i__)
374 {
375 i__2 = i__ + (j - 1) * b_dim1;
376 temp.r = b[i__2].r;
377 temp.i = b[i__2].i; // , expr subst
378 i__2 = i__ + (j - 1) * b_dim1;
379 i__3 = i__ + j * b_dim1;
380 z__2.r = x.r * temp.r - x.i * temp.i;
381 z__2.i = x.r * temp.i + x.i * temp.r; // , expr subst
382 z__1.r = b[i__3].r - z__2.r;
383 z__1.i = b[i__3].i - z__2.i; // , expr subst
384 b[i__2].r = z__1.r;
385 b[i__2].i = z__1.i; // , expr subst
386 i__2 = i__ + j * b_dim1;
387 b[i__2].r = temp.r;
388 b[i__2].i = temp.i; // , expr subst
389 /* L70: */
390 }
391 }
392 else
393 {
394 /* Eliminate without interchange. */
395 i__1 = j + j * b_dim1;
396 if (b[i__1].r == 0. && b[i__1].i == 0.)
397 {
398 i__2 = j + j * b_dim1;
399 b[i__2].r = *eps3;
400 b[i__2].i = 0.; // , expr subst
401 }
402 zladiv_(&z__1, &ej, &b[j + j * b_dim1]);
403 x.r = z__1.r;
404 x.i = z__1.i; // , expr subst
405 if (x.r != 0. || x.i != 0.)
406 {
407 i__1 = j - 1;
408 for (i__ = 1;
409 i__ <= i__1;
410 ++i__)
411 {
412 i__2 = i__ + (j - 1) * b_dim1;
413 i__3 = i__ + (j - 1) * b_dim1;
414 i__4 = i__ + j * b_dim1;
415 z__2.r = x.r * b[i__4].r - x.i * b[i__4].i;
416 z__2.i = x.r * b[i__4].i + x.i * b[i__4].r; // , expr subst
417 z__1.r = b[i__3].r - z__2.r;
418 z__1.i = b[i__3].i - z__2.i; // , expr subst
419 b[i__2].r = z__1.r;
420 b[i__2].i = z__1.i; // , expr subst
421 /* L80: */
422 }
423 }
424 }
425 /* L90: */
426 }
427 i__1 = b_dim1 + 1;
428 if (b[i__1].r == 0. && b[i__1].i == 0.)
429 {
430 i__2 = b_dim1 + 1;
431 b[i__2].r = *eps3;
432 b[i__2].i = 0.; // , expr subst
433 }
434 *(unsigned char *)trans = 'C';
435 }
436 *(unsigned char *)normin = 'N';
437 i__1 = *n;
438 for (its = 1;
439 its <= i__1;
440 ++its)
441 {
442 /* Solve U*x = scale*v for a right eigenvector */
443 /* or U**H *x = scale*v for a left eigenvector, */
444 /* overwriting x on v. */
445 zlatrs_("Upper", trans, "Nonunit", normin, n, &b[b_offset], ldb, &v[1] , &scale, &rwork[1], &ierr);
446 *(unsigned char *)normin = 'Y';
447 /* Test for sufficient growth in the norm of v. */
448 vnorm = dzasum_(n, &v[1], &c__1);
449 if (vnorm >= growto * scale)
450 {
451 goto L120;
452 }
453 /* Choose new orthogonal starting vector and try again. */
454 rtemp = *eps3 / (rootn + 1.);
455 v[1].r = *eps3;
456 v[1].i = 0.; // , expr subst
457 i__2 = *n;
458 for (i__ = 2;
459 i__ <= i__2;
460 ++i__)
461 {
462 i__3 = i__;
463 v[i__3].r = rtemp;
464 v[i__3].i = 0.; // , expr subst
465 /* L100: */
466 }
467 i__2 = *n - its + 1;
468 i__3 = *n - its + 1;
469 d__1 = *eps3 * rootn;
470 z__1.r = v[i__3].r - d__1;
471 z__1.i = v[i__3].i; // , expr subst
472 v[i__2].r = z__1.r;
473 v[i__2].i = z__1.i; // , expr subst
474 /* L110: */
475 }
476 /* Failure to find eigenvector in N iterations. */
477 *info = 1;
478 L120: /* Normalize eigenvector. */
479 i__ = izamax_(n, &v[1], &c__1);
480 i__1 = i__;
481 d__3 = 1. / ((d__1 = v[i__1].r, f2c_abs(d__1)) + (d__2 = d_imag(&v[i__]), f2c_abs( d__2)));
482 zdscal_(n, &d__3, &v[1], &c__1);
483 return 0;
484 /* End of ZLAEIN */
485 }
486 /* zlaein_ */
487