dtrevc.c

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

		}

/*              Solve the upper quasi-triangular system:   
                   (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK. */

		jnxt = ki - 1;
		for (j = ki - 1; j >= 1; --j) {
		    if (j > jnxt) {
			goto L60;
		    }
		    j1 = j;
		    j2 = j;
		    jnxt = j - 1;
		    if (j > 1) {
			if (t_ref(j, j - 1) != 0.) {
			    j1 = j - 1;
			    jnxt = j - 2;
			}
		    }

		    if (j1 == j2) {

/*                    1-by-1 diagonal block */

			dlaln2_(&c_false, &c__1, &c__1, &smin, &c_b22, &t_ref(
				j, j), ldt, &c_b22, &c_b22, &work[j + *n], n, 
				&wr, &c_b25, x, &c__2, &scale, &xnorm, &ierr);

/*                    Scale X(1,1) to avoid overflow when updating   
                      the right-hand side. */

			if (xnorm > 1.) {
			    if (work[j] > bignum / xnorm) {
				x_ref(1, 1) = x_ref(1, 1) / xnorm;
				scale /= xnorm;
			    }
			}

/*                    Scale if necessary */

			if (scale != 1.) {
			    dscal_(&ki, &scale, &work[*n + 1], &c__1);
			}
			work[j + *n] = x_ref(1, 1);

/*                    Update right-hand side */

			i__1 = j - 1;
			d__1 = -x_ref(1, 1);
			daxpy_(&i__1, &d__1, &t_ref(1, j), &c__1, &work[*n + 
				1], &c__1);
/* **   
                      Increment op count, ignoring the possible scaling */
			opst += (j - 1 << 1) + 6;
/* ** */

		    } else {

/*                    2-by-2 diagonal block */

			dlaln2_(&c_false, &c__2, &c__1, &smin, &c_b22, &t_ref(
				j - 1, j - 1), ldt, &c_b22, &c_b22, &work[j - 
				1 + *n], n, &wr, &c_b25, x, &c__2, &scale, &
				xnorm, &ierr);

/*                    Scale X(1,1) and X(2,1) to avoid overflow when   
                      updating the right-hand side. */

			if (xnorm > 1.) {
/* Computing MAX */
			    d__1 = work[j - 1], d__2 = work[j];
			    beta = max(d__1,d__2);
			    if (beta > bignum / xnorm) {
				x_ref(1, 1) = x_ref(1, 1) / xnorm;
				x_ref(2, 1) = x_ref(2, 1) / xnorm;
				scale /= xnorm;
			    }
			}

/*                    Scale if necessary */

			if (scale != 1.) {
			    dscal_(&ki, &scale, &work[*n + 1], &c__1);
			}
			work[j - 1 + *n] = x_ref(1, 1);
			work[j + *n] = x_ref(2, 1);

/*                    Update right-hand side */

			i__1 = j - 2;
			d__1 = -x_ref(1, 1);
			daxpy_(&i__1, &d__1, &t_ref(1, j - 1), &c__1, &work[*
				n + 1], &c__1);
			i__1 = j - 2;
			d__1 = -x_ref(2, 1);
			daxpy_(&i__1, &d__1, &t_ref(1, j), &c__1, &work[*n + 
				1], &c__1);
/* **   
                      Increment op count, ignoring the possible scaling */
			opst += (j - 2 << 2) + 24;
/* ** */
		    }
L60:
		    ;
		}

