ctrevc.c

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

    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	i__2 = i__ + *n;
	i__3 = t_subscr(i__, i__);
	work[i__2].r = t[i__3].r, work[i__2].i = t[i__3].i;
/* L20: */
    }

/*     Compute 1-norm of each column of strictly upper triangular   
       part of T to control overflow in triangular solver. */

    rwork[1] = 0.f;
    i__1 = *n;
    for (j = 2; j <= i__1; ++j) {
	i__2 = j - 1;
	rwork[j] = scasum_(&i__2, &t_ref(1, j), &c__1);
/* L30: */
    }
/* ** */
    latime_1.ops += *n * (*n - 1);
/* ** */

    if (rightv) {

/*        Compute right eigenvectors. */

	is = *m;
	for (ki = *n; ki >= 1; --ki) {

	    if (somev) {
		if (! select[ki]) {
		    goto L80;
		}
	    }
/* Computing MAX */
	    i__1 = t_subscr(ki, ki);
	    r__3 = ulp * ((r__1 = t[i__1].r, dabs(r__1)) + (r__2 = r_imag(&
		    t_ref(ki, ki)), dabs(r__2)));
	    smin = dmax(r__3,smlnum);

	    work[1].r = 1.f, work[1].i = 0.f;

/*           Form right-hand side. */

	    i__1 = ki - 1;
	    for (k = 1; k <= i__1; ++k) {
		i__2 = k;
		i__3 = t_subscr(k, ki);
		q__1.r = -t[i__3].r, q__1.i = -t[i__3].i;
		work[i__2].r = q__1.r, work[i__2].i = q__1.i;
/* L40: */
	    }

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

	    i__1 = ki - 1;
	    for (k = 1; k <= i__1; ++k) {
		i__2 = t_subscr(k, k);
		i__3 = t_subscr(k, k);
		i__4 = t_subscr(ki, ki);
		q__1.r = t[i__3].r - t[i__4].r, q__1.i = t[i__3].i - t[i__4]
			.i;
		t[i__2].r = q__1.r, t[i__2].i = q__1.i;
		i__2 = t_subscr(k, k);
		if ((r__1 = t[i__2].r, dabs(r__1)) + (r__2 = r_imag(&t_ref(k, 
			k)), dabs(r__2)) < smin) {
		    i__3 = t_subscr(k, k);
		    t[i__3].r = smin, t[i__3].i = 0.f;
		}
/* L50: */
	    }
/* ** */
	    opst += ki - 1 << 1;
/* ** */

	    if (ki > 1) {
		i__1 = ki - 1;
		clatrs_("Upper", "No transpose", "Non-unit", "Y", &i__1, &t[
			t_offset], ldt, &work[1], &scale, &rwork[1], info);
		i__1 = ki;
		work[i__1].r = scale, work[i__1].i = 0.f;
	    }
/* **   
             Increment opcount for triangular solver, assuming that   
             ops CLATRS = ops CTRSV, with no scaling in CLATRS. */
	    latime_1.ops += (ki << 2) * (ki - 1);
/* **   

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

	    if (! over) {
		ccopy_(&ki, &work[1], &c__1, &vr_ref(1, is), &c__1);

		ii = icamax_(&ki, &vr_ref(1, is), &c__1);
		i__1 = vr_subscr(ii, is);
		remax = 1.f / ((r__1 = vr[i__1].r, dabs(r__1)) + (r__2 = 
			r_imag(&vr_ref(ii, is)), dabs(r__2)));
		csscal_(&ki, &remax, &vr_ref(1, is), &c__1);
/* ** */
		opst += (ki << 2) + 3;
/* ** */

		i__1 = *n;
		for (k = ki + 1; k <= i__1; ++k) {
		    i__2 = vr_subscr(k, is);
		    vr[i__2].r = 0.f, vr[i__2].i = 0.f;
/* L60: */
		}
	    } else {
		if (ki > 1) {
		    i__1 = ki - 1;
		    q__1.r = scale, q__1.i = 0.f;
		    cgemv_("N", n, &i__1, &c_b2, &vr[vr_offset], ldvr, &work[
			    1], &c__1, &q__1, &vr_ref(1, ki), &c__1);
		}

		ii = icamax_(n, &vr_ref(1, ki), &c__1);
		i__1 = vr_subscr(ii, ki);
		remax = 1.f / ((r__1 = vr[i__1].r, dabs(r__1)) + (r__2 = 
			r_imag(&vr_ref(ii, ki)), dabs(r__2)));
		csscal_(n, &remax, &vr_ref(1, ki), &c__1);
/* ** */
		latime_1.ops += (*n << 3) * (ki - 1) + *n * 10 + 3;
/* ** */
	    }

/*           Set back the original diagonal elements of T. */

	    i__1 = ki - 1;
	    for (k = 1; k <= i__1; ++k) {
		i__2 = t_subscr(k, k);
		i__3 = k + *n;
		t[i__2].r = work[i__3].r, t[i__2].i = work[i__3].i;
/* L70: */
	    }

	    --is;
L80:
	    ;
	}
    }

    if (leftv) {

/*        Compute left eigenvectors. */

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

	    if (somev) {
		if (! select[ki]) {
		    goto L130;
		}
	    }
/* Computing MAX */
	    i__2 = t_subscr(ki, ki);
	    r__3 = ulp * ((r__1 = t[i__2].r, dabs(r__1)) + (r__2 = r_imag(&
		    t_ref(ki, ki)), dabs(r__2)));
	    smin = dmax(r__3,smlnum);

	    i__2 = *n;
	    work[i__2].r = 1.f, work[i__2].i = 0.f;

/*           Form right-hand side. */

	    i__2 = *n;
	    for (k = ki + 1; k <= i__2; ++k) {
		i__3 = k;
		r_cnjg(&q__2, &t_ref(ki, k));
		q__1.r = -q__2.r, q__1.i = -q__2.i;
		work[i__3].r = q__1.r, work[i__3].i = q__1.i;
/* L90: */
	    }

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

	    i__2 = *n;
	    for (k = ki + 1; k <= i__2; ++k) {
		i__3 = t_subscr(k, k);
		i__4 = t_subscr(k, k);
		i__5 = t_subscr(ki, ki);
		q__1.r = t[i__4].r - t[i__5].r, q__1.i = t[i__4].i - t[i__5]
			.i;
		t[i__3].r = q__1.r, t[i__3].i = q__1.i;
		i__3 = t_subscr(k, k);
		if ((r__1 = t[i__3].r, dabs(r__1)) + (r__2 = r_imag(&t_ref(k, 
			k)), dabs(r__2)) < smin) {
		    i__4 = t_subscr(k, k);
		    t[i__4].r = smin, t[i__4].i = 0.f;
		}
/* L100: */
	    }
/* ** */
	    opst += *n - ki << 1;
/* ** */

	    if (ki < *n) {
		i__2 = *n - ki;
		clatrs_("Upper", "Conjugate transpose", "Non-unit", "Y", &
			i__2, &t_ref(ki + 1, ki + 1), ldt, &work[ki + 1], &
			scale, &rwork[1], info);
		i__2 = ki;
		work[i__2].r = scale, work[i__2].i = 0.f;
	    }
/* **   
             Increment opcount for triangular solver, assuming that   
             ops CLATRS = ops CTRSV, with no scaling in CLATRS. */
	    latime_1.ops += (*n - ki << 2) * (*n - ki + 1);
/* **   

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

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

		i__2 = *n - ki + 1;
		ii = icamax_(&i__2, &vl_ref(ki, is), &c__1) + ki - 1;
		i__2 = vl_subscr(ii, is);
		remax = 1.f / ((r__1 = vl[i__2].r, dabs(r__1)) + (r__2 = 
			r_imag(&vl_ref(ii, is)), dabs(r__2)));
		i__2 = *n - ki + 1;
		csscal_(&i__2, &remax, &vl_ref(ki, is), &c__1);
/* ** */
		opst += (*n - ki + 1 << 2) + 3;
/* ** */

		i__2 = ki - 1;
		for (k = 1; k <= i__2; ++k) {
		    i__3 = vl_subscr(k, is);
		    vl[i__3].r = 0.f, vl[i__3].i = 0.f;
/* L110: */
		}
	    } else {
		if (ki < *n) {
		    i__2 = *n - ki;
		    q__1.r = scale, q__1.i = 0.f;
		    cgemv_("N", n, &i__2, &c_b2, &vl_ref(1, ki + 1), ldvl, &
			    work[ki + 1], &c__1, &q__1, &vl_ref(1, ki), &c__1);
		}

		ii = icamax_(n, &vl_ref(1, ki), &c__1);
		i__2 = vl_subscr(ii, ki);
		remax = 1.f / ((r__1 = vl[i__2].r, dabs(r__1)) + (r__2 = 
			r_imag(&vl_ref(ii, ki)), dabs(r__2)));
		csscal_(n, &remax, &vl_ref(1, ki), &c__1);
/* ** */
		latime_1.ops += (*n << 3) * (*n - ki) + *n * 10 + 3;
/* ** */
	    }

/*           Set back the original diagonal elements of T. */

	    i__2 = *n;
	    for (k = ki + 1; k <= i__2; ++k) {
		i__3 = t_subscr(k, k);
		i__4 = k + *n;
		t[i__3].r = work[i__4].r, t[i__3].i = work[i__4].i;
/* L120: */
	    }

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

    return 0;

/*     End of CTREVC */

} /* ctrevc_ */

#undef vr_ref
#undef vr_subscr
#undef vl_ref
#undef vl_subscr
#undef t_ref
#undef t_subscr