dggbal.c

Go to the documentation of this file.

/* dggbal.f -- translated by f2c (version 20061008).
   You must link the resulting object file with libf2c:
        on Microsoft Windows system, link with libf2c.lib;
        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
*/

#include "f2c.h"
#include "blaswrap.h"

/* Table of constant values */

static integer c__1 = 1;
static doublereal c_b35 = 10.;
static doublereal c_b71 = .5;

/* Subroutine */ int dggbal_(char *job, integer *n, doublereal *a, integer *
        lda, doublereal *b, integer *ldb, integer *ilo, integer *ihi, 
        doublereal *lscale, doublereal *rscale, doublereal *work, integer *
        info)
{
    /* System generated locals */
    integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
    doublereal d__1, d__2, d__3;

    /* Builtin functions */
    double d_lg10(doublereal *), d_sign(doublereal *, doublereal *), pow_di(
            doublereal *, integer *);

    /* Local variables */
    integer i__, j, k, l, m;
    doublereal t;
    integer jc;
    doublereal ta, tb, tc;
    integer ir;
    doublereal ew;
    integer it, nr, ip1, jp1, lm1;
    doublereal cab, rab, ewc, cor, sum;
    integer nrp2, icab, lcab;
    doublereal beta, coef;
    integer irab, lrab;
    doublereal basl, cmax;
    extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, 
            integer *);
    doublereal coef2, coef5, gamma, alpha;
    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 
            integer *);
    extern logical lsame_(char *, char *);
    doublereal sfmin, sfmax;
    extern /* Subroutine */ int dswap_(integer *, doublereal *, integer *, 
            doublereal *, integer *);
    integer iflow;
    extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *, 
            integer *, doublereal *, integer *);
    integer kount;
    extern doublereal dlamch_(char *);
    doublereal pgamma;
    extern integer idamax_(integer *, doublereal *, integer *);
    extern /* Subroutine */ int xerbla_(char *, integer *);
    integer lsfmin, lsfmax;


/*  -- LAPACK routine (version 3.2) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/*     November 2006 */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */

/*  Purpose */
/*  ======= */

/*  DGGBAL balances a pair of general real matrices (A,B).  This */
/*  involves, first, permuting A and B by similarity transformations to */
/*  isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N */
/*  elements on the diagonal; and second, applying a diagonal similarity */
/*  transformation to rows and columns ILO to IHI to make the rows */
/*  and columns as close in norm as possible. Both steps are optional. */

/*  Balancing may reduce the 1-norm of the matrices, and improve the */
/*  accuracy of the computed eigenvalues and/or eigenvectors in the */
/*  generalized eigenvalue problem A*x = lambda*B*x. */

/*  Arguments */
/*  ========= */

/*  JOB     (input) CHARACTER*1 */
/*          Specifies the operations to be performed on A and B: */
/*          = 'N':  none:  simply set ILO = 1, IHI = N, LSCALE(I) = 1.0 */
/*                  and RSCALE(I) = 1.0 for i = 1,...,N. */
/*          = 'P':  permute only; */
/*          = 'S':  scale only; */
/*          = 'B':  both permute and scale. */

/*  N       (input) INTEGER */
/*          The order of the matrices A and B.  N >= 0. */

/*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
/*          On entry, the input matrix A. */
/*          On exit,  A is overwritten by the balanced matrix. */
/*          If JOB = 'N', A is not referenced. */

/*  LDA     (input) INTEGER */
/*          The leading dimension of the array A. LDA >= max(1,N). */

/*  B       (input/output) DOUBLE PRECISION array, dimension (LDB,N) */
/*          On entry, the input matrix B. */
/*          On exit,  B is overwritten by the balanced matrix. */
/*          If JOB = 'N', B is not referenced. */

/*  LDB     (input) INTEGER */
/*          The leading dimension of the array B. LDB >= max(1,N). */

