zgegv.c

来自「算断裂的」· C语言 代码 · 共 705 行 · 第 1/2 页

	    return 0;
	}
    }

/*     Scale B */

    bnrm = zlange_("M", n, n, &B(1,1), ldb, &RWORK(1));
    bnrm1 = bnrm;
    bnrm2 = 1.;
    if (bnrm < 1.) {
	if (safmax * bnrm < 1.) {
	    bnrm1 = safmin;
	    bnrm2 = safmax * bnrm;
	}
    }

    if (bnrm > 0.) {
	zlascl_("G", &c_n1, &c_n1, &bnrm, &c_b16, n, n, &B(1,1), ldb, &
		iinfo);
	if (iinfo != 0) {
	    *info = *n + 10;
	    return 0;
	}
    }

/*     Permute the matrix to make it more nearly triangular   
       Also "balance" the matrix. */

    ileft = 1;
    iright = *n + 1;
    irwork = iright + *n;
    zggbal_("B", n, &A(1,1), lda, &B(1,1), ldb, &ilo, &ihi, &RWORK(
	    ileft), &RWORK(iright), &RWORK(irwork), &iinfo);
    if (iinfo != 0) {
	*info = *n + 1;
	goto L80;
    }

/*     Reduce B to triangular form, and initialize VL and/or VR */

    irows = ihi + 1 - ilo;
    if (ilv) {
	icols = *n + 1 - ilo;
    } else {
	icols = irows;
    }
    itau = 1;
    iwork = itau + irows;
    i__1 = *lwork + 1 - iwork;
    zgeqrf_(&irows, &icols, &B(ilo,ilo), ldb, &WORK(itau), &WORK(
	    iwork), &i__1, &iinfo);
    if (iinfo >= 0) {
/* Computing MAX */
	i__3 = iwork;
	i__1 = lwkopt, i__2 = (integer) WORK(iwork).r + iwork - 1;
	lwkopt = max(i__1,i__2);
    }
    if (iinfo != 0) {
	*info = *n + 2;
	goto L80;
    }

    i__1 = *lwork + 1 - iwork;
    zunmqr_("L", "C", &irows, &icols, &irows, &B(ilo,ilo), ldb, &
	    WORK(itau), &A(ilo,ilo), lda, &WORK(iwork), &i__1, &
	    iinfo);
    if (iinfo >= 0) {
/* Computing MAX */
	i__3 = iwork;
	i__1 = lwkopt, i__2 = (integer) WORK(iwork).r + iwork - 1;
	lwkopt = max(i__1,i__2);
    }
    if (iinfo != 0) {
	*info = *n + 3;
	goto L80;
    }

    if (ilvl) {
	zlaset_("Full", n, n, &c_b1, &c_b2, &VL(1,1), ldvl);
	i__1 = irows - 1;
	i__2 = irows - 1;
	zlacpy_("L", &i__1, &i__2, &B(ilo+1,ilo), ldb, &VL(ilo+1,ilo), ldvl);
	i__1 = *lwork + 1 - iwork;
	zungqr_(&irows, &irows, &irows, &VL(ilo,ilo), ldvl, &WORK(
		itau), &WORK(iwork), &i__1, &iinfo);
	if (iinfo >= 0) {
/* Computing MAX */
	    i__3 = iwork;
	    i__1 = lwkopt, i__2 = (integer) WORK(iwork).r + iwork - 1;
	    lwkopt = max(i__1,i__2);
	}
	if (iinfo != 0) {
	    *info = *n + 4;
	    goto L80;
	}
    }

    if (ilvr) {
	zlaset_("Full", n, n, &c_b1, &c_b2, &VR(1,1), ldvr);
    }

/*     Reduce to generalized Hessenberg form */

    if (ilv) {

/*        Eigenvectors requested -- work on whole matrix. */

	zgghrd_(jobvl, jobvr, n, &ilo, &ihi, &A(1,1), lda, &B(1,1), 
		ldb, &VL(1,1), ldvl, &VR(1,1), ldvr, &iinfo);
    } else {
	zgghrd_("N", "N", &irows, &c__1, &irows, &A(ilo,ilo), lda, 
		&B(ilo,ilo), ldb, &VL(1,1), ldvl, &VR(1,1), ldvr, &iinfo);
    }
    if (iinfo != 0) {
	*info = *n + 5;
	goto L80;
    }

