dgges.c

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

/*<       LOGICAL            LSAME >*/
/*<       INTEGER            ILAENV >*/
/*<       DOUBLE PRECISION   DLAMCH, DLANGE >*/
/*<       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE >*/
/*     .. */
/*     .. Intrinsic Functions .. */
/*<       INTRINSIC          ABS, MAX, SQRT >*/
/*     .. */
/*     .. Executable Statements .. */

/*     Decode the input arguments */

/*<       IF( LSAME( JOBVSL, 'N' ) ) THEN >*/
    /* Parameter adjustments */
    a_dim1 = *lda;
    a_offset = 1 + a_dim1;
    a -= a_offset;
    b_dim1 = *ldb;
    b_offset = 1 + b_dim1;
    b -= b_offset;
    --alphar;
    --alphai;
    --beta;
    vsl_dim1 = *ldvsl;
    vsl_offset = 1 + vsl_dim1;
    vsl -= vsl_offset;
    vsr_dim1 = *ldvsr;
    vsr_offset = 1 + vsr_dim1;
    vsr -= vsr_offset;
    --work;
    --bwork;

    /* Function Body */
    if (lsame_(jobvsl, "N", (ftnlen)1, (ftnlen)1)) {
/*<          IJOBVL = 1 >*/
        ijobvl = 1;
/*<          ILVSL = .FALSE. >*/
        ilvsl = FALSE_;
/*<       ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN >*/
    } else if (lsame_(jobvsl, "V", (ftnlen)1, (ftnlen)1)) {
/*<          IJOBVL = 2 >*/
        ijobvl = 2;
/*<          ILVSL = .TRUE. >*/
        ilvsl = TRUE_;
/*<       ELSE >*/
    } else {
/*<          IJOBVL = -1 >*/
        ijobvl = -1;
/*<          ILVSL = .FALSE. >*/
        ilvsl = FALSE_;
/*<       END IF >*/
    }

/*<       IF( LSAME( JOBVSR, 'N' ) ) THEN >*/
    if (lsame_(jobvsr, "N", (ftnlen)1, (ftnlen)1)) {
/*<          IJOBVR = 1 >*/
        ijobvr = 1;
/*<          ILVSR = .FALSE. >*/
        ilvsr = FALSE_;
/*<       ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN >*/
    } else if (lsame_(jobvsr, "V", (ftnlen)1, (ftnlen)1)) {
/*<          IJOBVR = 2 >*/
        ijobvr = 2;
/*<          ILVSR = .TRUE. >*/
        ilvsr = TRUE_;
/*<       ELSE >*/
    } else {
/*<          IJOBVR = -1 >*/
        ijobvr = -1;
/*<          ILVSR = .FALSE. >*/
        ilvsr = FALSE_;
/*<       END IF >*/
    }

/*<       WANTST = LSAME( SORT, 'S' ) >*/
    wantst = lsame_(sort, "S", (ftnlen)1, (ftnlen)1);

/*     Test the input arguments */

/*<       INFO = 0 >*/
    *info = 0;
/*<       LQUERY = ( LWORK.EQ.-1 ) >*/
    lquery = *lwork == -1;
/*<       IF( IJOBVL.LE.0 ) THEN >*/
    if (ijobvl <= 0) {
/*<          INFO = -1 >*/
        *info = -1;
/*<       ELSE IF( IJOBVR.LE.0 ) THEN >*/
    } else if (ijobvr <= 0) {
/*<          INFO = -2 >*/
        *info = -2;
/*<       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN >*/
    } else if (! wantst && ! lsame_(sort, "N", (ftnlen)1, (ftnlen)1)) {
/*<          INFO = -3 >*/
        *info = -3;
/*<       ELSE IF( N.LT.0 ) THEN >*/
    } else if (*n < 0) {
/*<          INFO = -5 >*/
        *info = -5;
/*<       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN >*/
    } else if (*lda < max(1,*n)) {
/*<          INFO = -7 >*/
        *info = -7;
/*<       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN >*/
    } else if (*ldb < max(1,*n)) {
/*<          INFO = -9 >*/
        *info = -9;
/*<       ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN >*/
    } else if (*ldvsl < 1 || (ilvsl && *ldvsl < *n)) {
/*<          INFO = -15 >*/
        *info = -15;
/*<       ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN >*/
    } else if (*ldvsr < 1 || (ilvsr && *ldvsr < *n)) {
/*<          INFO = -17 >*/
        *info = -17;
/*<       END IF >*/
    }