/*  ILO     (output) INTEGER */
/*  IHI     (output) INTEGER */
/*          ILO and IHI are set to integers such that on exit */
/*          A(i,j) = 0 and B(i,j) = 0 if i > j and */
/*          j = 1,...,ILO-1 or i = IHI+1,...,N. */
/*          If JOB = 'N' or 'S', ILO = 1 and IHI = N. */

/*  LSCALE  (output) DOUBLE PRECISION array, dimension (N) */
/*          Details of the permutations and scaling factors applied */
/*          to the left side of A and B.  If P(j) is the index of the */
/*          row interchanged with row j, and D(j) */
/*          is the scaling factor applied to row j, then */
/*            LSCALE(j) = P(j)    for J = 1,...,ILO-1 */
/*                      = D(j)    for J = ILO,...,IHI */
/*                      = P(j)    for J = IHI+1,...,N. */
/*          The order in which the interchanges are made is N to IHI+1, */
/*          then 1 to ILO-1. */

/*  RSCALE  (output) DOUBLE PRECISION array, dimension (N) */
/*          Details of the permutations and scaling factors applied */
/*          to the right side of A and B.  If P(j) is the index of the */
/*          column interchanged with column j, and D(j) */
/*          is the scaling factor applied to column j, then */
/*            LSCALE(j) = P(j)    for J = 1,...,ILO-1 */
/*                      = D(j)    for J = ILO,...,IHI */
/*                      = P(j)    for J = IHI+1,...,N. */
/*          The order in which the interchanges are made is N to IHI+1, */
/*          then 1 to ILO-1. */

/*  WORK    (workspace) REAL array, dimension (lwork) */
/*          lwork must be at least max(1,6*N) when JOB = 'S' or 'B', and */
/*          at least 1 when JOB = 'N' or 'P'. */

/*  INFO    (output) INTEGER */
/*          = 0:  successful exit */
/*          < 0:  if INFO = -i, the i-th argument had an illegal value. */

/*  Further Details */
/*  =============== */

/*  See R.C. WARD, Balancing the generalized eigenvalue problem, */
/*                 SIAM J. Sci. Stat. Comp. 2 (1981), 141-152. */

/*  ===================================================================== */

/*     .. Parameters .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. External Functions .. */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */

/*     Test the input parameters */

    /* Parameter adjustments */
    a_dim1 = *lda;
    a_offset = 1 + a_dim1;
    a -= a_offset;
    b_dim1 = *ldb;
    b_offset = 1 + b_dim1;
    b -= b_offset;
    --lscale;
    --rscale;
    --work;

    /* Function Body */
    *info = 0;
    if (! lsame_(job, "N") && ! lsame_(job, "P") && ! lsame_(job, "S") 
            && ! lsame_(job, "B")) {
        *info = -1;
    } else if (*n < 0) {
        *info = -2;
    } else if (*lda < max(1,*n)) {
        *info = -4;
    } else if (*ldb < max(1,*n)) {
        *info = -6;
    }
    if (*info != 0) {
        i__1 = -(*info);
        xerbla_("DGGBAL", &i__1);
        return 0;
    }

/*     Quick return if possible */

    if (*n == 0) {
        *ilo = 1;
        *ihi = *n;
        return 0;
    }

    if (*n == 1) {
        *ilo = 1;
        *ihi = *n;
        lscale[1] = 1.;
        rscale[1] = 1.;
        return 0;
    }

    if (lsame_(job, "N")) {
        *ilo = 1;
        *ihi = *n;
        i__1 = *n;
        for (i__ = 1; i__ <= i__1; ++i__) {
            lscale[i__] = 1.;
            rscale[i__] = 1.;
/* L10: */
        }
        return 0;
    }

    k = 1;
    l = *n;
    if (lsame_(job, "S")) {
        goto L190;
    }

    goto L30;

/*     Permute the matrices A and B to isolate the eigenvalues. */

/*     Find row with one nonzero in columns 1 through L */

L20:
    l = lm1;
    if (l != 1) {
        goto L30;
    }

    rscale[1] = 1.;
    lscale[1] = 1.;
    goto L190;

