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