/*              Copy the vector x or Q*x to VR and normalize. */

		if (! over) {
		    dcopy_(&ki, &work[*n + 1], &c__1, &vr_ref(1, is), &c__1);

		    ii = idamax_(&ki, &vr_ref(1, is), &c__1);
		    remax = 1. / (d__1 = vr_ref(ii, is), abs(d__1));
		    dscal_(&ki, &remax, &vr_ref(1, is), &c__1);
/* ** */
		    opst += (ki << 1) + 1;
/* ** */

		    i__1 = *n;
		    for (k = ki + 1; k <= i__1; ++k) {
			vr_ref(k, is) = 0.;
/* L70: */
		    }
		} else {
		    if (ki > 1) {
			i__1 = ki - 1;
			dgemv_("N", n, &i__1, &c_b22, &vr[vr_offset], ldvr, &
				work[*n + 1], &c__1, &work[ki + *n], &vr_ref(
				1, ki), &c__1);
		    }

		    ii = idamax_(n, &vr_ref(1, ki), &c__1);
		    remax = 1. / (d__1 = vr_ref(ii, ki), abs(d__1));
		    dscal_(n, &remax, &vr_ref(1, ki), &c__1);
/* ** */
		    latime_1.ops += (*n << 1) * ki + 1;
/* ** */
		}

	    } else {

/*              Complex right eigenvector.   

                Initial solve   
                  [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.   
                  [ (T(KI,KI-1)   T(KI,KI)   )               ] */

		if ((d__1 = t_ref(ki - 1, ki), abs(d__1)) >= (d__2 = t_ref(ki,
			 ki - 1), abs(d__2))) {
		    work[ki - 1 + *n] = 1.;
		    work[ki + n2] = wi / t_ref(ki - 1, ki);
		} else {
		    work[ki - 1 + *n] = -wi / t_ref(ki, ki - 1);
		    work[ki + n2] = 1.;
		}
		work[ki + *n] = 0.;
		work[ki - 1 + n2] = 0.;

/*              Form right-hand side */

		i__1 = ki - 2;
		for (k = 1; k <= i__1; ++k) {
		    work[k + *n] = -work[ki - 1 + *n] * t_ref(k, ki - 1);
		    work[k + n2] = -work[ki + n2] * t_ref(k, ki);
/* L80: */
		}
/* ** */
		opst += ki - 2 << 1;
/* **   

                Solve upper quasi-triangular system:   
                (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2) */

		jnxt = ki - 2;
		for (j = ki - 2; j >= 1; --j) {
		    if (j > jnxt) {
			goto L90;
		    }
		    j1 = j;
		    j2 = j;
		    jnxt = j - 1;
		    if (j > 1) {
			if (t_ref(j, j - 1) != 0.) {
			    j1 = j - 1;
			    jnxt = j - 2;
			}
		    }

		    if (j1 == j2) {

/*                    1-by-1 diagonal block */

			dlaln2_(&c_false, &c__1, &c__2, &smin, &c_b22, &t_ref(
				j, j), ldt, &c_b22, &c_b22, &work[j + *n], n, 
				&wr, &wi, x, &c__2, &scale, &xnorm, &ierr);

/*                    Scale X(1,1) and X(1,2) to avoid overflow when   
                      updating the right-hand side. */

			if (xnorm > 1.) {
			    if (work[j] > bignum / xnorm) {
				x_ref(1, 1) = x_ref(1, 1) / xnorm;
				x_ref(1, 2) = x_ref(1, 2) / xnorm;
				scale /= xnorm;
			    }
			}

/*                    Scale if necessary */

			if (scale != 1.) {
			    dscal_(&ki, &scale, &work[*n + 1], &c__1);
			    dscal_(&ki, &scale, &work[n2 + 1], &c__1);
			}
			work[j + *n] = x_ref(1, 1);
			work[j + n2] = x_ref(1, 2);

/*                    Update the right-hand side */

			i__1 = j - 1;
			d__1 = -x_ref(1, 1);
			daxpy_(&i__1, &d__1, &t_ref(1, j), &c__1, &work[*n + 
				1], &c__1);
			i__1 = j - 1;
			d__1 = -x_ref(1, 2);
			daxpy_(&i__1, &d__1, &t_ref(1, j), &c__1, &work[n2 + 
				1], &c__1);
/* **   
                      Increment op count, ignoring the possible scaling */
			opst += (j - 1 << 2) + 24;
/* ** */

		    } else {

/*                    2-by-2 diagonal block */

			dlaln2_(&c_false, &c__2, &c__2, &smin, &c_b22, &t_ref(
				j - 1, j - 1), ldt, &c_b22, &c_b22, &work[j - 
				1 + *n], n, &wr, &wi, x, &c__2, &scale, &
				xnorm, &ierr);

/*                    Scale X to avoid overflow when updating   
                      the right-hand side. */

			if (xnorm > 1.) {
/* Computing MAX */
			    d__1 = work[j - 1], d__2 = work[j];
			    beta = max(d__1,d__2);
			    if (beta > bignum / xnorm) {
				rec = 1. / xnorm;
				x_ref(1, 1) = x_ref(1, 1) * rec;
				x_ref(1, 2) = x_ref(1, 2) * rec;
				x_ref(2, 1) = x_ref(2, 1) * rec;
				x_ref(2, 2) = x_ref(2, 2) * rec;
				scale *= rec;
			    }
			}

/*                    Scale if necessary */

			if (scale != 1.) {
			    dscal_(&ki, &scale, &work[*n + 1], &c__1);
			    dscal_(&ki, &scale, &work[n2 + 1], &c__1);
			}
			work[j - 1 + *n] = x_ref(1, 1);
			work[j + *n] = x_ref(2, 1);
			work[j - 1 + n2] = x_ref(1, 2);
			work[j + n2] = x_ref(2, 2);

/*                    Update the right-hand side */

			i__1 = j - 2;
			d__1 = -x_ref(1, 1);
			daxpy_(&i__1, &d__1, &t_ref(1, j - 1), &c__1, &work[*
				n + 1], &c__1);
			i__1 = j - 2;
			d__1 = -x_ref(2, 1);
			daxpy_(&i__1, &d__1, &t_ref(1, j), &c__1, &work[*n + 
				1], &c__1);
			i__1 = j - 2;
			d__1 = -x_ref(1, 2);
			daxpy_(&i__1, &d__1, &t_ref(1, j - 1), &c__1, &work[
				n2 + 1], &c__1);
			i__1 = j - 2;
			d__1 = -x_ref(2, 2);
			daxpy_(&i__1, &d__1, &t_ref(1, j), &c__1, &work[n2 + 
				1], &c__1);
/* **   
                      Increment op count, ignoring the possible scaling */
			opst += (j - 2 << 3) + 64;
/* ** */
		    }
L90:
		    ;
		}

/*              Copy the vector x or Q*x to VR and normalize. */

		if (! over) {
		    dcopy_(&ki, &work[*n + 1], &c__1, &vr_ref(1, is - 1), &
			    c__1);
		    dcopy_(&ki, &work[n2 + 1], &c__1, &vr_ref(1, is), &c__1);

		    emax = 0.;
		    i__1 = ki;
		    for (k = 1; k <= i__1; ++k) {
/* Computing MAX */
			d__3 = emax, d__4 = (d__1 = vr_ref(k, is - 1), abs(
				d__1)) + (d__2 = vr_ref(k, is), abs(d__2));
			emax = max(d__3,d__4);
/* L100: */
		    }

		    remax = 1. / emax;
		    dscal_(&ki, &remax, &vr_ref(1, is - 1), &c__1);
		    dscal_(&ki, &remax, &vr_ref(1, is), &c__1);
/* ** */
		    opst += (ki << 2) + 1;
/* ** */

		    i__1 = *n;
		    for (k = ki + 1; k <= i__1; ++k) {
			vr_ref(k, is - 1) = 0.;
			vr_ref(k, is) = 0.;
/* L110: */
		    }

		} else {

		    if (ki > 2) {
			i__1 = ki - 2;
			dgemv_("N", n, &i__1, &c_b22, &vr[vr_offset], ldvr, &
				work[*n + 1], &c__1, &work[ki - 1 + *n], &
				vr_ref(1, ki - 1), &c__1);
			i__1 = ki - 2;
			dgemv_("N", n, &i__1, &c_b22, &vr[vr_offset], ldvr, &
				work[n2 + 1], &c__1, &work[ki + n2], &vr_ref(
				1, ki), &c__1);
		    } else {
			dscal_(n, &work[ki - 1 + *n], &vr_ref(1, ki - 1), &
				c__1);
			dscal_(n, &work[ki + n2], &vr_ref(1, ki), &c__1);
		    }

		    emax = 0.;
		    i__1 = *n;
		    for (k = 1; k <= i__1; ++k) {
/* Computing MAX */
			d__3 = emax, d__4 = (d__1 = vr_ref(k, ki - 1), abs(
				d__1)) + (d__2 = vr_ref(k, ki), abs(d__2));
			emax = max(d__3,d__4);
/* L120: */
		    }
		    remax = 1. / emax;
		    dscal_(n, &remax, &vr_ref(1, ki - 1), &c__1);
		    dscal_(n, &remax, &vr_ref(1, ki), &c__1);
/* ** */
		    latime_1.ops += (*n << 2) * (ki - 2) + *n * 6 + 1;
/* ** */
		}
	    }

	    --is;
	    if (ip != 0) {
		--is;
	    }
