dggbal.c

来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 515 行 · 第 1/2 页

#include "f2c.h"
#include "netlib.h"

/* Table of constant values */
static integer c__1 = 1;
static doublereal c_b34 = 10.;
static doublereal c_b70 = .5;

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

    /* Local variables */
    static integer lcab;
    static doublereal beta, coef;
    static integer irab, lrab;
    static doublereal basl, cmax;
    static doublereal coef2, coef5;
    static integer i, j, k, l, m;
    static doublereal gamma, t, alpha;
    static doublereal sfmin, sfmax;
    static integer iflow;
    static integer kount, jc;
    static doublereal ta, tb, tc;
    static integer ir, it;
    static doublereal ew;
    static integer nr;
    static doublereal pgamma;
    static integer lsfmin, lsfmax, ip1, jp1, lm1;
    static doublereal cab, rab, ewc, cor, sum;
    static integer nrp2, icab;


/*  -- LAPACK routine (version 3.0) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
/*     Courant Institute, Argonne National Lab, and Rice University */
/*     September 30, 1994 */

/*  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) DOUBLE PRECISION array, dimension (6*N)           */
/*                                                                        */
/*  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.            */
/*                                                                        */
/*  ===================================================================== */

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

/*     Test the input parameters */

    *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 = -5;
    }
    if (*info != 0) {
        i__1 = -(*info);
        xerbla_("DGGBAL", &i__1);
        return;
    }

    k = 1;
    l = *n;

/*     Quick return if possible */

    if (*n == 0) {
        return;
    }

    if (lsame_(job, "N")) {
        *ilo = 1;
        *ihi = *n;
        for (i = 1; i <= *n; ++i) {
            lscale[i] = 1.;
            rscale[i] = 1.;
        }
        return;
    }

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

    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) {
        for (j = 1; j <= lm1; ++j) {
            jp1 = j + 1;
            if (a[i + j * a_dim1] != 0. || b[i + j * b_dim1] != 0.) {
                goto L50;
            }
        }
        j = l;
        goto L70;

L50:
        for (j = jp1; j <= l; ++j) {
            if (a[i + j * a_dim1] != 0. || b[i + j * b_dim1] != 0.) {
                goto L80;
            }
        }
        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:
    for (j = k; j <= l; ++j) {
        for (i = k; i <= lm1; ++i) {
            ip1 = i + 1;
            if (a[i + j * a_dim1] != 0. || b[i + j * b_dim1] != 0.) {
                goto L120;
            }
        }
        i = l;
        goto L140;
L120:
        for (i = ip1; i <= l; ++i) {