strevc.c

提供矩阵类的函数库

		    }
		    j1 = j;
		    j2 = j;
		    jnxt = j + 1;
		    if (j < *n) {
			if (t_ref(j + 1, j) != 0.f) {
			    j2 = j + 1;
			    jnxt = j + 2;
			}
		    }

		    if (j1 == j2) {

/*                    1-by-1 diagonal block   

                      Scale if necessary to avoid overflow when forming   
                      the right-hand side. */

			if (work[j] > vcrit) {
			    rec = 1.f / vmax;
			    i__3 = *n - ki + 1;
			    sscal_(&i__3, &rec, &work[ki + *n], &c__1);
			    vmax = 1.f;
			    vcrit = bignum;
			}

			i__3 = j - ki - 1;
			work[j + *n] -= sdot_(&i__3, &t_ref(ki + 1, j), &c__1,
				 &work[ki + 1 + *n], &c__1);

/*                    Solve (T(J,J)-WR)'*X = WORK */

			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 if necessary */

			if (scale != 1.f) {
			    i__3 = *n - ki + 1;
			    sscal_(&i__3, &scale, &work[ki + *n], &c__1);
			}
			work[j + *n] = x_ref(1, 1);
/* Computing MAX */
			r__2 = (r__1 = work[j + *n], dabs(r__1));
			vmax = dmax(r__2,vmax);
			vcrit = bignum / vmax;
/* **   
                      Increment op count, ignoring the possible scaling */
			opst += (j - ki - 1 << 1) + 6;
/* ** */

		    } else {

/*                    2-by-2 diagonal block   

                      Scale if necessary to avoid overflow when forming   
                      the right-hand side.   

   Computing MAX */
			r__1 = work[j], r__2 = work[j + 1];
			beta = dmax(r__1,r__2);
			if (beta > vcrit) {
			    rec = 1.f / vmax;
			    i__3 = *n - ki + 1;
			    sscal_(&i__3, &rec, &work[ki + *n], &c__1);
			    vmax = 1.f;
			    vcrit = bignum;
			}

			i__3 = j - ki - 1;
			work[j + *n] -= sdot_(&i__3, &t_ref(ki + 1, j), &c__1,
				 &work[ki + 1 + *n], &c__1);

			i__3 = j - ki - 1;
			work[j + 1 + *n] -= sdot_(&i__3, &t_ref(ki + 1, j + 1)
				, &c__1, &work[ki + 1 + *n], &c__1);

/*                    Solve   
                        [T(J,J)-WR   T(J,J+1)     ]'* X = SCALE*( WORK1 )   
                        [T(J+1,J)    T(J+1,J+1)-WR]             ( WORK2 ) */

			slaln2_(&c_true, &c__2, &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 if necessary */

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

/* Computing MAX */
			r__3 = (r__1 = work[j + *n], dabs(r__1)), r__4 = (
				r__2 = work[j + 1 + *n], dabs(r__2)), r__3 = 
				max(r__3,r__4);
			vmax = dmax(r__3,vmax);
			vcrit = bignum / vmax;
/* **   
                      Increment op count, ignoring the possible scaling */
			opst += (j - ki - 1 << 2) + 24;
/* ** */

		    }
L170:
		    ;
		}

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

		if (! over) {
		    i__2 = *n - ki + 1;
		    scopy_(&i__2, &work[ki + *n], &c__1, &vl_ref(ki, is), &
			    c__1);

		    i__2 = *n - ki + 1;
		    ii = isamax_(&i__2, &vl_ref(ki, is), &c__1) + ki - 1;
		    remax = 1.f / (r__1 = vl_ref(ii, is), dabs(r__1));
		    i__2 = *n - ki + 1;
		    sscal_(&i__2, &remax, &vl_ref(ki, is), &c__1);
/* ** */
		    opst += (*n - ki + 1 << 1) + 1;
/* ** */

		    i__2 = ki - 1;
		    for (k = 1; k <= i__2; ++k) {
			vl_ref(k, is) = 0.f;
/* L180: */
		    }

		} else {

		    if (ki < *n) {
			i__2 = *n - ki;
			sgemv_("N", n, &i__2, &c_b22, &vl_ref(1, ki + 1), 
				ldvl, &work[ki + 1 + *n], &c__1, &work[ki + *
				n], &vl_ref(1, ki), &c__1);
		    }

		    ii = isamax_(n, &vl_ref(1, ki), &c__1);
		    remax = 1.f / (r__1 = vl_ref(ii, ki), dabs(r__1));
		    sscal_(n, &remax, &vl_ref(1, ki), &c__1);
/* ** */
		    latime_1.ops += (*n << 1) * (*n - ki + 1) + 1;
/* ** */

		}

	    } else {

/*              Complex left eigenvector.   

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

		if ((r__1 = t_ref(ki, ki + 1), dabs(r__1)) >= (r__2 = t_ref(
			ki + 1, ki), dabs(r__2))) {
		    work[ki + *n] = wi / t_ref(ki, ki + 1);
		    work[ki + 1 + n2] = 1.f;
		} else {
		    work[ki + *n] = 1.f;
		    work[ki + 1 + n2] = -wi / t_ref(ki + 1, ki);
		}
		work[ki + 1 + *n] = 0.f;
		work[ki + n2] = 0.f;

/*              Form right-hand side */

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

                Solve complex quasi-triangular system:   
                ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2 */

		vmax = 1.f;
		vcrit = bignum;

		jnxt = ki + 2;
		i__2 = *n;
		for (j = ki + 2; j <= i__2; ++j) {
		    if (j < jnxt) {
			goto L200;
		    }
		    j1 = j;
		    j2 = j;
		    jnxt = j + 1;
		    if (j < *n) {
			if (t_ref(j + 1, j) != 0.f) {
			    j2 = j + 1;
			    jnxt = j + 2;
			}
		    }

		    if (j1 == j2) {

/*                    1-by-1 diagonal block   

                      Scale if necessary to avoid overflow when   
                      forming the right-hand side elements. */

			if (work[j] > vcrit) {
			    rec = 1.f / vmax;
			    i__3 = *n - ki + 1;
			    sscal_(&i__3, &rec, &work[ki + *n], &c__1);
			    i__3 = *n - ki + 1;
			    sscal_(&i__3, &rec, &work[ki + n2], &c__1);
			    vmax = 1.f;
			    vcrit = bignum;
			}

			i__3 = j - ki - 2;
			work[j + *n] -= sdot_(&i__3, &t_ref(ki + 2, j), &c__1,
				 &work[ki + 2 + *n], &c__1);
			i__3 = j - ki - 2;
			work[j + n2] -= sdot_(&i__3, &t_ref(ki + 2, j), &c__1,
				 &work[ki + 2 + n2], &c__1);

/*                    Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2 */

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

/*                    Scale if necessary */

			if (scale != 1.f) {
			    i__3 = *n - ki + 1;
			    sscal_(&i__3, &scale, &work[ki + *n], &c__1);
			    i__3 = *n - ki + 1;
			    sscal_(&i__3, &scale, &work[ki + n2], &c__1);
			}
			work[j + *n] = x_ref(1, 1);
			work[j + n2] = x_ref(1, 2);
/* Computing MAX */
			r__3 = (r__1 = work[j + *n], dabs(r__1)), r__4 = (
				r__2 = work[j + n2], dabs(r__2)), r__3 = max(
				r__3,r__4);
			vmax = dmax(r__3,vmax);
			vcrit = bignum / vmax;
/* **   
                      Increment op count, ignoring the possible scaling */
			opst += (j - ki - 2 << 2) + 24;
/* ** */

		    } else {

/*                    2-by-2 diagonal block   

                      Scale if necessary to avoid overflow when forming   
                      the right-hand side elements.   

   Computing MAX */
			r__1 = work[j], r__2 = work[j + 1];
			beta = dmax(r__1,r__2);
			if (beta > vcrit) {
			    rec = 1.f / vmax;
			    i__3 = *n - ki + 1;
			    sscal_(&i__3, &rec, &work[ki + *n], &c__1);
			    i__3 = *n - ki + 1;
			    sscal_(&i__3, &rec, &work[ki + n2], &c__1);
			    vmax = 1.f;
			    vcrit = bignum;
			}

			i__3 = j - ki - 2;
			work[j + *n] -= sdot_(&i__3, &t_ref(ki + 2, j), &c__1,
				 &work[ki + 2 + *n], &c__1);

			i__3 = j - ki - 2;
			work[j + n2] -= sdot_(&i__3, &t_ref(ki + 2, j), &c__1,
				 &work[ki + 2 + n2], &c__1);

			i__3 = j - ki - 2;
			work[j + 1 + *n] -= sdot_(&i__3, &t_ref(ki + 2, j + 1)
				, &c__1, &work[ki + 2 + *n], &c__1);

			i__3 = j - ki - 2;
			work[j + 1 + n2] -= sdot_(&i__3, &t_ref(ki + 2, j + 1)
				, &c__1, &work[ki + 2 + n2], &c__1);

/*                    Solve 2-by-2 complex linear equation   
                        ([T(j,j)   T(j,j+1)  ]'-(wr-i*wi)*I)*X = SCALE*B   
                        ([T(j+1,j) T(j+1,j+1)]             ) */

			r__1 = -wi;
			slaln2_(&c_true, &c__2, &c__2, &smin, &c_b22, &t_ref(
				j, j), ldt, &c_b22, &c_b22, &work[j + *n], n, 
				&wr, &r__1, x, &c__2, &scale, &xnorm, &ierr);

/*                    Scale if necessary */

			if (scale != 1.f) {
			    i__3 = *n - ki + 1;
			    sscal_(&i__3, &scale, &work[ki + *n], &c__1);
			    i__3 = *n - ki + 1;
			    sscal_(&i__3, &scale, &work[ki + n2], &c__1);
			}
			work[j + *n] = x_ref(1, 1);
			work[j + n2] = x_ref(1, 2);
			work[j + 1 + *n] = x_ref(2, 1);
			work[j + 1 + n2] = x_ref(2, 2);
/* Computing MAX */
			r__5 = (r__1 = x_ref(1, 1), dabs(r__1)), r__6 = (r__2 
				= x_ref(1, 2), dabs(r__2)), r__5 = max(r__5,
				r__6), r__6 = (r__3 = x_ref(2, 1), dabs(r__3))
				, r__5 = max(r__5,r__6), r__6 = (r__4 = x_ref(
				2, 2), dabs(r__4)), r__5 = max(r__5,r__6);
			vmax = dmax(r__5,vmax);
			vcrit = bignum / vmax;
/* **   
                      Increment op count, ignoring the possible scaling */
			opst += (j - ki - 2 << 3) + 64;
/* ** */

		    }
L200:
		    ;
		}

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

   L210: */
		if (! over) {
		    i__2 = *n - ki + 1;
		    scopy_(&i__2, &work[ki + *n], &c__1, &vl_ref(ki, is), &
			    c__1);
		    i__2 = *n - ki + 1;
		    scopy_(&i__2, &work[ki + n2], &c__1, &vl_ref(ki, is + 1), 
			    &c__1);

		    emax = 0.f;
		    i__2 = *n;
		    for (k = ki; k <= i__2; ++k) {
/* Computing MAX */
			r__3 = emax, r__4 = (r__1 = vl_ref(k, is), dabs(r__1))
				 + (r__2 = vl_ref(k, is + 1), dabs(r__2));
			emax = dmax(r__3,r__4);
/* L220: */
		    }
		    remax = 1.f / emax;
		    i__2 = *n - ki + 1;
		    sscal_(&i__2, &remax, &vl_ref(ki, is), &c__1);
		    i__2 = *n - ki + 1;
		    sscal_(&i__2, &remax, &vl_ref(ki, is + 1), &c__1);
/* ** */
		    opst += (*n - ki + 1 << 2) + 1;
/* ** */

		    i__2 = ki - 1;
		    for (k = 1; k <= i__2; ++k) {
			vl_ref(k, is) = 0.f;
			vl_ref(k, is + 1) = 0.f;
/* L230: */
		    }
		} else {
		    if (ki < *n - 1) {
			i__2 = *n - ki - 1;
			sgemv_("N", n, &i__2, &c_b22, &vl_ref(1, ki + 2), 
				ldvl, &work[ki + 2 + *n], &c__1, &work[ki + *
				n], &vl_ref(1, ki), &c__1);
			i__2 = *n - ki - 1;
			sgemv_("N", n, &i__2, &c_b22, &vl_ref(1, ki + 2), 
				ldvl, &work[ki + 2 + n2], &c__1, &work[ki + 1 
				+ n2], &vl_ref(1, ki + 1), &c__1);
		    } else {
			sscal_(n, &work[ki + *n], &vl_ref(1, ki), &c__1);
			sscal_(n, &work[ki + 1 + n2], &vl_ref(1, ki + 1), &
				c__1);
		    }

		    emax = 0.f;
		    i__2 = *n;
		    for (k = 1; k <= i__2; ++k) {
/* Computing MAX */
			r__3 = emax, r__4 = (r__1 = vl_ref(k, ki), dabs(r__1))
				 + (r__2 = vl_ref(k, ki + 1), dabs(r__2));
			emax = dmax(r__3,r__4);
/* L240: */
		    }
		    remax = 1.f / emax;
		    sscal_(n, &remax, &vl_ref(1, ki), &c__1);
		    sscal_(n, &remax, &vl_ref(1, ki + 1), &c__1);
/* ** */
		    latime_1.ops += (*n << 2) * (*n - ki - 1) + *n * 6 + 1;
/* ** */

		}

	    }

	    ++is;
	    if (ip != 0) {
		++is;
	    }
L250:
	    if (ip == -1) {
		ip = 0;
	    }
	    if (ip == 1) {
		ip = -1;
	    }

/* L260: */
	}

    }
/* **   
       Compute final op count */
    latime_1.ops += opst;
/* ** */

    return 0;

/*     End of STREVC */

} /* strevc_ */

#undef vr_ref
#undef vl_ref
#undef x_ref
#undef t_ref