strevc.c

提供矩阵类的函数库

		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.f) {
			    j1 = j - 1;
			    jnxt = j - 2;
			}
		    }

		    if (j1 == j2) {

/*                    1-by-1 diagonal block */

			slaln2_(&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.f) {
			    if (work[j] > bignum / xnorm) {
				x_ref(1, 1) = x_ref(1, 1) / xnorm;
				scale /= xnorm;
			    }
			}

/*                    Scale if necessary */

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

/*                    Update right-hand side */

			i__1 = j - 1;
			r__1 = -x_ref(1, 1);
			saxpy_(&i__1, &r__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 */

			slaln2_(&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.f) {
/* Computing MAX */
			    r__1 = work[j - 1], r__2 = work[j];
			    beta = dmax(r__1,r__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.f) {
			    sscal_(&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;
			r__1 = -x_ref(1, 1);
			saxpy_(&i__1, &r__1, &t_ref(1, j - 1), &c__1, &work[*
				n + 1], &c__1);
			i__1 = j - 2;
			r__1 = -x_ref(2, 1);
			saxpy_(&i__1, &r__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) {
		    scopy_(&ki, &work[*n + 1], &c__1, &vr_ref(1, is), &c__1);

		    ii = isamax_(&ki, &vr_ref(1, is), &c__1);
		    remax = 1.f / (r__1 = vr_ref(ii, is), dabs(r__1));
		    sscal_(&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.f;
/* L70: */
		    }
		} else {
		    if (ki > 1) {
			i__1 = ki - 1;
			sgemv_("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 = isamax_(n, &vr_ref(1, ki), &c__1);
		    remax = 1.f / (r__1 = vr_ref(ii, ki), dabs(r__1));
		    sscal_(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 ((r__1 = t_ref(ki - 1, ki), dabs(r__1)) >= (r__2 = t_ref(
			ki, ki - 1), dabs(r__2))) {
		    work[ki - 1 + *n] = 1.f;
		    work[ki + n2] = wi / t_ref(ki - 1, ki);
		} else {
		    work[ki - 1 + *n] = -wi / t_ref(ki, ki - 1);
		    work[ki + n2] = 1.f;
		}
		work[ki + *n] = 0.f;
		work[ki - 1 + n2] = 0.f;

/*              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.f) {
			    j1 = j - 1;
			    jnxt = j - 2;
			}
		    }

		    if (j1 == j2) {

/*                    1-by-1 diagonal block */

			slaln2_(&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.f) {
			    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.f) {
			    sscal_(&ki, &scale, &work[*n + 1], &c__1);
			    sscal_(&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;
			r__1 = -x_ref(1, 1);
			saxpy_(&i__1, &r__1, &t_ref(1, j), &c__1, &work[*n + 
				1], &c__1);
			i__1 = j - 1;
			r__1 = -x_ref(1, 2);
			saxpy_(&i__1, &r__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 */

			slaln2_(&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.f) {
/* Computing MAX */
			    r__1 = work[j - 1], r__2 = work[j];
			    beta = dmax(r__1,r__2);
			    if (beta > bignum / xnorm) {
				rec = 1.f / 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.f) {
			    sscal_(&ki, &scale, &work[*n + 1], &c__1);
			    sscal_(&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;
			r__1 = -x_ref(1, 1);
			saxpy_(&i__1, &r__1, &t_ref(1, j - 1), &c__1, &work[*
				n + 1], &c__1);
			i__1 = j - 2;
			r__1 = -x_ref(2, 1);
			saxpy_(&i__1, &r__1, &t_ref(1, j), &c__1, &work[*n + 
				1], &c__1);
			i__1 = j - 2;
			r__1 = -x_ref(1, 2);
			saxpy_(&i__1, &r__1, &t_ref(1, j - 1), &c__1, &work[
				n2 + 1], &c__1);
			i__1 = j - 2;
			r__1 = -x_ref(2, 2);
			saxpy_(&i__1, &r__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) {
		    scopy_(&ki, &work[*n + 1], &c__1, &vr_ref(1, is - 1), &
			    c__1);
		    scopy_(&ki, &work[n2 + 1], &c__1, &vr_ref(1, is), &c__1);

		    emax = 0.f;
		    i__1 = ki;
		    for (k = 1; k <= i__1; ++k) {
/* Computing MAX */
			r__3 = emax, r__4 = (r__1 = vr_ref(k, is - 1), dabs(
				r__1)) + (r__2 = vr_ref(k, is), dabs(r__2));
			emax = dmax(r__3,r__4);
/* L100: */
		    }

		    remax = 1.f / emax;
		    sscal_(&ki, &remax, &vr_ref(1, is - 1), &c__1);
		    sscal_(&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.f;
			vr_ref(k, is) = 0.f;
/* L110: */
		    }

		} else {

		    if (ki > 2) {
			i__1 = ki - 2;
			sgemv_("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;
			sgemv_("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 {
			sscal_(n, &work[ki - 1 + *n], &vr_ref(1, ki - 1), &
				c__1);
			sscal_(n, &work[ki + n2], &vr_ref(1, ki), &c__1);
		    }

		    emax = 0.f;
		    i__1 = *n;
		    for (k = 1; k <= i__1; ++k) {
/* Computing MAX */
			r__3 = emax, r__4 = (r__1 = vr_ref(k, ki - 1), dabs(
				r__1)) + (r__2 = vr_ref(k, ki), dabs(r__2));
			emax = dmax(r__3,r__4);
/* L120: */
		    }
		    remax = 1.f / emax;
		    sscal_(n, &remax, &vr_ref(1, ki - 1), &c__1);
		    sscal_(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.f) {
		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.f;
	    if (ip != 0) {
		wi = sqrt((r__1 = t_ref(ki, ki + 1), dabs(r__1))) * sqrt((
			r__2 = t_ref(ki + 1, ki), dabs(r__2)));
	    }
/* Computing MAX */
	    r__1 = ulp * (dabs(wr) + dabs(wi));
	    smin = dmax(r__1,smlnum);

	    if (ip == 0) {

/*              Real left eigenvector. */

		work[ki + *n] = 1.f;

/*              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.f;
		vcrit = bignum;

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