/*     Compute workspace */
/*      (Note: Comments in the code beginning "Workspace:" describe the */
/*       minimal amount of workspace needed at that point in the code, */
/*       as well as the preferred amount for good performance. */
/*       NB refers to the optimal block size for the immediately */
/*       following subroutine, as returned by ILAENV.) */

/*<       MINWRK = 1 >*/
    minwrk = 1;
/*<       IF( INFO.EQ.0 .AND. ( LWORK.GE.1 .OR. LQUERY ) ) THEN >*/
    if (*info == 0 && (*lwork >= 1 || lquery)) {
/*<          MINWRK = 7*( N+1 ) + 16 >*/
        minwrk = (*n + 1) * 7 + 16;
/*<    >*/
        maxwrk = (*n + 1) * 7 + *n * ilaenv_(&c__1, "DGEQRF", " ", n, &c__1, 
                n, &c__0, (ftnlen)6, (ftnlen)1) + 16;
/*<          IF( ILVSL ) THEN >*/
        if (ilvsl) {
/*<    >*/
/* Computing MAX */
            i__1 = maxwrk, i__2 = (*n + 1) * 7 + *n * ilaenv_(&c__1, "DORGQR",
                     " ", n, &c__1, n, &c_n1, (ftnlen)6, (ftnlen)1);
            maxwrk = max(i__1,i__2);
/*<          END IF >*/
        }
/*<          WORK( 1 ) = MAXWRK >*/
        work[1] = (doublereal) maxwrk;
/*<       END IF >*/
    }

/*<    >*/
    if (*lwork < minwrk && ! lquery) {
        *info = -19;
    }
/*<       IF( INFO.NE.0 ) THEN >*/
    if (*info != 0) {
/*<          CALL XERBLA( 'DGGES ', -INFO ) >*/
        i__1 = -(*info);
        xerbla_("DGGES ", &i__1, (ftnlen)6);
/*<          RETURN >*/
        return 0;
/*<       ELSE IF( LQUERY ) THEN >*/
    } else if (lquery) {
/*<          RETURN >*/
        return 0;
/*<       END IF >*/
    }

/*     Quick return if possible */

/*<       IF( N.EQ.0 ) THEN >*/
    if (*n == 0) {
/*<          SDIM = 0 >*/
        *sdim = 0;
/*<          RETURN >*/
        return 0;
/*<       END IF >*/
    }

/*     Get machine constants */

/*<       EPS = DLAMCH( 'P' ) >*/
    eps = dlamch_("P", (ftnlen)1);
/*<       SAFMIN = DLAMCH( 'S' ) >*/
    safmin = dlamch_("S", (ftnlen)1);
/*<       SAFMAX = ONE / SAFMIN >*/
    safmax = 1. / safmin;
/*<       CALL DLABAD( SAFMIN, SAFMAX ) >*/
    dlabad_(&safmin, &safmax);
/*<       SMLNUM = SQRT( SAFMIN ) / EPS >*/
    smlnum = sqrt(safmin) / eps;
/*<       BIGNUM = ONE / SMLNUM >*/
    bignum = 1. / smlnum;

/*     Scale A if max element outside range [SMLNUM,BIGNUM] */

/*<       ANRM = DLANGE( 'M', N, N, A, LDA, WORK ) >*/
    anrm = dlange_("M", n, n, &a[a_offset], lda, &work[1], (ftnlen)1);
/*<       ILASCL = .FALSE. >*/
    ilascl = FALSE_;
/*<       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN >*/
    if (anrm > 0. && anrm < smlnum) {
/*<          ANRMTO = SMLNUM >*/
        anrmto = smlnum;
/*<          ILASCL = .TRUE. >*/
        ilascl = TRUE_;
/*<       ELSE IF( ANRM.GT.BIGNUM ) THEN >*/
    } else if (anrm > bignum) {
/*<          ANRMTO = BIGNUM >*/
        anrmto = bignum;
/*<          ILASCL = .TRUE. >*/
        ilascl = TRUE_;
/*<       END IF >*/
    }
