dtgsen.c

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

/*  bounded as */

/*   max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2)) */
/*   max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2)) */

/*  See LAPACK User's Guide section 4.11 or the following references */
/*  for more information. */

/*  Note that if the default method for computing the Frobenius-norm- */
/*  based estimate DIF is not wanted (see DLATDF), then the parameter */
/*  IDIFJB (see below) should be changed from 3 to 4 (routine DLATDF */
/*  (IJOB = 2 will be used)). See DTGSYL for more details. */

/*  Based on contributions by */
/*     Bo Kagstrom and Peter Poromaa, Department of Computing Science, */
/*     Umea University, S-901 87 Umea, Sweden. */

/*  References */
/*  ========== */

/*  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the */
/*      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in */
/*      M.S. Moonen et al (eds), Linear Algebra for Large Scale and */
/*      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. */

/*  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified */
/*      Eigenvalues of a Regular Matrix Pair (A, B) and Condition */
/*      Estimation: Theory, Algorithms and Software, */
/*      Report UMINF - 94.04, Department of Computing Science, Umea */
/*      University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working */
/*      Note 87. To appear in Numerical Algorithms, 1996. */

/*  [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software */
/*      for Solving the Generalized Sylvester Equation and Estimating the */
/*      Separation between Regular Matrix Pairs, Report UMINF - 93.23, */
/*      Department of Computing Science, Umea University, S-901 87 Umea, */
/*      Sweden, December 1993, Revised April 1994, Also as LAPACK Working */
/*      Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1, */
/*      1996. */

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

/*     .. Parameters .. */
/*<       INTEGER            IDIFJB >*/
/*<       PARAMETER          ( IDIFJB = 3 ) >*/
/*<       DOUBLE PRECISION   ZERO, ONE >*/
/*<       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 ) >*/
/*     .. */
/*     .. Local Scalars .. */
/*<    >*/
/*<    >*/
/*<       DOUBLE PRECISION   DSCALE, DSUM, EPS, RDSCAL, SMLNUM >*/
/*     .. */
/*     .. External Subroutines .. */
/*<    >*/
/*     .. */
/*     .. External Functions .. */
/*<       DOUBLE PRECISION   DLAMCH >*/
/*<       EXTERNAL           DLAMCH >*/
/*     .. */
/*     .. Intrinsic Functions .. */
/*<       INTRINSIC          MAX, SIGN, SQRT >*/
/*     .. */
/*     .. Executable Statements .. */

/*     Decode and test the input parameters */

/*<       INFO = 0 >*/
    /* Parameter adjustments */
    --select;
    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;
    q_dim1 = *ldq;
    q_offset = 1 + q_dim1;
    q -= q_offset;
    z_dim1 = *ldz;
    z_offset = 1 + z_dim1;
    z__ -= z_offset;
    --dif;
    --work;
    --iwork;

    /* Function Body */
    *info = 0;
/*<       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) >*/
    lquery = *lwork == -1 || *liwork == -1;

/*<       IF( IJOB.LT.0 .OR. IJOB.GT.5 ) THEN >*/
    if (*ijob < 0 || *ijob > 5) {
/*<          INFO = -1 >*/
        *info = -1;
/*<       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( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN >*/
    } else if (*ldq < 1 || (*wantq && *ldq < *n)) {
/*<          INFO = -14 >*/
        *info = -14;
/*<       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN >*/
    } else if (*ldz < 1 || (*wantz && *ldz < *n)) {
/*<          INFO = -16 >*/
        *info = -16;
/*<       END IF >*/
    }

/*<       IF( INFO.NE.0 ) THEN >*/
    if (*info != 0) {
/*<          CALL XERBLA( 'DTGSEN', -INFO ) >*/
        i__1 = -(*info);
        xerbla_("DTGSEN", &i__1, (ftnlen)6);
/*<          RETURN >*/
        return 0;
/*<       END IF >*/
    }

/*     Get machine constants */

/*<       EPS = DLAMCH( 'P' ) >*/
    eps = dlamch_("P", (ftnlen)1);
