#include "blaswrap.h" #include "f2c.h" /* Subroutine */ 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) { /* -- LAPACK auxiliary routine (version 3.1) -- Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. November 2006 Purpose ======= CLAEIN uses inverse iteration to find a right or left eigenvector corresponding to the eigenvalue W of a complex upper Hessenberg matrix H. Arguments ========= RIGHTV (input) LOGICAL = .TRUE. : compute right eigenvector; = .FALSE.: compute left eigenvector. NOINIT (input) LOGICAL = .TRUE. : no initial vector supplied in V = .FALSE.: initial vector supplied in V. N (input) INTEGER The order of the matrix H. N >= 0. H (input) COMPLEX array, dimension (LDH,N) The upper Hessenberg matrix H. LDH (input) INTEGER The leading dimension of the array H. LDH >= max(1,N). W (input) COMPLEX The eigenvalue of H whose corresponding right or left eigenvector is to be computed. V (input/output) COMPLEX array, dimension (N) On entry, if NOINIT = .FALSE., V must contain a starting vector for inverse iteration; otherwise V need not be set. On exit, V contains the computed eigenvector, normalized so that the component of largest magnitude has magnitude 1; here the magnitude of a complex number (x,y) is taken to be |x| + |y|. B (workspace) COMPLEX array, dimension (LDB,N) LDB (input) INTEGER The leading dimension of the array B. LDB >= max(1,N). RWORK (workspace) REAL array, dimension (N) EPS3 (input) REAL A small machine-dependent value which is used to perturb close eigenvalues, and to replace zero pivots. SMLNUM (input) REAL A machine-dependent value close to the underflow threshold. INFO (output) INTEGER = 0: successful exit = 1: inverse iteration did not converge; V is set to the last iterate. ===================================================================== Parameter adjustments */ /* Table of constant values */ static integer c__1 = 1; /* System generated locals */ integer b_dim1, b_offset, h_dim1, h_offset, i__1, i__2, i__3, i__4, i__5; real r__1, r__2, r__3, r__4; complex q__1, q__2; /* Builtin functions */ double sqrt(doublereal), r_imag(complex *); /* Local variables */ static integer i__, j; static complex x, ei, ej; static integer its, ierr; static complex temp; static real scale; static char trans[1]; static real rtemp, rootn, vnorm; extern doublereal scnrm2_(integer *, complex *, integer *); extern integer icamax_(integer *, complex *, integer *); extern /* Complex */ VOID cladiv_(complex *, complex *, complex *); extern /* Subroutine */ int csscal_(integer *, real *, complex *, integer *), clatrs_(char *, char *, char *, char *, integer *, complex *, integer *, complex *, real *, real *, integer *); extern doublereal scasum_(integer *, complex *, integer *); static char normin[1]; static real nrmsml, growto; h_dim1 = *ldh; h_offset = 1 + h_dim1; h__ -= h_offset; --v; b_dim1 = *ldb; b_offset = 1 + b_dim1; b -= b_offset; --rwork; /* Function Body */ *info = 0; /* GROWTO is the threshold used in the acceptance test for an eigenvector. */ rootn = sqrt((real) (*n)); growto = .1f / rootn; /* Computing MAX */ r__1 = 1.f, r__2 = *eps3 * rootn; nrmsml = dmax(r__1,r__2) * *smlnum; /* Form B = H - W*I (except that the subdiagonal elements are not stored). */ i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = j - 1; for (i__ = 1; i__ <= i__2; ++i__) { i__3 = i__ + j * b_dim1; i__4 = i__ + j * h_dim1; b[i__3].r = h__[i__4].r, b[i__3].i = h__[i__4].i; /* L10: */ } i__2 = j + j * b_dim1; i__3 = j + j * h_dim1; q__1.r = h__[i__3].r - w->r, q__1.i = h__[i__3].i - w->i; b[i__2].r = q__1.r, b[i__2].i = q__1.i; /* L20: */ } if (*noinit) { /* Initialize V. */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = i__; v[i__2].r = *eps3, v[i__2].i = 0.f; /* L30: */ } } else { /* Scale supplied initial vector. */ vnorm = scnrm2_(n, &v[1], &c__1); r__1 = *eps3 * rootn / dmax(vnorm,nrmsml); csscal_(n, &r__1, &v[1], &c__1); } if (*rightv) { /* LU decomposition with partial pivoting of B, replacing zero pivots by EPS3. */ i__1 = *n - 1; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = i__ + 1 + i__ * h_dim1; ei.r = h__[i__2].r, ei.i = h__[i__2].i; i__2 = i__ + i__ * b_dim1; if ((r__1 = b[i__2].r, dabs(r__1)) + (r__2 = r_imag(&b[i__ + i__ * b_dim1]), dabs(r__2)) < (r__3 = ei.r, dabs(r__3)) + ( r__4 = r_imag(&ei), dabs(r__4))) { /* Interchange rows and eliminate. */ cladiv_(&q__1, &b[i__ + i__ * b_dim1], &ei); x.r = q__1.r, x.i = q__1.i; i__2 = i__ + i__ * b_dim1; b[i__2].r = ei.r, b[i__2].i = ei.i; i__2 = *n; for (j = i__ + 1; j <= i__2; ++j) { i__3 = i__ + 1 + j * b_dim1; temp.r = b[i__3].r, temp.i = b[i__3].i; i__3 = i__ + 1 + j * b_dim1; i__4 = i__ + j * b_dim1; q__2.r = x.r * temp.r - x.i * temp.i, q__2.i = x.r * temp.i + x.i * temp.r; q__1.r = b[i__4].r - q__2.r, q__1.i = b[i__4].i - q__2.i; b[i__3].r = q__1.r, b[i__3].i = q__1.i; i__3 = i__ + j * b_dim1; b[i__3].r = temp.r, b[i__3].i = temp.i; /* L40: */ } } else { /* Eliminate without interchange. */ i__2 = i__ + i__ * b_dim1; if (b[i__2].r == 0.f && b[i__2].i == 0.f) { i__3 = i__ + i__ * b_dim1; b[i__3].r = *eps3, b[i__3].i = 0.f; } cladiv_(&q__1, &ei, &b[i__ + i__ * b_dim1]); x.r = q__1.r, x.i = q__1.i; if (x.r != 0.f || x.i != 0.f) { i__2 = *n; for (j = i__ + 1; j <= i__2; ++j) { i__3 = i__ + 1 + j * b_dim1; i__4 = i__ + 1 + j * b_dim1; i__5 = i__ + j * b_dim1; q__2.r = x.r * b[i__5].r - x.i * b[i__5].i, q__2.i = x.r * b[i__5].i + x.i * b[i__5].r; q__1.r = b[i__4].r - q__2.r, q__1.i = b[i__4].i - q__2.i; b[i__3].r = q__1.r, b[i__3].i = q__1.i; /* L50: */ } } } /* L60: */ } i__1 = *n + *n * b_dim1; if (b[i__1].r == 0.f && b[i__1].i == 0.f) { i__2 = *n + *n * b_dim1; b[i__2].r = *eps3, b[i__2].i = 0.f; } *(unsigned char *)trans = 'N'; } else { /* UL decomposition with partial pivoting of B, replacing zero pivots by EPS3. */ for (j = *n; j >= 2; --j) { i__1 = j + (j - 1) * h_dim1; ej.r = h__[i__1].r, ej.i = h__[i__1].i; i__1 = j + j * b_dim1; if ((r__1 = b[i__1].r, dabs(r__1)) + (r__2 = r_imag(&b[j + j * b_dim1]), dabs(r__2)) < (r__3 = ej.r, dabs(r__3)) + (r__4 = r_imag(&ej), dabs(r__4))) { /* Interchange columns and eliminate. */ cladiv_(&q__1, &b[j + j * b_dim1], &ej); x.r = q__1.r, x.i = q__1.i; i__1 = j + j * b_dim1; b[i__1].r = ej.r, b[i__1].i = ej.i; i__1 = j - 1; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = i__ + (j - 1) * b_dim1; temp.r = b[i__2].r, temp.i = b[i__2].i; i__2 = i__ + (j - 1) * b_dim1; i__3 = i__ + j * b_dim1; q__2.r = x.r * temp.r - x.i * temp.i, q__2.i = x.r * temp.i + x.i * temp.r; q__1.r = b[i__3].r - q__2.r, q__1.i = b[i__3].i - q__2.i; b[i__2].r = q__1.r, b[i__2].i = q__1.i; i__2 = i__ + j * b_dim1; b[i__2].r = temp.r, b[i__2].i = temp.i; /* L70: */ } } else { /* Eliminate without interchange. */ i__1 = j + j * b_dim1; if (b[i__1].r == 0.f && b[i__1].i == 0.f) { i__2 = j + j * b_dim1; b[i__2].r = *eps3, b[i__2].i = 0.f; } cladiv_(&q__1, &ej, &b[j + j * b_dim1]); x.r = q__1.r, x.i = q__1.i; if (x.r != 0.f || x.i != 0.f) { i__1 = j - 1; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = i__ + (j - 1) * b_dim1; i__3 = i__ + (j - 1) * b_dim1; i__4 = i__ + j * b_dim1; q__2.r = x.r * b[i__4].r - x.i * b[i__4].i, q__2.i = x.r * b[i__4].i + x.i * b[i__4].r; q__1.r = b[i__3].r - q__2.r, q__1.i = b[i__3].i - q__2.i; b[i__2].r = q__1.r, b[i__2].i = q__1.i; /* L80: */ } } } /* L90: */ } i__1 = b_dim1 + 1; if (b[i__1].r == 0.f && b[i__1].i == 0.f) { i__2 = b_dim1 + 1; b[i__2].r = *eps3, b[i__2].i = 0.f; } *(unsigned char *)trans = 'C'; } *(unsigned char *)normin = 'N'; i__1 = *n; for (its = 1; its <= i__1; ++its) { /* Solve U*x = scale*v for a right eigenvector or U'*x = scale*v for a left eigenvector, overwriting x on v. */ clatrs_("Upper", trans, "Nonunit", normin, n, &b[b_offset], ldb, &v[1] , &scale, &rwork[1], &ierr); *(unsigned char *)normin = 'Y'; /* Test for sufficient growth in the norm of v. */ vnorm = scasum_(n, &v[1], &c__1); if (vnorm >= growto * scale) { goto L120; } /* Choose new orthogonal starting vector and try again. */ rtemp = *eps3 / (rootn + 1.f); v[1].r = *eps3, v[1].i = 0.f; i__2 = *n; for (i__ = 2; i__ <= i__2; ++i__) { i__3 = i__; v[i__3].r = rtemp, v[i__3].i = 0.f; /* L100: */ } i__2 = *n - its + 1; i__3 = *n - its + 1; r__1 = *eps3 * rootn; q__1.r = v[i__3].r - r__1, q__1.i = v[i__3].i; v[i__2].r = q__1.r, v[i__2].i = q__1.i; /* L110: */ } /* Failure to find eigenvector in N iterations. */ *info = 1; L120: /* Normalize eigenvector. */ i__ = icamax_(n, &v[1], &c__1); i__1 = i__; r__3 = 1.f / ((r__1 = v[i__1].r, dabs(r__1)) + (r__2 = r_imag(&v[i__]), dabs(r__2))); csscal_(n, &r__3, &v[1], &c__1); return 0; /* End of CLAEIN */ } /* claein_ */