cstegr.c

提供矩阵类的函数库

       Test the input parameters.   

       Parameter adjustments */
    --d__;
    --e;
    --w;
    z_dim1 = *ldz;
    z_offset = 1 + z_dim1 * 1;
    z__ -= z_offset;
    --isuppz;
    --work;
    --iwork;

    /* Function Body */
    wantz = lsame_(jobz, "V");
    alleig = lsame_(range, "A");
    valeig = lsame_(range, "V");
    indeig = lsame_(range, "I");

    lquery = *lwork == -1 || *liwork == -1;
    lwmin = *n * 18;
    liwmin = *n * 10;

    *info = 0;
    if (! (wantz || lsame_(jobz, "N"))) {
	*info = -1;
    } else if (! (alleig || valeig || indeig)) {
	*info = -2;

/*     The following two lines need to be removed once the   
       RANGE = 'V' and RANGE = 'I' options are provided. */

    } else if (valeig || indeig) {
	*info = -2;
    } else if (*n < 0) {
	*info = -3;
    } else if (valeig && *n > 0 && *vu <= *vl) {
	*info = -7;
    } else if (indeig && *il < 1) {
	*info = -8;
/*     The following change should be made in DSTEVX also, otherwise   
       IL can be specified as N+1 and IU as N.   
       ELSE IF( INDEIG .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) THEN */
    } else if (indeig && (*iu < *il || *iu > *n)) {
	*info = -9;
    } else if (*ldz < 1 || wantz && *ldz < *n) {
	*info = -14;
    } else if (*lwork < lwmin && ! lquery) {
	*info = -17;
    } else if (*liwork < liwmin && ! lquery) {
	*info = -19;
    }
    if (*info == 0) {
	work[1] = (real) lwmin;
	iwork[1] = liwmin;
    }

    if (*info != 0) {
	i__1 = -(*info);
	xerbla_("SSTEGR", &i__1);
	return 0;
    } else if (lquery) {
	return 0;
    }

/*     Quick return if possible */

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

    if (*n == 1) {
	if (alleig || indeig) {
	    *m = 1;
	    w[1] = d__[1];
	} else {
	    if (*vl < d__[1] && *vu >= d__[1]) {
		*m = 1;
		w[1] = d__[1];
	    }
	}
	if (wantz) {
	    i__1 = z___subscr(1, 1);
	    z__[i__1].r = 1.f, z__[i__1].i = 0.f;
	}
	return 0;
    }

/*     Get machine constants. */

    latime_1.ops += 7.f;
    safmin = slamch_("Safe minimum");
    eps = slamch_("Precision");
    smlnum = safmin / eps;
    bignum = 1.f / smlnum;
    rmin = sqrt(smlnum);
/* Computing MIN */
    r__1 = sqrt(bignum), r__2 = 1.f / sqrt(sqrt(safmin));
    rmax = dmin(r__1,r__2);

/*     Scale matrix to allowable range, if necessary. */

    scale = 1.f;
    tnrm = slanst_("M", n, &d__[1], &e[1]);
    if (tnrm > 0.f && tnrm < rmin) {
	latime_1.ops += 1.f;
	scale = rmin / tnrm;
    } else if (tnrm > rmax) {
	latime_1.ops += 1.f;
	scale = rmax / tnrm;
    }
    if (scale != 1.f) {
	latime_1.ops += (real) (*n << 1);
	sscal_(n, &scale, &d__[1], &c__1);
	i__1 = *n - 1;
	sscal_(&i__1, &scale, &e[1], &c__1);
	tnrm *= scale;
    }
    indgrs = 1;
    indwof = (*n << 1) + 1;
    indwrk = *n * 3 + 1;

    iinspl = 1;
    iindbl = *n + 1;
    iindwk = (*n << 1) + 1;

    claset_("Full", n, n, &c_b1, &c_b1, &z__[z_offset], ldz);

/*     Compute the desired eigenvalues of the tridiagonal after splitting   
       into smaller subblocks if the corresponding of-diagonal elements   
       are small */

    latime_1.ops += 1.f;
    thresh = eps * tnrm;
    slarre_(n, &d__[1], &e[1], &thresh, &nsplit, &iwork[iinspl], m, &w[1], &
	    work[indwof], &work[indgrs], &work[indwrk], &iinfo);
    if (iinfo != 0) {
	*info = 1;
	return 0;
    }

    if (wantz) {

/*        Compute the desired eigenvectors corresponding to the computed   
          eigenvalues */

	latime_1.ops += 1.f;
/* Computing MAX */
	r__1 = *abstol, r__2 = (real) (*n) * thresh;
	tol = dmax(r__1,r__2);
	ibegin = 1;
	i__1 = nsplit;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    iend = iwork[iinspl + i__ - 1];
	    i__2 = iend;
	    for (j = ibegin; j <= i__2; ++j) {
		iwork[iindbl + j - 1] = i__;
/* L10: */
	    }
	    ibegin = iend + 1;
/* L20: */
	}

	clarrv_(n, &d__[1], &e[1], &iwork[iinspl], m, &w[1], &iwork[iindbl], &
		work[indgrs], &tol, &z__[z_offset], ldz, &isuppz[1], &work[
		indwrk], &iwork[iindwk], &iinfo);
	if (iinfo != 0) {
	    *info = 2;
	    return 0;
	}

    }

    ibegin = 1;
    i__1 = nsplit;
    for (i__ = 1; i__ <= i__1; ++i__) {
	iend = iwork[iinspl + i__ - 1];
	i__2 = iend;
	for (j = ibegin; j <= i__2; ++j) {
	    latime_1.ops += 1.f;
	    w[j] += work[indwof + i__ - 1];
/* L30: */
	}
	ibegin = iend + 1;
/* L40: */
    }

/*     If matrix was scaled, then rescale eigenvalues appropriately. */

    if (scale != 1.f) {
	r__1 = 1.f / scale;
	sscal_(m, &r__1, &w[1], &c__1);
    }

/*     If eigenvalues are not in order, then sort them, along with   
       eigenvectors. */

    if (nsplit > 1) {
	i__1 = *m - 1;
	for (j = 1; j <= i__1; ++j) {
	    i__ = 0;
	    tmp = w[j];
	    i__2 = *m;
	    for (jj = j + 1; jj <= i__2; ++jj) {
		if (w[jj] < tmp) {
		    i__ = jj;
		    tmp = w[jj];
		}
/* L50: */
	    }
	    if (i__ != 0) {
		w[i__] = w[j];
		w[j] = tmp;
		if (wantz) {
		    cswap_(n, &z___ref(1, i__), &c__1, &z___ref(1, j), &c__1);
		    itmp = isuppz[(i__ << 1) - 1];
		    isuppz[(i__ << 1) - 1] = isuppz[(j << 1) - 1];
		    isuppz[(j << 1) - 1] = itmp;
		    itmp = isuppz[i__ * 2];
		    isuppz[i__ * 2] = isuppz[j * 2];
		    isuppz[j * 2] = itmp;
		}
	    }
/* L60: */
	}
    }

    work[1] = (real) lwmin;
    iwork[1] = liwmin;
    return 0;

/*     End of CSTEGR */

} /* cstegr_ */

#undef z___ref
#undef z___subscr