/*<       SMLNUM = DLAMCH( 'S' ) / EPS >*/
    smlnum = dlamch_("S", (ftnlen)1) / eps;
/*<       IERR = 0 >*/
    ierr = 0;

/*<       WANTP = IJOB.EQ.1 .OR. IJOB.GE.4 >*/
    wantp = *ijob == 1 || *ijob >= 4;
/*<       WANTD1 = IJOB.EQ.2 .OR. IJOB.EQ.4 >*/
    wantd1 = *ijob == 2 || *ijob == 4;
/*<       WANTD2 = IJOB.EQ.3 .OR. IJOB.EQ.5 >*/
    wantd2 = *ijob == 3 || *ijob == 5;
/*<       WANTD = WANTD1 .OR. WANTD2 >*/
    wantd = wantd1 || wantd2;

/*     Set M to the dimension of the specified pair of deflating */
/*     subspaces. */

/*<       M = 0 >*/
    *m = 0;
/*<       PAIR = .FALSE. >*/
    pair = FALSE_;
/*<       DO 10 K = 1, N >*/
    i__1 = *n;
    for (k = 1; k <= i__1; ++k) {
/*<          IF( PAIR ) THEN >*/
        if (pair) {
/*<             PAIR = .FALSE. >*/
            pair = FALSE_;
/*<          ELSE >*/
        } else {
/*<             IF( K.LT.N ) THEN >*/
            if (k < *n) {
/*<                IF( A( K+1, K ).EQ.ZERO ) THEN >*/
                if (a[k + 1 + k * a_dim1] == 0.) {
/*<    >*/
                    if (select[k]) {
                        ++(*m);
                    }
/*<                ELSE >*/
                } else {
/*<                   PAIR = .TRUE. >*/
                    pair = TRUE_;
/*<    >*/
                    if (select[k] || select[k + 1]) {
                        *m += 2;
                    }
/*<                END IF >*/
                }
/*<             ELSE >*/
            } else {
/*<    >*/
                if (select[*n]) {
                    ++(*m);
                }
/*<             END IF >*/
            }
/*<          END IF >*/
        }
/*<    10 CONTINUE >*/
/* L10: */
    }

/*<       IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN >*/
    if (*ijob == 1 || *ijob == 2 || *ijob == 4) {
/*<          LWMIN = MAX( 1, 4*N+16, 2*M*( N-M ) ) >*/
/* Computing MAX */
        i__1 = 1, i__2 = (*n << 2) + 16, i__1 = max(i__1,i__2), i__2 = (*m << 
                1) * (*n - *m);
        lwmin = max(i__1,i__2);
/*<          LIWMIN = MAX( 1, N+6 ) >*/
/* Computing MAX */
        i__1 = 1, i__2 = *n + 6;
        liwmin = max(i__1,i__2);
/*<       ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN >*/
    } else if (*ijob == 3 || *ijob == 5) {
/*<          LWMIN = MAX( 1, 4*N+16, 4*M*( N-M ) ) >*/
/* Computing MAX */
        i__1 = 1, i__2 = (*n << 2) + 16, i__1 = max(i__1,i__2), i__2 = (*m << 
                2) * (*n - *m);
        lwmin = max(i__1,i__2);
/*<          LIWMIN = MAX( 1, 2*M*( N-M ), N+6 ) >*/
/* Computing MAX */
        i__1 = 1, i__2 = (*m << 1) * (*n - *m), i__1 = max(i__1,i__2), i__2 = 
                *n + 6;
        liwmin = max(i__1,i__2);
/*<       ELSE >*/
    } else {
/*<          LWMIN = MAX( 1, 4*N+16 ) >*/
/* Computing MAX */
        i__1 = 1, i__2 = (*n << 2) + 16;
        lwmin = max(i__1,i__2);
/*<          LIWMIN = 1 >*/
        liwmin = 1;
/*<       END IF >*/
    }

/*<       WORK( 1 ) = LWMIN >*/
    work[1] = (doublereal) lwmin;
