swiftnav: claein.c Source File

Go to the documentation of this file.
00001 /* claein.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Table of constant values */
00017 
00018 static integer c__1 = 1;
00019 
00020 /* Subroutine */ int claein_(logical *rightv, logical *noinit, integer *n, 
00021         complex *h__, integer *ldh, complex *w, complex *v, complex *b, 
00022         integer *ldb, real *rwork, real *eps3, real *smlnum, integer *info)
00023 {
00024     /* System generated locals */
00025     integer b_dim1, b_offset, h_dim1, h_offset, i__1, i__2, i__3, i__4, i__5;
00026     real r__1, r__2, r__3, r__4;
00027     complex q__1, q__2;
00028 
00029     /* Builtin functions */
00030     double sqrt(doublereal), r_imag(complex *);
00031 
00032     /* Local variables */
00033     integer i__, j;
00034     complex x, ei, ej;
00035     integer its, ierr;
00036     complex temp;
00037     real scale;
00038     char trans[1];
00039     real rtemp, rootn, vnorm;
00040     extern doublereal scnrm2_(integer *, complex *, integer *);
00041     extern integer icamax_(integer *, complex *, integer *);
00042     extern /* Complex */ VOID cladiv_(complex *, complex *, complex *);
00043     extern /* Subroutine */ int csscal_(integer *, real *, complex *, integer 
00044             *), clatrs_(char *, char *, char *, char *, integer *, complex *, 
00045             integer *, complex *, real *, real *, integer *);
00046     extern doublereal scasum_(integer *, complex *, integer *);
00047     char normin[1];
00048     real nrmsml, growto;
00049 
00050 
00051 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00052 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00053 /*     November 2006 */
00054 
00055 /*     .. Scalar Arguments .. */
00056 /*     .. */
00057 /*     .. Array Arguments .. */
00058 /*     .. */
00059 
00060 /*  Purpose */
00061 /*  ======= */
00062 
00063 /*  CLAEIN uses inverse iteration to find a right or left eigenvector */
00064 /*  corresponding to the eigenvalue W of a complex upper Hessenberg */
00065 /*  matrix H. */
00066 
00067 /*  Arguments */
00068 /*  ========= */
00069 
00070 /*  RIGHTV   (input) LOGICAL */
00071 /*          = .TRUE. : compute right eigenvector; */
00072 /*          = .FALSE.: compute left eigenvector. */
00073 
00074 /*  NOINIT   (input) LOGICAL */
00075 /*          = .TRUE. : no initial vector supplied in V */
00076 /*          = .FALSE.: initial vector supplied in V. */
00077 
00078 /*  N       (input) INTEGER */
00079 /*          The order of the matrix H.  N >= 0. */
00080 
00081 /*  H       (input) COMPLEX array, dimension (LDH,N) */
00082 /*          The upper Hessenberg matrix H. */
00083 
00084 /*  LDH     (input) INTEGER */
00085 /*          The leading dimension of the array H.  LDH >= max(1,N). */
00086 
00087 /*  W       (input) COMPLEX */
00088 /*          The eigenvalue of H whose corresponding right or left */
00089 /*          eigenvector is to be computed. */
00090 
00091 /*  V       (input/output) COMPLEX array, dimension (N) */
00092 /*          On entry, if NOINIT = .FALSE., V must contain a starting */
00093 /*          vector for inverse iteration; otherwise V need not be set. */
00094 /*          On exit, V contains the computed eigenvector, normalized so */
00095 /*          that the component of largest magnitude has magnitude 1; here */
00096 /*          the magnitude of a complex number (x,y) is taken to be */
00097 /*          |x| + |y|. */
00098 
00099 /*  B       (workspace) COMPLEX array, dimension (LDB,N) */
00100 
00101 /*  LDB     (input) INTEGER */
00102 /*          The leading dimension of the array B.  LDB >= max(1,N). */
00103 
00104 /*  RWORK   (workspace) REAL array, dimension (N) */
00105 
00106 /*  EPS3    (input) REAL */
00107 /*          A small machine-dependent value which is used to perturb */
00108 /*          close eigenvalues, and to replace zero pivots. */
00109 
00110 /*  SMLNUM  (input) REAL */
00111 /*          A machine-dependent value close to the underflow threshold. */
00112 
00113 /*  INFO    (output) INTEGER */
00114 /*          = 0:  successful exit */
00115 /*          = 1:  inverse iteration did not converge; V is set to the */
00116 /*                last iterate. */
00117 
00118 /*  ===================================================================== */
00119 
00120 /*     .. Parameters .. */
00121 /*     .. */
00122 /*     .. Local Scalars .. */
00123 /*     .. */
00124 /*     .. External Functions .. */
00125 /*     .. */
00126 /*     .. External Subroutines .. */
00127 /*     .. */
00128 /*     .. Intrinsic Functions .. */
00129 /*     .. */
00130 /*     .. Statement Functions .. */
00131 /*     .. */
00132 /*     .. Statement Function definitions .. */
00133 /*     .. */
00134 /*     .. Executable Statements .. */
00135 
00136     /* Parameter adjustments */
00137     h_dim1 = *ldh;
00138     h_offset = 1 + h_dim1;
00139     h__ -= h_offset;
00140     --v;
00141     b_dim1 = *ldb;
00142     b_offset = 1 + b_dim1;
00143     b -= b_offset;
00144     --rwork;
00145 
00146     /* Function Body */
00147     *info = 0;
00148 
00149 /*     GROWTO is the threshold used in the acceptance test for an */
00150 /*     eigenvector. */
00151 
00152     rootn = sqrt((real) (*n));
00153     growto = .1f / rootn;
00154 /* Computing MAX */
00155     r__1 = 1.f, r__2 = *eps3 * rootn;
00156     nrmsml = dmax(r__1,r__2) * *smlnum;
00157 
00158 /*     Form B = H - W*I (except that the subdiagonal elements are not */
00159 /*     stored). */
00160 
00161     i__1 = *n;
00162     for (j = 1; j <= i__1; ++j) {
00163         i__2 = j - 1;
00164         for (i__ = 1; i__ <= i__2; ++i__) {
00165             i__3 = i__ + j * b_dim1;
00166             i__4 = i__ + j * h_dim1;
00167             b[i__3].r = h__[i__4].r, b[i__3].i = h__[i__4].i;
00168 /* L10: */
00169         }
00170         i__2 = j + j * b_dim1;
00171         i__3 = j + j * h_dim1;
00172         q__1.r = h__[i__3].r - w->r, q__1.i = h__[i__3].i - w->i;
00173         b[i__2].r = q__1.r, b[i__2].i = q__1.i;
00174 /* L20: */
00175     }
00176 
00177     if (*noinit) {
00178 
00179 /*        Initialize V. */
00180 
00181         i__1 = *n;
00182         for (i__ = 1; i__ <= i__1; ++i__) {
00183             i__2 = i__;
00184             v[i__2].r = *eps3, v[i__2].i = 0.f;
00185 /* L30: */
00186         }
00187     } else {
00188 
00189 /*        Scale supplied initial vector. */
00190 
00191         vnorm = scnrm2_(n, &v[1], &c__1);
00192         r__1 = *eps3 * rootn / dmax(vnorm,nrmsml);
00193         csscal_(n, &r__1, &v[1], &c__1);
00194     }
00195 
00196     if (*rightv) {
00197 
00198 /*        LU decomposition with partial pivoting of B, replacing zero */
00199 /*        pivots by EPS3. */
00200 
00201         i__1 = *n - 1;
00202         for (i__ = 1; i__ <= i__1; ++i__) {
00203             i__2 = i__ + 1 + i__ * h_dim1;
00204             ei.r = h__[i__2].r, ei.i = h__[i__2].i;
00205             i__2 = i__ + i__ * b_dim1;
00206             if ((r__1 = b[i__2].r, dabs(r__1)) + (r__2 = r_imag(&b[i__ + i__ *
00207                      b_dim1]), dabs(r__2)) < (r__3 = ei.r, dabs(r__3)) + (
00208                     r__4 = r_imag(&ei), dabs(r__4))) {
00209 
00210 /*              Interchange rows and eliminate. */
00211 
00212                 cladiv_(&q__1, &b[i__ + i__ * b_dim1], &ei);
00213                 x.r = q__1.r, x.i = q__1.i;
00214                 i__2 = i__ + i__ * b_dim1;
00215                 b[i__2].r = ei.r, b[i__2].i = ei.i;
00216                 i__2 = *n;
00217                 for (j = i__ + 1; j <= i__2; ++j) {
00218                     i__3 = i__ + 1 + j * b_dim1;
00219                     temp.r = b[i__3].r, temp.i = b[i__3].i;
00220                     i__3 = i__ + 1 + j * b_dim1;
00221                     i__4 = i__ + j * b_dim1;
00222                     q__2.r = x.r * temp.r - x.i * temp.i, q__2.i = x.r * 
00223                             temp.i + x.i * temp.r;
00224                     q__1.r = b[i__4].r - q__2.r, q__1.i = b[i__4].i - q__2.i;
00225                     b[i__3].r = q__1.r, b[i__3].i = q__1.i;
00226                     i__3 = i__ + j * b_dim1;
00227                     b[i__3].r = temp.r, b[i__3].i = temp.i;
00228 /* L40: */
00229                 }
00230             } else {
00231 
00232 /*              Eliminate without interchange. */
00233 
00234                 i__2 = i__ + i__ * b_dim1;
00235                 if (b[i__2].r == 0.f && b[i__2].i == 0.f) {
00236                     i__3 = i__ + i__ * b_dim1;
00237                     b[i__3].r = *eps3, b[i__3].i = 0.f;
00238                 }
00239                 cladiv_(&q__1, &ei, &b[i__ + i__ * b_dim1]);
00240                 x.r = q__1.r, x.i = q__1.i;
00241                 if (x.r != 0.f || x.i != 0.f) {
00242                     i__2 = *n;
00243                     for (j = i__ + 1; j <= i__2; ++j) {
00244                         i__3 = i__ + 1 + j * b_dim1;
00245                         i__4 = i__ + 1 + j * b_dim1;
00246                         i__5 = i__ + j * b_dim1;
00247                         q__2.r = x.r * b[i__5].r - x.i * b[i__5].i, q__2.i = 
00248                                 x.r * b[i__5].i + x.i * b[i__5].r;
00249                         q__1.r = b[i__4].r - q__2.r, q__1.i = b[i__4].i - 
00250                                 q__2.i;
00251                         b[i__3].r = q__1.r, b[i__3].i = q__1.i;
00252 /* L50: */
00253                     }
00254                 }
00255             }
00256 /* L60: */
00257         }
00258         i__1 = *n + *n * b_dim1;
00259         if (b[i__1].r == 0.f && b[i__1].i == 0.f) {
00260             i__2 = *n + *n * b_dim1;
00261             b[i__2].r = *eps3, b[i__2].i = 0.f;
00262         }
00263 
00264         *(unsigned char *)trans = 'N';
00265 
00266     } else {
00267 
00268 /*        UL decomposition with partial pivoting of B, replacing zero */
00269 /*        pivots by EPS3. */
00270 
00271         for (j = *n; j >= 2; --j) {
00272             i__1 = j + (j - 1) * h_dim1;
00273             ej.r = h__[i__1].r, ej.i = h__[i__1].i;
00274             i__1 = j + j * b_dim1;
00275             if ((r__1 = b[i__1].r, dabs(r__1)) + (r__2 = r_imag(&b[j + j * 
00276                     b_dim1]), dabs(r__2)) < (r__3 = ej.r, dabs(r__3)) + (r__4 
00277                     = r_imag(&ej), dabs(r__4))) {
00278 
00279 /*              Interchange columns and eliminate. */
00280 
00281                 cladiv_(&q__1, &b[j + j * b_dim1], &ej);
00282                 x.r = q__1.r, x.i = q__1.i;
00283                 i__1 = j + j * b_dim1;
00284                 b[i__1].r = ej.r, b[i__1].i = ej.i;
00285                 i__1 = j - 1;
00286                 for (i__ = 1; i__ <= i__1; ++i__) {
00287                     i__2 = i__ + (j - 1) * b_dim1;
00288                     temp.r = b[i__2].r, temp.i = b[i__2].i;
00289                     i__2 = i__ + (j - 1) * b_dim1;
00290                     i__3 = i__ + j * b_dim1;
00291                     q__2.r = x.r * temp.r - x.i * temp.i, q__2.i = x.r * 
00292                             temp.i + x.i * temp.r;
00293                     q__1.r = b[i__3].r - q__2.r, q__1.i = b[i__3].i - q__2.i;
00294                     b[i__2].r = q__1.r, b[i__2].i = q__1.i;
00295                     i__2 = i__ + j * b_dim1;
00296                     b[i__2].r = temp.r, b[i__2].i = temp.i;
00297 /* L70: */
00298                 }
00299             } else {
00300 
00301 /*              Eliminate without interchange. */
00302 
00303                 i__1 = j + j * b_dim1;
00304                 if (b[i__1].r == 0.f && b[i__1].i == 0.f) {
00305                     i__2 = j + j * b_dim1;
00306                     b[i__2].r = *eps3, b[i__2].i = 0.f;
00307                 }
00308                 cladiv_(&q__1, &ej, &b[j + j * b_dim1]);
00309                 x.r = q__1.r, x.i = q__1.i;
00310                 if (x.r != 0.f || x.i != 0.f) {
00311                     i__1 = j - 1;
00312                     for (i__ = 1; i__ <= i__1; ++i__) {
00313                         i__2 = i__ + (j - 1) * b_dim1;
00314                         i__3 = i__ + (j - 1) * b_dim1;
00315                         i__4 = i__ + j * b_dim1;
00316                         q__2.r = x.r * b[i__4].r - x.i * b[i__4].i, q__2.i = 
00317                                 x.r * b[i__4].i + x.i * b[i__4].r;
00318                         q__1.r = b[i__3].r - q__2.r, q__1.i = b[i__3].i - 
00319                                 q__2.i;
00320                         b[i__2].r = q__1.r, b[i__2].i = q__1.i;
00321 /* L80: */
00322                     }
00323                 }
00324             }
00325 /* L90: */
00326         }
00327         i__1 = b_dim1 + 1;
00328         if (b[i__1].r == 0.f && b[i__1].i == 0.f) {
00329             i__2 = b_dim1 + 1;
00330             b[i__2].r = *eps3, b[i__2].i = 0.f;
00331         }
00332 
00333         *(unsigned char *)trans = 'C';
00334 
00335     }
00336 
00337     *(unsigned char *)normin = 'N';
00338     i__1 = *n;
00339     for (its = 1; its <= i__1; ++its) {
00340 
00341 /*        Solve U*x = scale*v for a right eigenvector */
00342 /*          or U'*x = scale*v for a left eigenvector, */
00343 /*        overwriting x on v. */
00344 
00345         clatrs_("Upper", trans, "Nonunit", normin, n, &b[b_offset], ldb, &v[1]
00346 , &scale, &rwork[1], &ierr);
00347         *(unsigned char *)normin = 'Y';
00348 
00349 /*        Test for sufficient growth in the norm of v. */
00350 
00351         vnorm = scasum_(n, &v[1], &c__1);
00352         if (vnorm >= growto * scale) {
00353             goto L120;
00354         }
00355 
00356 /*        Choose new orthogonal starting vector and try again. */
00357 
00358         rtemp = *eps3 / (rootn + 1.f);
00359         v[1].r = *eps3, v[1].i = 0.f;
00360         i__2 = *n;
00361         for (i__ = 2; i__ <= i__2; ++i__) {
00362             i__3 = i__;
00363             v[i__3].r = rtemp, v[i__3].i = 0.f;
00364 /* L100: */
00365         }
00366         i__2 = *n - its + 1;
00367         i__3 = *n - its + 1;
00368         r__1 = *eps3 * rootn;
00369         q__1.r = v[i__3].r - r__1, q__1.i = v[i__3].i;
00370         v[i__2].r = q__1.r, v[i__2].i = q__1.i;
00371 /* L110: */
00372     }
00373 
00374 /*     Failure to find eigenvector in N iterations. */
00375 
00376     *info = 1;
00377 
00378 L120:
00379 
00380 /*     Normalize eigenvector. */
00381 
00382     i__ = icamax_(n, &v[1], &c__1);
00383     i__1 = i__;
00384     r__3 = 1.f / ((r__1 = v[i__1].r, dabs(r__1)) + (r__2 = r_imag(&v[i__]), 
00385             dabs(r__2)));
00386     csscal_(n, &r__3, &v[1], &c__1);
00387 
00388     return 0;
00389 
00390 /*     End of CLAEIN */
00391 
00392 } /* claein_ */