cgebal.c
Go to the documentation of this file.
00001 /* cgebal.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 cgebal_(char *job, integer *n, complex *a, integer *lda, 
00021         integer *ilo, integer *ihi, real *scale, integer *info)
00022 {
00023     /* System generated locals */
00024     integer a_dim1, a_offset, i__1, i__2, i__3;
00025     real r__1, r__2;
00026 
00027     /* Builtin functions */
00028     double r_imag(complex *), c_abs(complex *);
00029 
00030     /* Local variables */
00031     real c__, f, g;
00032     integer i__, j, k, l, m;
00033     real r__, s, ca, ra;
00034     integer ica, ira, iexc;
00035     extern logical lsame_(char *, char *);
00036     extern /* Subroutine */ int cswap_(integer *, complex *, integer *, 
00037             complex *, integer *);
00038     real sfmin1, sfmin2, sfmax1, sfmax2;
00039     extern integer icamax_(integer *, complex *, integer *);
00040     extern doublereal slamch_(char *);
00041     extern /* Subroutine */ int csscal_(integer *, real *, complex *, integer 
00042             *), xerbla_(char *, integer *);
00043     logical noconv;
00044 
00045 
00046 /*  -- LAPACK routine (version 3.2) -- */
00047 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00048 /*     November 2006 */
00049 
00050 /*     .. Scalar Arguments .. */
00051 /*     .. */
00052 /*     .. Array Arguments .. */
00053 /*     .. */
00054 
00055 /*  Purpose */
00056 /*  ======= */
00057 
00058 /*  CGEBAL balances a general complex matrix A.  This involves, first, */
00059 /*  permuting A by a similarity transformation to isolate eigenvalues */
00060 /*  in the first 1 to ILO-1 and last IHI+1 to N elements on the */
00061 /*  diagonal; and second, applying a diagonal similarity transformation */
00062 /*  to rows and columns ILO to IHI to make the rows and columns as */
00063 /*  close in norm as possible.  Both steps are optional. */
00064 
00065 /*  Balancing may reduce the 1-norm of the matrix, and improve the */
00066 /*  accuracy of the computed eigenvalues and/or eigenvectors. */
00067 
00068 /*  Arguments */
00069 /*  ========= */
00070 
00071 /*  JOB     (input) CHARACTER*1 */
00072 /*          Specifies the operations to be performed on A: */
00073 /*          = 'N':  none:  simply set ILO = 1, IHI = N, SCALE(I) = 1.0 */
00074 /*                  for i = 1,...,N; */
00075 /*          = 'P':  permute only; */
00076 /*          = 'S':  scale only; */
00077 /*          = 'B':  both permute and scale. */
00078 
00079 /*  N       (input) INTEGER */
00080 /*          The order of the matrix A.  N >= 0. */
00081 
00082 /*  A       (input/output) COMPLEX array, dimension (LDA,N) */
00083 /*          On entry, the input matrix A. */
00084 /*          On exit,  A is overwritten by the balanced matrix. */
00085 /*          If JOB = 'N', A is not referenced. */
00086 /*          See Further Details. */
00087 
00088 /*  LDA     (input) INTEGER */
00089 /*          The leading dimension of the array A.  LDA >= max(1,N). */
00090 
00091 /*  ILO     (output) INTEGER */
00092 /*  IHI     (output) INTEGER */
00093 /*          ILO and IHI are set to integers such that on exit */
00094 /*          A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N. */
00095 /*          If JOB = 'N' or 'S', ILO = 1 and IHI = N. */
00096 
00097 /*  SCALE   (output) REAL array, dimension (N) */
00098 /*          Details of the permutations and scaling factors applied to */
00099 /*          A.  If P(j) is the index of the row and column interchanged */
00100 /*          with row and column j and D(j) is the scaling factor */
00101 /*          applied to row and column j, then */
00102 /*          SCALE(j) = P(j)    for j = 1,...,ILO-1 */
00103 /*                   = D(j)    for j = ILO,...,IHI */
00104 /*                   = P(j)    for j = IHI+1,...,N. */
00105 /*          The order in which the interchanges are made is N to IHI+1, */
00106 /*          then 1 to ILO-1. */
00107 
00108 /*  INFO    (output) INTEGER */
00109 /*          = 0:  successful exit. */
00110 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
00111 
00112 /*  Further Details */
00113 /*  =============== */
00114 
00115 /*  The permutations consist of row and column interchanges which put */
00116 /*  the matrix in the form */
00117 
00118 /*             ( T1   X   Y  ) */
00119 /*     P A P = (  0   B   Z  ) */
00120 /*             (  0   0   T2 ) */
00121 
00122 /*  where T1 and T2 are upper triangular matrices whose eigenvalues lie */
00123 /*  along the diagonal.  The column indices ILO and IHI mark the starting */
00124 /*  and ending columns of the submatrix B. Balancing consists of applying */
00125 /*  a diagonal similarity transformation inv(D) * B * D to make the */
00126 /*  1-norms of each row of B and its corresponding column nearly equal. */
00127 /*  The output matrix is */
00128 
00129 /*     ( T1     X*D          Y    ) */
00130 /*     (  0  inv(D)*B*D  inv(D)*Z ). */
00131 /*     (  0      0           T2   ) */
00132 
00133 /*  Information about the permutations P and the diagonal matrix D is */
00134 /*  returned in the vector SCALE. */
00135 
00136 /*  This subroutine is based on the EISPACK routine CBAL. */
00137 
00138 /*  Modified by Tzu-Yi Chen, Computer Science Division, University of */
00139 /*    California at Berkeley, USA */
00140 
00141 /*  ===================================================================== */
00142 
00143 /*     .. Parameters .. */
00144 /*     .. */
00145 /*     .. Local Scalars .. */
00146 /*     .. */
00147 /*     .. External Functions .. */
00148 /*     .. */
00149 /*     .. External Subroutines .. */
00150 /*     .. */
00151 /*     .. Intrinsic Functions .. */
00152 /*     .. */
00153 /*     .. Statement Functions .. */
00154 /*     .. */
00155 /*     .. Statement Function definitions .. */
00156 /*     .. */
00157 /*     .. Executable Statements .. */
00158 
00159 /*     Test the input parameters */
00160 
00161     /* Parameter adjustments */
00162     a_dim1 = *lda;
00163     a_offset = 1 + a_dim1;
00164     a -= a_offset;
00165     --scale;
00166 
00167     /* Function Body */
00168     *info = 0;
00169     if (! lsame_(job, "N") && ! lsame_(job, "P") && ! lsame_(job, "S") 
00170             && ! lsame_(job, "B")) {
00171         *info = -1;
00172     } else if (*n < 0) {
00173         *info = -2;
00174     } else if (*lda < max(1,*n)) {
00175         *info = -4;
00176     }
00177     if (*info != 0) {
00178         i__1 = -(*info);
00179         xerbla_("CGEBAL", &i__1);
00180         return 0;
00181     }
00182 
00183     k = 1;
00184     l = *n;
00185 
00186     if (*n == 0) {
00187         goto L210;
00188     }
00189 
00190     if (lsame_(job, "N")) {
00191         i__1 = *n;
00192         for (i__ = 1; i__ <= i__1; ++i__) {
00193             scale[i__] = 1.f;
00194 /* L10: */
00195         }
00196         goto L210;
00197     }
00198 
00199     if (lsame_(job, "S")) {
00200         goto L120;
00201     }
00202 
00203 /*     Permutation to isolate eigenvalues if possible */
00204 
00205     goto L50;
00206 
00207 /*     Row and column exchange. */
00208 
00209 L20:
00210     scale[m] = (real) j;
00211     if (j == m) {
00212         goto L30;
00213     }
00214 
00215     cswap_(&l, &a[j * a_dim1 + 1], &c__1, &a[m * a_dim1 + 1], &c__1);
00216     i__1 = *n - k + 1;
00217     cswap_(&i__1, &a[j + k * a_dim1], lda, &a[m + k * a_dim1], lda);
00218 
00219 L30:
00220     switch (iexc) {
00221         case 1:  goto L40;
00222         case 2:  goto L80;
00223     }
00224 
00225 /*     Search for rows isolating an eigenvalue and push them down. */
00226 
00227 L40:
00228     if (l == 1) {
00229         goto L210;
00230     }
00231     --l;
00232 
00233 L50:
00234     for (j = l; j >= 1; --j) {
00235 
00236         i__1 = l;
00237         for (i__ = 1; i__ <= i__1; ++i__) {
00238             if (i__ == j) {
00239                 goto L60;
00240             }
00241             i__2 = j + i__ * a_dim1;
00242             if (a[i__2].r != 0.f || r_imag(&a[j + i__ * a_dim1]) != 0.f) {
00243                 goto L70;
00244             }
00245 L60:
00246             ;
00247         }
00248 
00249         m = l;
00250         iexc = 1;
00251         goto L20;
00252 L70:
00253         ;
00254     }
00255 
00256     goto L90;
00257 
00258 /*     Search for columns isolating an eigenvalue and push them left. */
00259 
00260 L80:
00261     ++k;
00262 
00263 L90:
00264     i__1 = l;
00265     for (j = k; j <= i__1; ++j) {
00266 
00267         i__2 = l;
00268         for (i__ = k; i__ <= i__2; ++i__) {
00269             if (i__ == j) {
00270                 goto L100;
00271             }
00272             i__3 = i__ + j * a_dim1;
00273             if (a[i__3].r != 0.f || r_imag(&a[i__ + j * a_dim1]) != 0.f) {
00274                 goto L110;
00275             }
00276 L100:
00277             ;
00278         }
00279 
00280         m = k;
00281         iexc = 2;
00282         goto L20;
00283 L110:
00284         ;
00285     }
00286 
00287 L120:
00288     i__1 = l;
00289     for (i__ = k; i__ <= i__1; ++i__) {
00290         scale[i__] = 1.f;
00291 /* L130: */
00292     }
00293 
00294     if (lsame_(job, "P")) {
00295         goto L210;
00296     }
00297 
00298 /*     Balance the submatrix in rows K to L. */
00299 
00300 /*     Iterative loop for norm reduction */
00301 
00302     sfmin1 = slamch_("S") / slamch_("P");
00303     sfmax1 = 1.f / sfmin1;
00304     sfmin2 = sfmin1 * 2.f;
00305     sfmax2 = 1.f / sfmin2;
00306 L140:
00307     noconv = FALSE_;
00308 
00309     i__1 = l;
00310     for (i__ = k; i__ <= i__1; ++i__) {
00311         c__ = 0.f;
00312         r__ = 0.f;
00313 
00314         i__2 = l;
00315         for (j = k; j <= i__2; ++j) {
00316             if (j == i__) {
00317                 goto L150;
00318             }
00319             i__3 = j + i__ * a_dim1;
00320             c__ += (r__1 = a[i__3].r, dabs(r__1)) + (r__2 = r_imag(&a[j + i__ 
00321                     * a_dim1]), dabs(r__2));
00322             i__3 = i__ + j * a_dim1;
00323             r__ += (r__1 = a[i__3].r, dabs(r__1)) + (r__2 = r_imag(&a[i__ + j 
00324                     * a_dim1]), dabs(r__2));
00325 L150:
00326             ;
00327         }
00328         ica = icamax_(&l, &a[i__ * a_dim1 + 1], &c__1);
00329         ca = c_abs(&a[ica + i__ * a_dim1]);
00330         i__2 = *n - k + 1;
00331         ira = icamax_(&i__2, &a[i__ + k * a_dim1], lda);
00332         ra = c_abs(&a[i__ + (ira + k - 1) * a_dim1]);
00333 
00334 /*        Guard against zero C or R due to underflow. */
00335 
00336         if (c__ == 0.f || r__ == 0.f) {
00337             goto L200;
00338         }
00339         g = r__ / 2.f;
00340         f = 1.f;
00341         s = c__ + r__;
00342 L160:
00343 /* Computing MAX */
00344         r__1 = max(f,c__);
00345 /* Computing MIN */
00346         r__2 = min(r__,g);
00347         if (c__ >= g || dmax(r__1,ca) >= sfmax2 || dmin(r__2,ra) <= sfmin2) {
00348             goto L170;
00349         }
00350         f *= 2.f;
00351         c__ *= 2.f;
00352         ca *= 2.f;
00353         r__ /= 2.f;
00354         g /= 2.f;
00355         ra /= 2.f;
00356         goto L160;
00357 
00358 L170:
00359         g = c__ / 2.f;
00360 L180:
00361 /* Computing MIN */
00362         r__1 = min(f,c__), r__1 = min(r__1,g);
00363         if (g < r__ || dmax(r__,ra) >= sfmax2 || dmin(r__1,ca) <= sfmin2) {
00364             goto L190;
00365         }
00366         f /= 2.f;
00367         c__ /= 2.f;
00368         g /= 2.f;
00369         ca /= 2.f;
00370         r__ *= 2.f;
00371         ra *= 2.f;
00372         goto L180;
00373 
00374 /*        Now balance. */
00375 
00376 L190:
00377         if (c__ + r__ >= s * .95f) {
00378             goto L200;
00379         }
00380         if (f < 1.f && scale[i__] < 1.f) {
00381             if (f * scale[i__] <= sfmin1) {
00382                 goto L200;
00383             }
00384         }
00385         if (f > 1.f && scale[i__] > 1.f) {
00386             if (scale[i__] >= sfmax1 / f) {
00387                 goto L200;
00388             }
00389         }
00390         g = 1.f / f;
00391         scale[i__] *= f;
00392         noconv = TRUE_;
00393 
00394         i__2 = *n - k + 1;
00395         csscal_(&i__2, &g, &a[i__ + k * a_dim1], lda);
00396         csscal_(&l, &f, &a[i__ * a_dim1 + 1], &c__1);
00397 
00398 L200:
00399         ;
00400     }
00401 
00402     if (noconv) {
00403         goto L140;
00404     }
00405 
00406 L210:
00407     *ilo = k;
00408     *ihi = l;
00409 
00410     return 0;
00411 
00412 /*     End of CGEBAL */
00413 
00414 } /* cgebal_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:26