dtgevc.c

著名的LAPACK矩阵计算软件包, 是比较新的版本, 一般用到矩阵分解的朋友也许会用到

	    if (ilcplx) {
		ilcplx = FALSE_;
		goto L10;
	    }
	    if (j < *n) {
		if (a_ref(j + 1, j) != 0.) {
		    ilcplx = TRUE_;
		}
	    }
	    if (ilcplx) {
		if (select[j] || select[j + 1]) {
		    im += 2;
		}
	    } else {
		if (select[j]) {
		    ++im;
		}
	    }
L10:
	    ;
	}
    } else {
	im = *n;
    }

/*     Check 2-by-2 diagonal blocks of A, B */

    ilabad = FALSE_;
    ilbbad = FALSE_;
    i__1 = *n - 1;
    for (j = 1; j <= i__1; ++j) {
	if (a_ref(j + 1, j) != 0.) {
	    if (b_ref(j, j) == 0. || b_ref(j + 1, j + 1) == 0. || b_ref(j, j 
		    + 1) != 0.) {
		ilbbad = TRUE_;
	    }
	    if (j < *n - 1) {
		if (a_ref(j + 2, j + 1) != 0.) {
		    ilabad = TRUE_;
		}
	    }
	}
/* L20: */
    }

    if (ilabad) {
	*info = -5;
    } else if (ilbbad) {
	*info = -7;
    } else if (compl && *ldvl < *n || *ldvl < 1) {
	*info = -10;
    } else if (compr && *ldvr < *n || *ldvr < 1) {
	*info = -12;
    } else if (*mm < im) {
	*info = -13;
    }
    if (*info != 0) {
	i__1 = -(*info);
	xerbla_("DTGEVC", &i__1);
	return 0;
    }

/*     Quick return if possible */

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

/*     Machine Constants */

    safmin = dlamch_("Safe minimum");
    big = 1. / safmin;
    dlabad_(&safmin, &big);
    ulp = dlamch_("Epsilon") * dlamch_("Base");
    small = safmin * *n / ulp;
    big = 1. / small;
    bignum = 1. / (safmin * *n);

/*     Compute the 1-norm of each column of the strictly upper triangular   
       part (i.e., excluding all elements belonging to the diagonal   
       blocks) of A and B to check for possible overflow in the   
       triangular solver. */

    anorm = (d__1 = a_ref(1, 1), abs(d__1));
    if (*n > 1) {
	anorm += (d__1 = a_ref(2, 1), abs(d__1));
    }
    bnorm = (d__1 = b_ref(1, 1), abs(d__1));
    work[1] = 0.;
    work[*n + 1] = 0.;

    i__1 = *n;
    for (j = 2; j <= i__1; ++j) {
	temp = 0.;
	temp2 = 0.;
	if (a_ref(j, j - 1) == 0.) {
	    iend = j - 1;
	} else {
	    iend = j - 2;
	}
	i__2 = iend;
	for (i__ = 1; i__ <= i__2; ++i__) {
	    temp += (d__1 = a_ref(i__, j), abs(d__1));
	    temp2 += (d__1 = b_ref(i__, j), abs(d__1));
/* L30: */
	}
	work[j] = temp;
	work[*n + j] = temp2;
/* Computing MIN */
	i__3 = j + 1;
	i__2 = min(i__3,*n);
	for (i__ = iend + 1; i__ <= i__2; ++i__) {
	    temp += (d__1 = a_ref(i__, j), abs(d__1));
	    temp2 += (d__1 = b_ref(i__, j), abs(d__1));
/* L40: */
	}
	anorm = max(anorm,temp);
	bnorm = max(bnorm,temp2);
/* L50: */
    }

    ascale = 1. / max(anorm,safmin);
    bscale = 1. / max(bnorm,safmin);

/*     ---------------------- Begin Timing Code -------------------------   
   Computing 2nd power */
    i__1 = *n;
    latime_1.ops += (doublereal) (i__1 * i__1 + *n * 3 + 6);