/*     Perform QZ algorithm */

    iwork = itau;
    if (ilv) {
	*(unsigned char *)chtemp = 'S';
    } else {
	*(unsigned char *)chtemp = 'E';
    }
    i__1 = *lwork + 1 - iwork;
    zhgeqz_(chtemp, jobvl, jobvr, n, &ilo, &ihi, &A(1,1), lda, &B(1,1), ldb, &ALPHA(1), &BETA(1), &VL(1,1), ldvl, &VR(1,1), ldvr, &WORK(iwork), &i__1, &RWORK(irwork), &iinfo);
    if (iinfo >= 0) {
/* Computing MAX */
	i__3 = iwork;
	i__1 = lwkopt, i__2 = (integer) WORK(iwork).r + iwork - 1;
	lwkopt = max(i__1,i__2);
    }
    if (iinfo != 0) {
	if (iinfo > 0 && iinfo <= *n) {
	    *info = iinfo;
	} else if (iinfo > *n && iinfo <= *n << 1) {
	    *info = iinfo - *n;
	} else {
	    *info = *n + 6;
	}
	goto L80;
    }

    if (ilv) {

/*        Compute Eigenvectors */

	if (ilvl) {
	    if (ilvr) {
		*(unsigned char *)chtemp = 'B';
	    } else {
		*(unsigned char *)chtemp = 'L';
	    }
	} else {
	    *(unsigned char *)chtemp = 'R';
	}

	ztgevc_(chtemp, "B", ldumma, n, &A(1,1), lda, &B(1,1), ldb, 
		&VL(1,1), ldvl, &VR(1,1), ldvr, n, &in, &WORK(
		iwork), &RWORK(irwork), &iinfo);
	if (iinfo != 0) {
	    *info = *n + 7;
	    goto L80;
	}

/*        Undo balancing on VL and VR, rescale */

	if (ilvl) {
	    zggbak_("B", "L", n, &ilo, &ihi, &RWORK(ileft), &RWORK(iright), n,
		     &VL(1,1), ldvl, &iinfo);
	    if (iinfo != 0) {
		*info = *n + 8;
		goto L80;
	    }
	    i__1 = *n;
	    for (jc = 1; jc <= *n; ++jc) {
		temp = 0.;
		i__2 = *n;
		for (jr = 1; jr <= *n; ++jr) {
/* Computing MAX */
		    i__3 = jr + jc * vl_dim1;
		    d__3 = temp, d__4 = (d__1 = VL(jr,jc).r, abs(d__1)) + (
			    d__2 = d_imag(&VL(jr,jc)), abs(d__2));
		    temp = max(d__3,d__4);
/* L10: */
		}
		if (temp < safmin) {
		    goto L30;
		}
		temp = 1. / temp;
		i__2 = *n;
		for (jr = 1; jr <= *n; ++jr) {
		    i__3 = jr + jc * vl_dim1;
		    i__4 = jr + jc * vl_dim1;
		    z__1.r = temp * VL(jr,jc).r, z__1.i = temp * VL(jr,jc).i;
		    VL(jr,jc).r = z__1.r, VL(jr,jc).i = z__1.i;
/* L20: */
		}
L30:
		;
	    }
	}
	if (ilvr) {
	    zggbak_("B", "R", n, &ilo, &ihi, &RWORK(ileft), &RWORK(iright), n,
		     &VR(1,1), ldvr, &iinfo);
	    if (iinfo != 0) {
		*info = *n + 9;
		goto L80;
	    }
	    i__1 = *n;
	    for (jc = 1; jc <= *n; ++jc) {
		temp = 0.;
		i__2 = *n;
		for (jr = 1; jr <= *n; ++jr) {
/* Computing MAX */
		    i__3 = jr + jc * vr_dim1;
		    d__3 = temp, d__4 = (d__1 = VR(jr,jc).r, abs(d__1)) + (
			    d__2 = d_imag(&VR(jr,jc)), abs(d__2));
		    temp = max(d__3,d__4);
/* L40: */
		}
		if (temp < safmin) {
		    goto L60;
		}
		temp = 1. / temp;
		i__2 = *n;
		for (jr = 1; jr <= *n; ++jr) {
		    i__3 = jr + jc * vr_dim1;
		    i__4 = jr + jc * vr_dim1;
		    z__1.r = temp * VR(jr,jc).r, z__1.i = temp * VR(jr,jc).i;
		    VR(jr,jc).r = z__1.r, VR(jr,jc).i = z__1.i;
/* L50: */
		}
L60:
		;
	    }
	}

/*        End of eigenvector calculation */

    }