L30:
    lm1 = l - 1;
    for (i__ = l; i__ >= 1; --i__) {
        i__1 = lm1;
        for (j = 1; j <= i__1; ++j) {
            jp1 = j + 1;
            if (a[i__ + j * a_dim1] != 0. || b[i__ + j * b_dim1] != 0.) {
                goto L50;
            }
/* L40: */
        }
        j = l;
        goto L70;

L50:
        i__1 = l;
        for (j = jp1; j <= i__1; ++j) {
            if (a[i__ + j * a_dim1] != 0. || b[i__ + j * b_dim1] != 0.) {
                goto L80;
            }
/* L60: */
        }
        j = jp1 - 1;

L70:
        m = l;
        iflow = 1;
        goto L160;
L80:
        ;
    }
    goto L100;

/*     Find column with one nonzero in rows K through N */

L90:
    ++k;

L100:
    i__1 = l;
    for (j = k; j <= i__1; ++j) {
        i__2 = lm1;
        for (i__ = k; i__ <= i__2; ++i__) {
            ip1 = i__ + 1;
            if (a[i__ + j * a_dim1] != 0. || b[i__ + j * b_dim1] != 0.) {
                goto L120;
            }
/* L110: */
        }
        i__ = l;
        goto L140;
L120:
        i__2 = l;
        for (i__ = ip1; i__ <= i__2; ++i__) {
            if (a[i__ + j * a_dim1] != 0. || b[i__ + j * b_dim1] != 0.) {
                goto L150;
            }
/* L130: */
        }
        i__ = ip1 - 1;
L140:
        m = k;
        iflow = 2;
        goto L160;
L150:
        ;
    }
    goto L190;

/*     Permute rows M and I */

L160:
    lscale[m] = (doublereal) i__;
    if (i__ == m) {
        goto L170;
    }
    i__1 = *n - k + 1;
    dswap_(&i__1, &a[i__ + k * a_dim1], lda, &a[m + k * a_dim1], lda);
    i__1 = *n - k + 1;
    dswap_(&i__1, &b[i__ + k * b_dim1], ldb, &b[m + k * b_dim1], ldb);

/*     Permute columns M and J */

L170:
    rscale[m] = (doublereal) j;
    if (j == m) {
        goto L180;
    }
    dswap_(&l, &a[j * a_dim1 + 1], &c__1, &a[m * a_dim1 + 1], &c__1);
    dswap_(&l, &b[j * b_dim1 + 1], &c__1, &b[m * b_dim1 + 1], &c__1);

L180:
    switch (iflow) {
        case 1:  goto L20;
        case 2:  goto L90;
    }

L190:
    *ilo = k;
    *ihi = l;

    if (lsame_(job, "P")) {
        i__1 = *ihi;
        for (i__ = *ilo; i__ <= i__1; ++i__) {
            lscale[i__] = 1.;
            rscale[i__] = 1.;
/* L195: */
        }
        return 0;
    }

    if (*ilo == *ihi) {
        return 0;
    }

/*     Balance the submatrix in rows ILO to IHI. */

    nr = *ihi - *ilo + 1;
    i__1 = *ihi;
    for (i__ = *ilo; i__ <= i__1; ++i__) {
        rscale[i__] = 0.;
        lscale[i__] = 0.;

        work[i__] = 0.;
        work[i__ + *n] = 0.;
        work[i__ + (*n << 1)] = 0.;
        work[i__ + *n * 3] = 0.;
        work[i__ + (*n << 2)] = 0.;
        work[i__ + *n * 5] = 0.;
/* L200: */
    }

/*     Compute right side vector in resulting linear equations */

    basl = d_lg10(&c_b35);
    i__1 = *ihi;
    for (i__ = *ilo; i__ <= i__1; ++i__) {
        i__2 = *ihi;
        for (j = *ilo; j <= i__2; ++j) {
            tb = b[i__ + j * b_dim1];
            ta = a[i__ + j * a_dim1];
            if (ta == 0.) {
                goto L210;
            }
            d__1 = abs(ta);
            ta = d_lg10(&d__1) / basl;
L210:
            if (tb == 0.) {
                goto L220;
            }
            d__1 = abs(tb);
            tb = d_lg10(&d__1) / basl;
L220:
            work[i__ + (*n << 2)] = work[i__ + (*n << 2)] - ta - tb;
            work[j + *n * 5] = work[j + *n * 5] - ta - tb;
/* L230: */
        }
/* L240: */
    }

    coef = 1. / (doublereal) (nr << 1);
    coef2 = coef * coef;
    coef5 = coef2 * .5;
    nrp2 = nr + 2;
    beta = 0.;
    it = 1;

