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