/*<       IWORK( 1 ) = LIWMIN >*/
    iwork[1] = liwmin;

/*<       IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN >*/
    if (*lwork < lwmin && ! lquery) {
/*<          INFO = -22 >*/
        *info = -22;
/*<       ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN >*/
    } else if (*liwork < liwmin && ! lquery) {
/*<          INFO = -24 >*/
        *info = -24;
/*<       END IF >*/
    }

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

/*     Quick return if possible. */

/*<       IF( M.EQ.N .OR. M.EQ.0 ) THEN >*/
    if (*m == *n || *m == 0) {
/*<          IF( WANTP ) THEN >*/
        if (wantp) {
/*<             PL = ONE >*/
            *pl = 1.;
/*<             PR = ONE >*/
            *pr = 1.;
/*<          END IF >*/
        }
/*<          IF( WANTD ) THEN >*/
        if (wantd) {
/*<             DSCALE = ZERO >*/
            dscale = 0.;
/*<             DSUM = ONE >*/
            dsum = 1.;
/*<             DO 20 I = 1, N >*/
            i__1 = *n;
            for (i__ = 1; i__ <= i__1; ++i__) {
/*<                CALL DLASSQ( N, A( 1, I ), 1, DSCALE, DSUM ) >*/
                dlassq_(n, &a[i__ * a_dim1 + 1], &c__1, &dscale, &dsum);
/*<                CALL DLASSQ( N, B( 1, I ), 1, DSCALE, DSUM ) >*/
                dlassq_(n, &b[i__ * b_dim1 + 1], &c__1, &dscale, &dsum);
/*<    20       CONTINUE >*/
/* L20: */
            }
/*<             DIF( 1 ) = DSCALE*SQRT( DSUM ) >*/
            dif[1] = dscale * sqrt(dsum);
/*<             DIF( 2 ) = DIF( 1 ) >*/
            dif[2] = dif[1];
/*<          END IF >*/
        }
/*<          GO TO 60 >*/
        goto L60;
/*<       END IF >*/
    }

/*     Collect the selected blocks at the top-left corner of (A, B). */

/*<       KS = 0 >*/
    ks = 0;
/*<       PAIR = .FALSE. >*/
    pair = FALSE_;
/*<       DO 30 K = 1, N >*/
    i__1 = *n;
    for (k = 1; k <= i__1; ++k) {
/*<          IF( PAIR ) THEN >*/
        if (pair) {
/*<             PAIR = .FALSE. >*/
            pair = FALSE_;
/*<          ELSE >*/
        } else {

/*<             SWAP = SELECT( K ) >*/
            swap = select[k];
/*<             IF( K.LT.N ) THEN >*/
            if (k < *n) {
/*<                IF( A( K+1, K ).NE.ZERO ) THEN >*/
                if (a[k + 1 + k * a_dim1] != 0.) {
/*<                   PAIR = .TRUE. >*/
                    pair = TRUE_;
/*<                   SWAP = SWAP .OR. SELECT( K+1 ) >*/
                    swap = swap || select[k + 1];
/*<                END IF >*/
                }
/*<             END IF >*/
            }

/*<             IF( SWAP ) THEN >*/
            if (swap) {
/*<                KS = KS + 1 >*/
                ++ks;

/*              Swap the K-th block to position KS. */
/*              Perform the reordering of diagonal blocks in (A, B) */
/*              by orthogonal transformation matrices and update */
/*              Q and Z accordingly (if requested): */

/*<                KK = K >*/
                kk = k;
/*<    >*/
                if (k != ks) {
                    dtgexc_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], 
                            ldb, &q[q_offset], ldq, &z__[z_offset], ldz, &kk, 
                            &ks, &work[1], lwork, &ierr);
                }

/*<                IF( IERR.GT.0 ) THEN >*/
                if (ierr > 0) {

/*                 Swap is rejected: exit. */

/*<                   INFO = 1 >*/
                    *info = 1;
/*<                   IF( WANTP ) THEN >*/
                    if (wantp) {