/*     ----------------------- End Timing Code --------------------------   

       Left eigenvectors */

    if (compl) {
	ieig = 0;

/*        Main loop over eigenvalues */

	ilcplx = FALSE_;
	i__1 = *n;
	for (je = 1; je <= i__1; ++je) {

/*           Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or   
             (b) this would be the second of a complex pair.   
             Check for complex eigenvalue, so as to be sure of which   
             entry(-ies) of SELECT to look at. */

	    if (ilcplx) {
		ilcplx = FALSE_;
		goto L220;
	    }
	    nw = 1;
	    if (je < *n) {
		if (a_ref(je + 1, je) != 0.) {
		    ilcplx = TRUE_;
		    nw = 2;
		}
	    }
	    if (ilall) {
		ilcomp = TRUE_;
	    } else if (ilcplx) {
		ilcomp = select[je] || select[je + 1];
	    } else {
		ilcomp = select[je];
	    }
	    if (! ilcomp) {
		goto L220;
	    }

/*           Decide if (a) singular pencil, (b) real eigenvalue, or   
             (c) complex eigenvalue. */

	    if (! ilcplx) {
		if ((d__1 = a_ref(je, je), abs(d__1)) <= safmin && (d__2 = 
			b_ref(je, je), abs(d__2)) <= safmin) {

/*                 Singular matrix pencil -- returns unit eigenvector */

		    ++ieig;
		    i__2 = *n;
		    for (jr = 1; jr <= i__2; ++jr) {
			vl_ref(jr, ieig) = 0.;
/* L60: */
		    }
		    vl_ref(ieig, ieig) = 1.;
		    goto L220;
		}
	    }

/*           Clear vector */

	    i__2 = nw * *n;
	    for (jr = 1; jr <= i__2; ++jr) {
		work[(*n << 1) + jr] = 0.;
/* L70: */
	    }
/*                                                 T   
             Compute coefficients in  ( a A - b B )  y = 0   
                a  is  ACOEF   
                b  is  BCOEFR + i*BCOEFI */

	    if (! ilcplx) {

/*              Real eigenvalue   

   Computing MAX */
		d__3 = (d__1 = a_ref(je, je), abs(d__1)) * ascale, d__4 = (
			d__2 = b_ref(je, je), abs(d__2)) * bscale, d__3 = max(
			d__3,d__4);
		temp = 1. / max(d__3,safmin);
		salfar = temp * a_ref(je, je) * ascale;
		sbeta = temp * b_ref(je, je) * bscale;
		acoef = sbeta * ascale;
		bcoefr = salfar * bscale;
		bcoefi = 0.;

/*              Scale to avoid underflow */

		scale = 1.;
		lsa = abs(sbeta) >= safmin && abs(acoef) < small;
		lsb = abs(salfar) >= safmin && abs(bcoefr) < small;
		if (lsa) {
		    scale = small / abs(sbeta) * min(anorm,big);
		}
		if (lsb) {
/* Computing MAX */
		    d__1 = scale, d__2 = small / abs(salfar) * min(bnorm,big);
		    scale = max(d__1,d__2);
		}
		if (lsa || lsb) {
/* Computing MIN   
   Computing MAX */
		    d__3 = 1., d__4 = abs(acoef), d__3 = max(d__3,d__4), d__4 
			    = abs(bcoefr);
		    d__1 = scale, d__2 = 1. / (safmin * max(d__3,d__4));
		    scale = min(d__1,d__2);
		    if (lsa) {
			acoef = ascale * (scale * sbeta);
		    } else {
			acoef = scale * acoef;
		    }
		    if (lsb) {
			bcoefr = bscale * (scale * salfar);
		    } else {
			bcoefr = scale * bcoefr;
		    }
		}
		acoefa = abs(acoef);
		bcoefa = abs(bcoefr);

/*              First component is 1 */

		work[(*n << 1) + je] = 1.;
		xmax = 1.;
	    } else {

/*              Complex eigenvalue */

		d__1 = safmin * 100.;
		dlag2_(&a_ref(je, je), lda, &b_ref(je, je), ldb, &d__1, &
			acoef, &temp, &bcoefr, &temp2, &bcoefi);
		bcoefi = -bcoefi;
		if (bcoefi == 0.) {
		    *info = je;
		    return 0;
		}

/*              Scale to avoid over/underflow */

		acoefa = abs(acoef);
		bcoefa = abs(bcoefr) + abs(bcoefi);
		scale = 1.;
		if (acoefa * ulp < safmin && acoefa >= safmin) {
		    scale = safmin / ulp / acoefa;
		}
		if (bcoefa * ulp < safmin && bcoefa >= safmin) {
/* Computing MAX */
		    d__1 = scale, d__2 = safmin / ulp / bcoefa;
		    scale = max(d__1,d__2);
		}
		if (safmin * acoefa > ascale) {
		    scale = ascale / (safmin * acoefa);
		}
		if (safmin * bcoefa > bscale) {
/* Computing MIN */
		    d__1 = scale, d__2 = bscale / (safmin * bcoefa);
		    scale = min(d__1,d__2);
		}
		if (scale != 1.) {
		    acoef = scale * acoef;
		    acoefa = abs(acoef);
		    bcoefr = scale * bcoefr;
		    bcoefi = scale * bcoefi;
		    bcoefa = abs(bcoefr) + abs(bcoefi);
		}

/*              Compute first two components of eigenvector */

		temp = acoef * a_ref(je + 1, je);
		temp2r = acoef * a_ref(je, je) - bcoefr * b_ref(je, je);
		temp2i = -bcoefi * b_ref(je, je);
		if (abs(temp) > abs(temp2r) + abs(temp2i)) {
		    work[(*n << 1) + je] = 1.;
		    work[*n * 3 + je] = 0.;
		    work[(*n << 1) + je + 1] = -temp2r / temp;
		    work[*n * 3 + je + 1] = -temp2i / temp;
		} else {
		    work[(*n << 1) + je + 1] = 1.;
		    work[*n * 3 + je + 1] = 0.;
		    temp = acoef * a_ref(je, je + 1);
		    work[(*n << 1) + je] = (bcoefr * b_ref(je + 1, je + 1) - 
			    acoef * a_ref(je + 1, je + 1)) / temp;
		    work[*n * 3 + je] = bcoefi * b_ref(je + 1, je + 1) / temp;
		}
/* Computing MAX */
		d__5 = (d__1 = work[(*n << 1) + je], abs(d__1)) + (d__2 = 
			work[*n * 3 + je], abs(d__2)), d__6 = (d__3 = work[(*
			n << 1) + je + 1], abs(d__3)) + (d__4 = work[*n * 3 + 
			je + 1], abs(d__4));
		xmax = max(d__5,d__6);
	    }

/* Computing MAX */
	    d__1 = ulp * acoefa * anorm, d__2 = ulp * bcoefa * bnorm, d__1 = 
		    max(d__1,d__2);
	    dmin__ = max(d__1,safmin);

/*                                           T   
             Triangular solve of  (a A - b B)  y = 0   

                                     T   
             (rowwise in  (a A - b B) , or columnwise in (a A - b B) ) */

	    il2by2 = FALSE_;
/*           ------------------- Begin Timing Code ---------------------- */
	    opst = 0.;
	    in2by2 = 0;
/*           -------------------- End Timing Code ----------------------- */

	    i__2 = *n;
	    for (j = je + nw; j <= i__2; ++j) {
/*              ------------------- Begin Timing Code ------------------- */
		opssca = (doublereal) (nw * (j - je) + 1);
/*              -------------------- End Timing Code -------------------- */
		if (il2by2) {
		    il2by2 = FALSE_;
		    goto L160;
		}

		na = 1;
		bdiag[0] = b_ref(j, j);
		if (j < *n) {
		    if (a_ref(j + 1, j) != 0.) {
			il2by2 = TRUE_;
			bdiag[1] = b_ref(j + 1, j + 1);
			na = 2;
/*                    ---------------- Begin Timing Code ---------------- */
			++in2by2;
/*                    ----------------- End Timing Code ----------------- */
		    }
		}

/*              Check whether scaling is necessary for dot products */

		xscale = 1. / max(1.,xmax);
/* Computing MAX */
		d__1 = work[j], d__2 = work[*n + j], d__1 = max(d__1,d__2), 
			d__2 = acoefa * work[j] + bcoefa * work[*n + j];
		temp = max(d__1,d__2);
		if (il2by2) {
/* Computing MAX */
		    d__1 = temp, d__2 = work[j + 1], d__1 = max(d__1,d__2), 
			    d__2 = work[*n + j + 1], d__1 = max(d__1,d__2), 
			    d__2 = acoefa * work[j + 1] + bcoefa * work[*n + 
			    j + 1];
		    temp = max(d__1,d__2);
		}
		if (temp > bignum * xscale) {
		    i__3 = nw - 1;
		    for (jw = 0; jw <= i__3; ++jw) {
			i__4 = j - 1;
			for (jr = je; jr <= i__4; ++jr) {
			    work[(jw + 2) * *n + jr] = xscale * work[(jw + 2) 
				    * *n + jr];