L130:
	    if (ip == 1) {
		ip = 0;
	    }
	    if (ip == -1) {
		ip = 1;
	    }
/* L140: */
	}
    }

    if (leftv) {

/*        Compute left eigenvectors. */

	ip = 0;
	is = 1;
	i__1 = *n;
	for (ki = 1; ki <= i__1; ++ki) {

	    if (ip == -1) {
		goto L250;
	    }
	    if (ki == *n) {
		goto L150;
	    }
	    if (t_ref(ki + 1, ki) == 0.) {
		goto L150;
	    }
	    ip = 1;

L150:
	    if (somev) {
		if (! select[ki]) {
		    goto L250;
		}
	    }

/*           Compute the KI-th eigenvalue (WR,WI). */

	    wr = t_ref(ki, ki);
	    wi = 0.;
	    if (ip != 0) {
		wi = sqrt((d__1 = t_ref(ki, ki + 1), abs(d__1))) * sqrt((d__2 
			= t_ref(ki + 1, ki), abs(d__2)));
	    }
/* Computing MAX */
	    d__1 = ulp * (abs(wr) + abs(wi));
	    smin = max(d__1,smlnum);

	    if (ip == 0) {

/*              Real left eigenvector. */

		work[ki + *n] = 1.;

/*              Form right-hand side */

		i__2 = *n;
		for (k = ki + 1; k <= i__2; ++k) {
		    work[k + *n] = -t_ref(ki, k);
/* L160: */
		}

/*              Solve the quasi-triangular system:   
                   (T(KI+1:N,KI+1:N) - WR)'*X = SCALE*WORK */

		vmax = 1.;
		vcrit = bignum;

		jnxt = ki + 1;
		i__2 = *n;
		for (j = ki + 1; j <= i__2; ++j) {