/*     Start generalized conjugate gradient iteration */

L250:

    gamma = ddot_(&nr, &work[*ilo + (*n << 2)], &c__1, &work[*ilo + (*n << 2)]
, &c__1) + ddot_(&nr, &work[*ilo + *n * 5], &c__1, &work[*ilo + *
            n * 5], &c__1);

    ew = 0.;
    ewc = 0.;
    i__1 = *ihi;
    for (i__ = *ilo; i__ <= i__1; ++i__) {
        ew += work[i__ + (*n << 2)];
        ewc += work[i__ + *n * 5];
/* L260: */
    }

/* Computing 2nd power */
    d__1 = ew;
/* Computing 2nd power */
    d__2 = ewc;
/* Computing 2nd power */
    d__3 = ew - ewc;
    gamma = coef * gamma - coef2 * (d__1 * d__1 + d__2 * d__2) - coef5 * (
            d__3 * d__3);
    if (gamma == 0.) {
        goto L350;
    }
    if (it != 1) {
        beta = gamma / pgamma;
    }
    t = coef5 * (ewc - ew * 3.);
    tc = coef5 * (ew - ewc * 3.);

    dscal_(&nr, &beta, &work[*ilo], &c__1);
    dscal_(&nr, &beta, &work[*ilo + *n], &c__1);

    daxpy_(&nr, &coef, &work[*ilo + (*n << 2)], &c__1, &work[*ilo + *n], &
            c__1);
    daxpy_(&nr, &coef, &work[*ilo + *n * 5], &c__1, &work[*ilo], &c__1);

    i__1 = *ihi;
    for (i__ = *ilo; i__ <= i__1; ++i__) {
        work[i__] += tc;
        work[i__ + *n] += t;
/* L270: */
    }

/*     Apply matrix to vector */

    i__1 = *ihi;
    for (i__ = *ilo; i__ <= i__1; ++i__) {
        kount = 0;
        sum = 0.;
        i__2 = *ihi;
        for (j = *ilo; j <= i__2; ++j) {
            if (a[i__ + j * a_dim1] == 0.) {
                goto L280;
            }
            ++kount;
            sum += work[j];
L280:
            if (b[i__ + j * b_dim1] == 0.) {
                goto L290;
            }
            ++kount;
            sum += work[j];
L290:
            ;
        }
        work[i__ + (*n << 1)] = (doublereal) kount * work[i__ + *n] + sum;
/* L300: */
    }

    i__1 = *ihi;
    for (j = *ilo; j <= i__1; ++j) {
        kount = 0;
        sum = 0.;
        i__2 = *ihi;
        for (i__ = *ilo; i__ <= i__2; ++i__) {
            if (a[i__ + j * a_dim1] == 0.) {
                goto L310;
            }
            ++kount;
            sum += work[i__ + *n];
L310:
            if (b[i__ + j * b_dim1] == 0.) {
                goto L320;
            }
            ++kount;
            sum += work[i__ + *n];
L320:
            ;
        }
        work[j + *n * 3] = (doublereal) kount * work[j] + sum;
/* L330: */
    }

    sum = ddot_(&nr, &work[*ilo + *n], &c__1, &work[*ilo + (*n << 1)], &c__1) 
            + ddot_(&nr, &work[*ilo], &c__1, &work[*ilo + *n * 3], &c__1);
    alpha = gamma / sum;