/*     Undo scaling in alpha, beta   

       Note: this does not give the alpha and beta for the unscaled   
       problem.   

       Un-scaling is limited to avoid underflow in alpha and beta   
       if they are significant. */

    i__1 = *n;
    for (jc = 1; jc <= *n; ++jc) {
	i__2 = jc;
	absar = (d__1 = ALPHA(jc).r, abs(d__1));
	absai = (d__1 = d_imag(&ALPHA(jc)), abs(d__1));
	i__2 = jc;
	absb = (d__1 = BETA(jc).r, abs(d__1));
	i__2 = jc;
	salfar = anrm * ALPHA(jc).r;
	salfai = anrm * d_imag(&ALPHA(jc));
	i__2 = jc;
	sbeta = bnrm * BETA(jc).r;
	ilimit = FALSE_;
	scale = 1.;

/*        Check for significant underflow in imaginary part of ALPHA 
  

   Computing MAX */
	d__1 = safmin, d__2 = eps * absar, d__1 = max(d__1,d__2), d__2 = eps *
		 absb;
	if (abs(salfai) < safmin && absai >= max(d__1,d__2)) {
	    ilimit = TRUE_;
/* Computing MAX */
	    d__1 = safmin, d__2 = anrm2 * absai;
	    scale = safmin / anrm1 / max(d__1,d__2);
	}

/*        Check for significant underflow in real part of ALPHA   

   Computing MAX */
	d__1 = safmin, d__2 = eps * absai, d__1 = max(d__1,d__2), d__2 = eps *
		 absb;
	if (abs(salfar) < safmin && absar >= max(d__1,d__2)) {
	    ilimit = TRUE_;
/* Computing MAX   
   Computing MAX */
	    d__3 = safmin, d__4 = anrm2 * absar;
	    d__1 = scale, d__2 = safmin / anrm1 / max(d__3,d__4);
	    scale = max(d__1,d__2);
	}

/*        Check for significant underflow in BETA   

   Computing MAX */
	d__1 = safmin, d__2 = eps * absar, d__1 = max(d__1,d__2), d__2 = eps *
		 absai;
	if (abs(sbeta) < safmin && absb >= max(d__1,d__2)) {
	    ilimit = TRUE_;
/* Computing MAX   
   Computing MAX */
	    d__3 = safmin, d__4 = bnrm2 * absb;
	    d__1 = scale, d__2 = safmin / bnrm1 / max(d__3,d__4);
	    scale = max(d__1,d__2);
	}

/*        Check for possible overflow when limiting scaling */

	if (ilimit) {
/* Computing MAX */
	    d__1 = abs(salfar), d__2 = abs(salfai), d__1 = max(d__1,d__2), 
		    d__2 = abs(sbeta);
	    temp = scale * safmin * max(d__1,d__2);
	    if (temp > 1.) {
		scale /= temp;
	    }
	    if (scale < 1.) {
		ilimit = FALSE_;
	    }
	}

/*        Recompute un-scaled ALPHA, BETA if necessary. */

	if (ilimit) {
	    i__2 = jc;
	    salfar = scale * ALPHA(jc).r * anrm;
	    salfai = scale * d_imag(&ALPHA(jc)) * anrm;
	    i__2 = jc;
	    z__2.r = scale * BETA(jc).r, z__2.i = scale * BETA(jc).i;
	    z__1.r = bnrm * z__2.r, z__1.i = bnrm * z__2.i;
	    sbeta = z__1.r;
	}
	i__2 = jc;
	z__1.r = salfar, z__1.i = salfai;
	ALPHA(jc).r = z__1.r, ALPHA(jc).i = z__1.i;
	i__2 = jc;
	BETA(jc).r = sbeta, BETA(jc).i = 0.;
/* L70: */
    }

L80:
    WORK(1).r = (doublereal) lwkopt, WORK(1).i = 0.;

    return 0;

/*     End of ZGEGV */

} /* zgegv_ */