/*<    >*/
    if (ilascl) {
        dlascl_("G", &c__0, &c__0, &anrm, &anrmto, n, n, &a[a_offset], lda, &
                ierr, (ftnlen)1);
    }

/*     Scale B if max element outside range [SMLNUM,BIGNUM] */

/*<       BNRM = DLANGE( 'M', N, N, B, LDB, WORK ) >*/
    bnrm = dlange_("M", n, n, &b[b_offset], ldb, &work[1], (ftnlen)1);
/*<       ILBSCL = .FALSE. >*/
    ilbscl = FALSE_;
/*<       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN >*/
    if (bnrm > 0. && bnrm < smlnum) {
/*<          BNRMTO = SMLNUM >*/
        bnrmto = smlnum;
/*<          ILBSCL = .TRUE. >*/
        ilbscl = TRUE_;
/*<       ELSE IF( BNRM.GT.BIGNUM ) THEN >*/
    } else if (bnrm > bignum) {
/*<          BNRMTO = BIGNUM >*/
        bnrmto = bignum;
/*<          ILBSCL = .TRUE. >*/
        ilbscl = TRUE_;
/*<       END IF >*/
    }
/*<    >*/
    if (ilbscl) {
        dlascl_("G", &c__0, &c__0, &bnrm, &bnrmto, n, n, &b[b_offset], ldb, &
                ierr, (ftnlen)1);
    }

/*     Permute the matrix to make it more nearly triangular */
/*     (Workspace: need 6*N + 2*N space for storing balancing factors) */

/*<       ILEFT = 1 >*/
    ileft = 1;
/*<       IRIGHT = N + 1 >*/
    iright = *n + 1;
/*<       IWRK = IRIGHT + N >*/
    iwrk = iright + *n;
/*<    >*/
    dggbal_("P", n, &a[a_offset], lda, &b[b_offset], ldb, &ilo, &ihi, &work[
            ileft], &work[iright], &work[iwrk], &ierr, (ftnlen)1);

/*     Reduce B to triangular form (QR decomposition of B) */
/*     (Workspace: need N, prefer N*NB) */

/*<       IROWS = IHI + 1 - ILO >*/
    irows = ihi + 1 - ilo;
/*<       ICOLS = N + 1 - ILO >*/
    icols = *n + 1 - ilo;
/*<       ITAU = IWRK >*/
    itau = iwrk;
/*<       IWRK = ITAU + IROWS >*/
    iwrk = itau + irows;
/*<    >*/
    i__1 = *lwork + 1 - iwrk;
    dgeqrf_(&irows, &icols, &b[ilo + ilo * b_dim1], ldb, &work[itau], &work[
            iwrk], &i__1, &ierr);

/*     Apply the orthogonal transformation to matrix A */
/*     (Workspace: need N, prefer N*NB) */

/*<    >*/
    i__1 = *lwork + 1 - iwrk;
    dormqr_("L", "T", &irows, &icols, &irows, &b[ilo + ilo * b_dim1], ldb, &
            work[itau], &a[ilo + ilo * a_dim1], lda, &work[iwrk], &i__1, &
            ierr, (ftnlen)1, (ftnlen)1);

/*     Initialize VSL */
/*     (Workspace: need N, prefer N*NB) */

/*<       IF( ILVSL ) THEN >*/
    if (ilvsl) {
/*<          CALL DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL ) >*/
        dlaset_("Full", n, n, &c_b33, &c_b34, &vsl[vsl_offset], ldvsl, (
                ftnlen)4);
/*<    >*/
        i__1 = irows - 1;
        i__2 = irows - 1;
        dlacpy_("L", &i__1, &i__2, &b[ilo + 1 + ilo * b_dim1], ldb, &vsl[ilo 
                + 1 + ilo * vsl_dim1], ldvsl, (ftnlen)1);
/*<    >*/
        i__1 = *lwork + 1 - iwrk;
        dorgqr_(&irows, &irows, &irows, &vsl[ilo + ilo * vsl_dim1], ldvsl, &
                work[itau], &work[iwrk], &i__1, &ierr);
/*<       END IF >*/
    }

/*     Initialize VSR */