/*     Determine correction to current iteration */

    cmax = 0.;
    i__1 = *ihi;
    for (i__ = *ilo; i__ <= i__1; ++i__) {
        cor = alpha * work[i__ + *n];
        if (abs(cor) > cmax) {
            cmax = abs(cor);
        }
        lscale[i__] += cor;
        cor = alpha * work[i__];
        if (abs(cor) > cmax) {
            cmax = abs(cor);
        }
        rscale[i__] += cor;
/* L340: */
    }
    if (cmax < .5) {
        goto L350;
    }

    d__1 = -alpha;
    daxpy_(&nr, &d__1, &work[*ilo + (*n << 1)], &c__1, &work[*ilo + (*n << 2)]
, &c__1);
    d__1 = -alpha;
    daxpy_(&nr, &d__1, &work[*ilo + *n * 3], &c__1, &work[*ilo + *n * 5], &
            c__1);

    pgamma = gamma;
    ++it;
    if (it <= nrp2) {
        goto L250;
    }

/*     End generalized conjugate gradient iteration */

L350:
    sfmin = dlamch_("S");
    sfmax = 1. / sfmin;
    lsfmin = (integer) (d_lg10(&sfmin) / basl + 1.);
    lsfmax = (integer) (d_lg10(&sfmax) / basl);
    i__1 = *ihi;
    for (i__ = *ilo; i__ <= i__1; ++i__) {
        i__2 = *n - *ilo + 1;
        irab = idamax_(&i__2, &a[i__ + *ilo * a_dim1], lda);
        rab = (d__1 = a[i__ + (irab + *ilo - 1) * a_dim1], abs(d__1));
        i__2 = *n - *ilo + 1;
        irab = idamax_(&i__2, &b[i__ + *ilo * b_dim1], ldb);
/* Computing MAX */
        d__2 = rab, d__3 = (d__1 = b[i__ + (irab + *ilo - 1) * b_dim1], abs(
                d__1));
        rab = max(d__2,d__3);
        d__1 = rab + sfmin;
        lrab = (integer) (d_lg10(&d__1) / basl + 1.);
        ir = (integer) (lscale[i__] + d_sign(&c_b71, &lscale[i__]));
/* Computing MIN */
        i__2 = max(ir,lsfmin), i__2 = min(i__2,lsfmax), i__3 = lsfmax - lrab;
        ir = min(i__2,i__3);
        lscale[i__] = pow_di(&c_b35, &ir);
        icab = idamax_(ihi, &a[i__ * a_dim1 + 1], &c__1);
        cab = (d__1 = a[icab + i__ * a_dim1], abs(d__1));
        icab = idamax_(ihi, &b[i__ * b_dim1 + 1], &c__1);
/* Computing MAX */
        d__2 = cab, d__3 = (d__1 = b[icab + i__ * b_dim1], abs(d__1));
        cab = max(d__2,d__3);
        d__1 = cab + sfmin;
        lcab = (integer) (d_lg10(&d__1) / basl + 1.);
        jc = (integer) (rscale[i__] + d_sign(&c_b71, &rscale[i__]));
/* Computing MIN */
        i__2 = max(jc,lsfmin), i__2 = min(i__2,lsfmax), i__3 = lsfmax - lcab;
        jc = min(i__2,i__3);
        rscale[i__] = pow_di(&c_b35, &jc);
/* L360: */
    }

/*     Row scaling of matrices A and B */

    i__1 = *ihi;
    for (i__ = *ilo; i__ <= i__1; ++i__) {
        i__2 = *n - *ilo + 1;
        dscal_(&i__2, &lscale[i__], &a[i__ + *ilo * a_dim1], lda);
        i__2 = *n - *ilo + 1;
        dscal_(&i__2, &lscale[i__], &b[i__ + *ilo * b_dim1], ldb);
/* L370: */
    }

/*     Column scaling of matrices A and B */

    i__1 = *ihi;
    for (j = *ilo; j <= i__1; ++j) {
        dscal_(ihi, &rscale[j], &a[j * a_dim1 + 1], &c__1);
        dscal_(ihi, &rscale[j], &b[j * b_dim1 + 1], &c__1);
/* L380: */
    }

    return 0;

/*     End of DGGBAL */

} /* dggbal_ */