/* Subroutine */ int cla_lin_berr__(integer *n, integer *nz, integer *nrhs, 
	complex *res, real *ayb, real *berr)
{
    /* System generated locals */
    integer ayb_dim1, ayb_offset, res_dim1, res_offset, i__1, i__2, i__3, 
	    i__4;
    real r__1, r__2, r__3;
    complex q__1, q__2, q__3;

    /* Local variables */
    integer i__, j;
    real tmp, safe1;

/*     -- LAPACK routine (version 3.2.1)                                 -- */
/*     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- */
/*     -- Jason Riedy of Univ. of California Berkeley.                 -- */
/*     -- April 2009                                                   -- */

/*     -- LAPACK is a software package provided by Univ. of Tennessee, -- */
/*     -- Univ. of California Berkeley and NAG Ltd.                    -- */

/*  Purpose */
/*  ======= */

/*     CLA_LIN_BERR computes componentwise relative backward error from */
/*     the formula */
/*         max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) ) */
/*     where abs(Z) is the componentwise absolute value of the matrix */
/*     or vector Z. */

/*     N       (input) INTEGER */
/*     The number of linear equations, i.e., the order of the */
/*     matrix A.  N >= 0. */

/*     NZ      (input) INTEGER */
/*     We add (NZ+1)*SLAMCH( 'Safe minimum' ) to R(i) in the numerator to */
/*     guard against spuriously zero residuals. Default value is N. */

/*     NRHS    (input) INTEGER */
/*     The number of right hand sides, i.e., the number of columns */
/*     of the matrices AYB, RES, and BERR.  NRHS >= 0. */

/*     RES    (input) DOUBLE PRECISION array, dimension (N,NRHS) */
/*     The residual matrix, i.e., the matrix R in the relative backward */
/*     error formula above. */

/*     AYB    (input) DOUBLE PRECISION array, dimension (N, NRHS) */
/*     The denominator in the relative backward error formula above, i.e., */
/*     the matrix abs(op(A_s))*abs(Y) + abs(B_s). The matrices A, Y, and B */
/*     are from iterative refinement (see cla_gerfsx_extended.f). */

/*     RES    (output) COMPLEX array, dimension (NRHS) */
/*     The componentwise relative backward error from the formula above. */

/*  ===================================================================== */

/*     Adding SAFE1 to the numerator guards against spuriously zero */
/*     residuals.  A similar safeguard is in the CLA_yyAMV routine used */
/*     to compute AYB. */

    /* Parameter adjustments */
    --berr;
    ayb_dim1 = *n;
    ayb_offset = 1 + ayb_dim1;
    ayb -= ayb_offset;
    res_dim1 = *n;
    res_offset = 1 + res_dim1;
    res -= res_offset;

    /* Function Body */
    safe1 = slamch_("Safe minimum");
    safe1 = (*nz + 1) * safe1;
    i__1 = *nrhs;
    for (j = 1; j <= i__1; ++j) {
	berr[j] = 0.f;
	i__2 = *n;
	for (i__ = 1; i__ <= i__2; ++i__) {
	    if (ayb[i__ + j * ayb_dim1] != 0.f) {
		i__3 = i__ + j * res_dim1;
		r__3 = (r__1 = res[i__3].r, dabs(r__1)) + (r__2 = r_imag(&res[
			i__ + j * res_dim1]), dabs(r__2));
		q__3.r = r__3, q__3.i = 0.f;
		q__2.r = safe1 + q__3.r, q__2.i = q__3.i;
		i__4 = i__ + j * ayb_dim1;
		q__1.r = q__2.r / ayb[i__4], q__1.i = q__2.i / ayb[i__4];
		tmp = q__1.r;
/* Computing MAX */
		r__1 = berr[j];
		berr[j] = dmax(r__1,tmp);
	    }

/*     If AYB is exactly 0.0 (and if computed by CLA_yyAMV), then we know */
/*     the true residual also must be exactly 0.0. */

	}
    }
    return 0;
} /* cla_lin_berr__ */

//.........这里部分代码省略.........
    b -= b_offset;
    x_dim1 = *ldx;
    x_offset = 1 + x_dim1;
    x -= x_offset;
    xact_dim1 = *ldxact;
    xact_offset = 1 + xact_dim1;
    xact -= xact_offset;
    --ferr;
    --berr;
    --reslts;

    /* Function Body */
    if (*n <= 0 || *nrhs <= 0) {
	reslts[1] = 0.f;
	reslts[2] = 0.f;
	return 0;
    }

    eps = slamch_("Epsilon");
    unfl = slamch_("Safe minimum");
    ovfl = 1.f / unfl;
    notran = lsame_(trans, "N");

/*     Test 1:  Compute the maximum of */
/*        norm(X - XACT) / ( norm(X) * FERR ) */
/*     over all the vectors X and XACT using the infinity-norm. */

    errbnd = 0.f;
    if (*chkferr) {
	i__1 = *nrhs;
	for (j = 1; j <= i__1; ++j) {
	    imax = isamax_(n, &x[j * x_dim1 + 1], &c__1);
/* Computing MAX */
	    r__2 = (r__1 = x[imax + j * x_dim1], dabs(r__1));
	    xnorm = dmax(r__2,unfl);
	    diff = 0.f;
	    i__2 = *n;
	    for (i__ = 1; i__ <= i__2; ++i__) {
/* Computing MAX */
		r__2 = diff, r__3 = (r__1 = x[i__ + j * x_dim1] - xact[i__ + 
			j * xact_dim1], dabs(r__1));
		diff = dmax(r__2,r__3);
/* L10: */
	    }

	    if (xnorm > 1.f) {
		goto L20;
	    } else if (diff <= ovfl * xnorm) {
		goto L20;
	    } else {
		errbnd = 1.f / eps;
		goto L30;
	    }

L20:
	    if (diff / xnorm <= ferr[j]) {
/* Computing MAX */
		r__1 = errbnd, r__2 = diff / xnorm / ferr[j];
		errbnd = dmax(r__1,r__2);
	    } else {
		errbnd = 1.f / eps;
	    }
L30:
	    ;
	}
    }

//.........这里部分代码省略.........
					lda, &work[1], lwork, &rwork[1], &dif,
					 &iinfo);
			    } else if (ijvt == 2) {
				cunt03_("R", &n, &mnmin, &n, &mnmin, &vtsav[
					vtsav_offset], ldvt, &vt[vt_offset], 
					ldvt, &work[1], lwork, &rwork[1], &
					dif, &iinfo);
			    } else if (ijvt == 3) {
				cunt03_("R", &n, &n, &n, &mnmin, &vtsav[
					vtsav_offset], ldvt, &vt[vt_offset], 
					ldvt, &work[1], lwork, &rwork[1], &
					dif, &iinfo);
			    }
			}
			result[5] = dmax(result[5],dif);

/*                    Compare S */

			dif = 0.f;
/* Computing MAX */
			r__1 = (real) mnmin * ulp * s[1], r__2 = slamch_(
				"Safe minimum");
			div = dmax(r__1,r__2);
			i__3 = mnmin - 1;
			for (i__ = 1; i__ <= i__3; ++i__) {
			    if (ssav[i__] < ssav[i__ + 1]) {
				dif = ulpinv;
			    }
			    if (ssav[i__] < 0.f) {
				dif = ulpinv;
			    }
/* Computing MAX */
			    r__2 = dif, r__3 = (r__1 = ssav[i__] - s[i__], 
				    dabs(r__1)) / div;
			    dif = dmax(r__2,r__3);
/* L80: */
			}
			result[6] = dmax(result[6],dif);
L90:
			;
		    }
/* L100: */
		}

/*              Test for CGESDD */

		iwtmp = (mnmin << 1) * mnmin + (mnmin << 1) + max(m,n);
		lswork = iwtmp + (iwspc - 1) * (*lwork - iwtmp) / 3;
		lswork = min(lswork,*lwork);
		lswork = max(lswork,1);
		if (iwspc == 4) {
		    lswork = *lwork;
		}

/*              Factorize A */

		clacpy_("F", &m, &n, &asav[asav_offset], lda, &a[a_offset], 
			lda);
		cgesdd_("A", &m, &n, &a[a_offset], lda, &ssav[1], &usav[
			usav_offset], ldu, &vtsav[vtsav_offset], ldvt, &work[
			1], &lswork, &rwork[1], &iwork[1], &iinfo);
		if (iinfo != 0) {
		    io___39.ciunit = *nounit;
		    s_wsfe(&io___39);
		    do_fio(&c__1, "GESDD", (ftnlen)5);
		    do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));

//.........这里部分代码省略.........
	for (id2 = 1; id2 <= 4; ++id2) {
	    d2 = vdd[id2 - 1];
	    for (ica = 1; ica <= 5; ++ica) {
		ca = vca[ica - 1];
		for (itrans = 0; itrans <= 1; ++itrans) {
		    for (ismin = 1; ismin <= 4; ++ismin) {
			smin = vsmin[ismin - 1];

			na = 1;
			nw = 1;
			for (ia = 1; ia <= 3; ++ia) {
			    a_ref(1, 1) = vab[ia - 1];
			    for (ib = 1; ib <= 3; ++ib) {
				b_ref(1, 1) = vab[ib - 1];
				for (iwr = 1; iwr <= 4; ++iwr) {
				    if (d1 == 1.f && d2 == 1.f && ca == 1.f) {
					wr = vwr[iwr - 1] * a_ref(1, 1);
				    } else {
					wr = vwr[iwr - 1];
				    }
				    wi = 0.f;
				    slaln2_(&ltrans[itrans], &na, &nw, &smin, 
					    &ca, a, &c__2, &d1, &d2, b, &c__2,
					     &wr, &wi, x, &c__2, &scale, &
					    xnorm, &info);
				    if (info < 0) {
					++ninfo[1];
				    }
				    if (info > 0) {
					++ninfo[2];
				    }
				    res = (r__1 = (ca * a_ref(1, 1) - wr * d1)
					     * x_ref(1, 1) - scale * b_ref(1, 
					    1), dabs(r__1));
				    if (info == 0) {
/* Computing MAX */
					r__2 = eps * (r__1 = (ca * a_ref(1, 1)
						 - wr * d1) * x_ref(1, 1), 
						dabs(r__1));
					den = dmax(r__2,smlnum);
				    } else {
/* Computing MAX */
					r__2 = smin * (r__1 = x_ref(1, 1), 
						dabs(r__1));
					den = dmax(r__2,smlnum);
				    }
				    res /= den;
				    if ((r__1 = x_ref(1, 1), dabs(r__1)) < 
					    unfl && (r__3 = b_ref(1, 1), dabs(
					    r__3)) <= smlnum * (r__2 = ca * 
					    a_ref(1, 1) - wr * d1, dabs(r__2))
					    ) {
					res = 0.f;
				    }
				    if (scale > 1.f) {
					res += 1.f / eps;
				    }
				    res += (r__2 = xnorm - (r__1 = x_ref(1, 1)
					    , dabs(r__1)), dabs(r__2)) / dmax(
					    smlnum,xnorm) / eps;
				    if (info != 0 && info != 1) {
					res += 1.f / eps;
				    }
				    ++(*knt);
				    if (res > *rmax) {
					*lmax = *knt;

//.........这里部分代码省略.........
					    anrmpv = dmax(r__1,r__2);
/* L60: */
					}
/* L70: */
				    }
/* Computing MIN */
				    i__7 = info - 1, i__6 = kl + ku;
				    i__8 = min(i__7,i__6);
/* Computing MAX */
				    i__9 = 1, i__10 = kl + ku + 2 - info;
				    rpvgrw = clantb_("M", "U", "N", &info, &
					    i__8, &afb[max(i__9, i__10)], &
					    ldafb, rdum);
				    if (rpvgrw == 0.f) {
					rpvgrw = 1.f;
				    } else {
					rpvgrw = anrmpv / rpvgrw;
				    }
				} else {
				    i__8 = kl + ku;
				    rpvgrw = clantb_("M", "U", "N", &n, &i__8, 
					     &afb[1], &ldafb, rdum);
				    if (rpvgrw == 0.f) {
					rpvgrw = 1.f;
				    } else {
					rpvgrw = clangb_("M", &n, &kl, &ku, &
						a[1], &lda, rdum) /
						 rpvgrw;
				    }
				}
/* Computing MAX */
				r__2 = rwork[(*nrhs << 1) + 1];
				result[6] = (r__1 = rpvgrw - rwork[(*nrhs << 
					1) + 1], dabs(r__1)) / dmax(r__2,
					rpvgrw) / slamch_("E");

				if (! prefac) {

/*                             Reconstruct matrix from factors and */
/*                             compute residual. */

				    cgbt01_(&n, &n, &kl, &ku, &a[1], &lda, &
					    afb[1], &ldafb, &iwork[1], &work[
					    1], result);
				    k1 = 1;
				} else {
				    k1 = 2;
				}

				if (info == 0) {
				    trfcon = FALSE_;

/*                             Compute residual of the computed solution. */

				    clacpy_("Full", &n, nrhs, &bsav[1], &ldb, 
					    &work[1], &ldb);
				    cgbt02_(trans, &n, &n, &kl, &ku, nrhs, &
					    asav[1], &lda, &x[1], &ldb, &work[
					    1], &ldb, &result[1]);

/*                             Check solution from generated exact */
/*                             solution. */

				    if (nofact || prefac && lsame_(equed, 
					    "N")) {
					cget04_(&n, nrhs, &x[1], &ldb, &xact[

//.........这里部分代码省略.........
    if (kase != 0) {
	if (kase == kase1) {

/*           Multiply by inv(L). */

	    if (lnoti) {
		i__1 = *n - 1;
		for (j = 1; j <= i__1; ++j) {
/* Computing MIN */
		    i__2 = *kl, i__3 = *n - j;
		    lm = min(i__2,i__3);
		    jp = ipiv[j];
		    i__2 = jp;
		    t.r = work[i__2].r, t.i = work[i__2].i;
		    if (jp != j) {
			i__2 = jp;
			i__3 = j;
			work[i__2].r = work[i__3].r, work[i__2].i = work[i__3]
				.i;
			i__2 = j;
			work[i__2].r = t.r, work[i__2].i = t.i;
		    }
		    q__1.r = -t.r, q__1.i = -t.i;
		    caxpy_(&lm, &q__1, &ab[kd + 1 + j * ab_dim1], &c__1, &
			    work[j + 1], &c__1);
/* L20: */
		}
	    }

/*           Multiply by inv(U). */

	    i__1 = *kl + *ku;
	    clatbs_("Upper", "No transpose", "Non-unit", normin, n, &i__1, &
		    ab[ab_offset], ldab, &work[1], &scale, &rwork[1], info);
	} else {

/*           Multiply by inv(U'). */

	    i__1 = *kl + *ku;
	    clatbs_("Upper", "Conjugate transpose", "Non-unit", normin, n, &
		    i__1, &ab[ab_offset], ldab, &work[1], &scale, &rwork[1], 
		    info);

/*           Multiply by inv(L'). */

	    if (lnoti) {
		for (j = *n - 1; j >= 1; --j) {
/* Computing MIN */
		    i__1 = *kl, i__2 = *n - j;
		    lm = min(i__1,i__2);
		    i__1 = j;
		    i__2 = j;
		    cdotc_(&q__2, &lm, &ab[kd + 1 + j * ab_dim1], &c__1, &
			    work[j + 1], &c__1);
		    q__1.r = work[i__2].r - q__2.r, q__1.i = work[i__2].i - 
			    q__2.i;
		    work[i__1].r = q__1.r, work[i__1].i = q__1.i;
		    jp = ipiv[j];
		    if (jp != j) {
			i__1 = jp;
			t.r = work[i__1].r, t.i = work[i__1].i;
			i__1 = jp;
			i__2 = j;
			work[i__1].r = work[i__2].r, work[i__1].i = work[i__2]
				.i;
			i__1 = j;
			work[i__1].r = t.r, work[i__1].i = t.i;
		    }
/* L30: */
		}
	    }
	}

/*        Divide X by 1/SCALE if doing so will not cause overflow. */

	*(unsigned char *)normin = 'Y';
	if (scale != 1.f) {
	    ix = icamax_(n, &work[1], &c__1);
	    i__1 = ix;
	    if (scale < ((r__1 = work[i__1].r, dabs(r__1)) + (r__2 = r_imag(&
		    work[ix]), dabs(r__2))) * smlnum || scale == 0.f) {
		goto L40;
	    }
	    csrscl_(n, &scale, &work[1], &c__1);
	}
	goto L10;
    }

/*     Compute the estimate of the reciprocal condition number. */

    if (ainvnm != 0.f) {
	*rcond = 1.f / ainvnm / *anorm;
    }

L40:
    return 0;

/*     End of CGBCON */

} /* cgbcon_ */

//.........这里部分代码省略.........
			    l + i__) * q_dim1 + 1], &c__1, &csq, &snq);
		}

/* L10: */
	    }
/* L20: */
	}

	if (! upper) {

/*           The matrices A13 and B13 were lower triangular at the start */
/*           of the cycle, and are now upper triangular. */

/*           Convergence test: test the parallelism of the corresponding */
/*           rows of A and B. */

	    error = 0.f;
/* Computing MIN */
	    i__2 = *l, i__3 = *m - *k;
	    i__1 = min(i__2,i__3);
	    for (i__ = 1; i__ <= i__1; ++i__) {
		i__2 = *l - i__ + 1;
		ccopy_(&i__2, &a[*k + i__ + (*n - *l + i__) * a_dim1], lda, &
			work[1], &c__1);
		i__2 = *l - i__ + 1;
		ccopy_(&i__2, &b[i__ + (*n - *l + i__) * b_dim1], ldb, &work[*
			l + 1], &c__1);
		i__2 = *l - i__ + 1;
		clapll_(&i__2, &work[1], &c__1, &work[*l + 1], &c__1, &ssmin);
		error = dmax(error,ssmin);
/* L30: */
	    }

	    if (dabs(error) <= dmin(*tola,*tolb)) {
		goto L50;
	    }
	}

/*        End of cycle loop */

/* L40: */
    }

/*     The algorithm has not converged after MAXIT cycles. */

    *info = 1;
    goto L100;

L50:

/*     If ERROR <= MIN(TOLA,TOLB), then the algorithm has converged. */
/*     Compute the generalized singular value pairs (ALPHA, BETA), and */
/*     set the triangular matrix R to array A. */

    i__1 = *k;
    for (i__ = 1; i__ <= i__1; ++i__) {
	alpha[i__] = 1.f;
	beta[i__] = 0.f;
/* L60: */
    }

/* Computing MIN */
    i__2 = *l, i__3 = *m - *k;
    i__1 = min(i__2,i__3);
    for (i__ = 1; i__ <= i__1; ++i__) {

//.........这里部分代码省略.........

    i__1 = *nrhs;
    for (j = 1; j <= *nrhs; ++j) {

	count = 1;
	lstres = 3.f;
L20:

/*        Loop until stopping criterion is satisfied.   

          Compute residual R = B - A * X */

	ccopy_(n, &B(1,j), &c__1, &WORK(1), &c__1);
	q__1.r = -1.f, q__1.i = 0.f;
	chbmv_(uplo, n, kd, &q__1, &AB(1,1), ldab, &X(1,j), &
		c__1, &c_b1, &WORK(1), &c__1);

/*        Compute componentwise relative backward error from formula 
  

          max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )   

          where abs(Z) is the componentwise absolute value of the matr
ix   
          or vector Z.  If the i-th component of the denominator is le
ss   
          than SAFE2, then SAFE1 is added to the i-th components of th
e   
          numerator and denominator before dividing. */

	i__2 = *n;
	for (i = 1; i <= *n; ++i) {
	    i__3 = i + j * b_dim1;
	    RWORK(i) = (r__1 = B(i,j).r, dabs(r__1)) + (r__2 = r_imag(&B(i,j)), dabs(r__2));
/* L30: */
	}

/*        Compute abs(A)*abs(X) + abs(B). */

	if (upper) {
	    i__2 = *n;
	    for (k = 1; k <= *n; ++k) {
		s = 0.f;
		i__3 = k + j * x_dim1;
		xk = (r__1 = X(k,j).r, dabs(r__1)) + (r__2 = r_imag(&X(k,j)), dabs(r__2));
		l = *kd + 1 - k;
/* Computing MAX */
		i__3 = 1, i__4 = k - *kd;
		i__5 = k - 1;
		for (i = max(1,k-*kd); i <= k-1; ++i) {
		    i__3 = l + i + k * ab_dim1;
		    RWORK(i) += ((r__1 = AB(l+i,k).r, dabs(r__1)) + (r__2 = 
			    r_imag(&AB(l+i,k)), dabs(r__2))) * 
			    xk;
		    i__3 = l + i + k * ab_dim1;
		    i__4 = i + j * x_dim1;
		    s += ((r__1 = AB(l+i,k).r, dabs(r__1)) + (r__2 = r_imag(&
			    AB(l+i,k)), dabs(r__2))) * ((r__3 = 
			    X(i,j).r, dabs(r__3)) + (r__4 = r_imag(&X(i,j)), dabs(r__4)));
/* L40: */
		}
		i__5 = *kd + 1 + k * ab_dim1;
		RWORK(k) = RWORK(k) + (r__1 = AB(*kd+1,k).r, dabs(r__1)) * xk + 
			s;
/* L50: */
	    }

real sasum_(integer *n, real *sx, integer *incx)
{


    /* System generated locals */
    integer i__1, i__2;
    real ret_val, r__1, r__2, r__3, r__4, r__5, r__6;

    /* Local variables */
    integer i, m, nincx;
    real stemp;
    integer mp1;


/*     takes the sum of the absolute values.   
       uses unrolled loops for increment equal to one.   
       jack dongarra, linpack, 3/11/78.   
       modified 3/93 to return if incx .le. 0.   
       modified 12/3/93, array(1) declarations changed to array(*)   


    
   Parameter adjustments   
       Function Body */
#define SX(I) sx[(I)-1]


    ret_val = 0.f;
    stemp = 0.f;
    if (*n <= 0 || *incx <= 0) {
	return ret_val;
    }
    if (*incx == 1) {
	goto L20;
    }

/*        code for increment not equal to 1 */

    nincx = *n * *incx;
    i__1 = nincx;
    i__2 = *incx;
    for (i = 1; *incx < 0 ? i >= nincx : i <= nincx; i += *incx) {
	stemp += (r__1 = SX(i), dabs(r__1));
/* L10: */
    }
    ret_val = stemp;
    return ret_val;

/*        code for increment equal to 1   


          clean-up loop */

L20:
    m = *n % 6;
    if (m == 0) {
	goto L40;
    }
    i__2 = m;
    for (i = 1; i <= m; ++i) {
	stemp += (r__1 = SX(i), dabs(r__1));
/* L30: */
    }
    if (*n < 6) {
	goto L60;
    }
L40:
    mp1 = m + 1;
    i__2 = *n;
    for (i = mp1; i <= *n; i += 6) {
	stemp = stemp + (r__1 = SX(i), dabs(r__1)) + (r__2 = SX(i + 1), dabs(
		r__2)) + (r__3 = SX(i + 2), dabs(r__3)) + (r__4 = SX(i + 3), 
		dabs(r__4)) + (r__5 = SX(i + 4), dabs(r__5)) + (r__6 = SX(i + 
		5), dabs(r__6));
/* L50: */
    }
L60:
    ret_val = stemp;
    return ret_val;
} /* sasum_ */

//.........这里部分代码省略.........
	    z__[(*n0 << 2) - *pp] = dmin(r__1,r__2);
/* Computing MAX */
	    r__1 = *qmax, r__2 = z__[(*i0 << 2) + *pp - 3], r__1 = max(r__1,
		    r__2), r__2 = z__[(*i0 << 2) + *pp + 1];
	    *qmax = dmax(r__1,r__2);
	    *dmin__ = -0.f;
	}
    }

/*     Choose a shift. */

    slasq4_(i0, n0, &z__[1], pp, &n0in, dmin__, dmin1, dmin2, dn, dn1, dn2, 
	    tau, ttype, g);

/*     Call dqds until DMIN > 0. */

L70:

    slasq5_(i0, n0, &z__[1], pp, tau, dmin__, dmin1, dmin2, dn, dn1, dn2, 
	    ieee);

    *ndiv += *n0 - *i0 + 2;
    ++(*iter);

/*     Check status. */

    if (*dmin__ >= 0.f && *dmin1 > 0.f) {

/*        Success. */

	goto L90;

    } else if (*dmin__ < 0.f && *dmin1 > 0.f && z__[(*n0 - 1 << 2) - *pp] < 
	    tol * (*sigma + *dn1) && dabs(*dn) < tol * *sigma) {

/*        Convergence hidden by negative DN. */

	z__[(*n0 - 1 << 2) - *pp + 2] = 0.f;
	*dmin__ = 0.f;
	goto L90;
    } else if (*dmin__ < 0.f) {

/*        TAU too big. Select new TAU and try again. */

	++(*nfail);
	if (*ttype < -22) {

/*           Failed twice. Play it safe. */

	    *tau = 0.f;
	} else if (*dmin1 > 0.f) {

/*           Late failure. Gives excellent shift. */

	    *tau = (*tau + *dmin__) * (1.f - eps * 2.f);
	    *ttype += -11;
	} else {

/*           Early failure. Divide by 4. */

	    *tau *= .25f;
	    *ttype += -12;
	}
	goto L70;
    } else if (sisnan_(dmin__)) {

//.........这里部分代码省略.........
	}
	t = c_abs(&a[k + (k + 1) * a_dim1]);
	i__2 = k + 1 + (k + 1) * a_dim1;
	d__ = d__ / t * a[i__2].r - t;
	goto L40;
L30:
	d__ = t;
	t = 0.f;
L40:
L50:

	if (noert) {
	    goto L60;
	}
	if (d__ > 0.f) {
	    ++inert[1];
	}
	if (d__ < 0.f) {
	    ++inert[2];
	}
	if (d__ == 0.f) {
	    ++inert[3];
	}
L60:

	if (nodet) {
	    goto L120;
	}
	det[1] = d__ * det[1];
	if (det[1] == 0.f) {
	    goto L110;
	}
L70:
	if (dabs(det[1]) >= 1.f) {
	    goto L80;
	}
	det[1] = ten * det[1];
	det[2] += -1.f;
	goto L70;
L80:
L90:
	if (dabs(det[1]) < ten) {
	    goto L100;
	}
	det[1] /= ten;
	det[2] += 1.f;
	goto L90;
L100:
L110:
L120:
/* L130: */
	;
    }
L140:

/*     COMPUTE INVERSE(A) */

    if (noinv) {
	goto L270;
    }
    k = 1;
L150:
    if (k > *n) {
	goto L260;
    }
    km1 = k - 1;

//.........这里部分代码省略.........

/*     Determine the unit roundoff for this environment. */

    eps = slamch_("E");
/* Computing 2nd power */
    r__1 = eps;
    eps2 = r__1 * r__1;
    safmin = slamch_("S");
    safmax = 1.f / safmin;
    ssfmax = sqrt(safmax) / 3.f;
    ssfmin = sqrt(safmin) / eps2;

/*     Compute the eigenvalues of the tridiagonal matrix. */

    nmaxit = *n * 30;
    sigma = 0.f;
    jtot = 0;

/*     Determine where the matrix splits and choose QL or QR iteration   
       for each block, according to whether top or bottom diagonal   
       element is smaller. */

    l1 = 1;

L10:
    if (l1 > *n) {
	goto L170;
    }
    if (l1 > 1) {
	e[l1 - 1] = 0.f;
    }
    i__1 = *n - 1;
    for (m = l1; m <= i__1; ++m) {
	if ((r__3 = e[m], dabs(r__3)) <= sqrt((r__1 = d__[m], dabs(r__1))) * 
		sqrt((r__2 = d__[m + 1], dabs(r__2))) * eps) {
	    e[m] = 0.f;
	    goto L30;
	}
/* L20: */
    }
    m = *n;

L30:
    l = l1;
    lsv = l;
    lend = m;
    lendsv = lend;
    l1 = m + 1;
    if (lend == l) {
	goto L10;
    }

/*     Scale submatrix in rows and columns L to LEND */

    i__1 = lend - l + 1;
    anorm = slanst_("I", &i__1, &d__[l], &e[l]);
    iscale = 0;
    if (anorm > ssfmax) {
	iscale = 1;
	i__1 = lend - l + 1;
	slascl_("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n, 
		info);
	i__1 = lend - l;
	slascl_("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n, 
		info);
    } else if (anorm < ssfmin) {

//.........这里部分代码省略.........
    a_offset = 1 + a_dim1;
    a -= a_offset;
    --work;
    --rwork;

    /* Function Body */
    *info = 0;
    onenrm = *(unsigned char *)norm == '1' || lsame_(norm, "O");
    if (! onenrm && ! lsame_(norm, "I")) {
	*info = -1;
    } else if (*n < 0) {
	*info = -2;
    } else if (*lda < max(1,*n)) {
	*info = -4;
    } else if (*anorm < 0.f) {
	*info = -5;
    }
    if (*info != 0) {
	i__1 = -(*info);
	xerbla_("CGECON", &i__1);
	return 0;
    }

/*     Quick return if possible */

    *rcond = 0.f;
    if (*n == 0) {
	*rcond = 1.f;
	return 0;
    } else if (*anorm == 0.f) {
	return 0;
    }

    smlnum = slamch_("Safe minimum");

/*     Estimate the norm of inv(A). */

    ainvnm = 0.f;
    *(unsigned char *)normin = 'N';
    if (onenrm) {
	kase1 = 1;
    } else {
	kase1 = 2;
    }
    kase = 0;
L10:
    clacn2_(n, &work[*n + 1], &work[1], &ainvnm, &kase, isave);
    if (kase != 0) {
	if (kase == kase1) {

/*           Multiply by inv(L). */

	    clatrs_("Lower", "No transpose", "Unit", normin, n, &a[a_offset], 
		    lda, &work[1], &sl, &rwork[1], info);

/*           Multiply by inv(U). */

	    clatrs_("Upper", "No transpose", "Non-unit", normin, n, &a[
		    a_offset], lda, &work[1], &su, &rwork[*n + 1], info);
	} else {

/*           Multiply by inv(U'). */

	    clatrs_("Upper", "Conjugate transpose", "Non-unit", normin, n, &a[
		    a_offset], lda, &work[1], &su, &rwork[*n + 1], info);

/*           Multiply by inv(L'). */

	    clatrs_("Lower", "Conjugate transpose", "Unit", normin, n, &a[
		    a_offset], lda, &work[1], &sl, &rwork[1], info);
	}

/*        Divide X by 1/(SL*SU) if doing so will not cause overflow. */

	scale = sl * su;
	*(unsigned char *)normin = 'Y';
	if (scale != 1.f) {
	    ix = icamax_(n, &work[1], &c__1);
	    i__1 = ix;
	    if (scale < ((r__1 = work[i__1].r, dabs(r__1)) + (r__2 = r_imag(&
		    work[ix]), dabs(r__2))) * smlnum || scale == 0.f) {
		goto L20;
	    }
	    csrscl_(n, &scale, &work[1], &c__1);
	}
	goto L10;
    }

/*     Compute the estimate of the reciprocal condition number. */

    if (ainvnm != 0.f) {
	*rcond = 1.f / ainvnm / *anorm;
    }

L20:
    return 0;

/*     End of CGECON */

} /* cgecon_ */

/* Subroutine */ int slassq_(integer *n, real *x, integer *incx, real *scale, 
	real *sumsq)
{
    /* System generated locals */
    integer i__1, i__2;
    real r__1;

    /* Local variables */
    integer ix;
    real absxi;

/*  -- LAPACK auxiliary routine (version 3.2) -- */
/*     November 2006 */

/*  Purpose */
/*  ======= */

/*  SLASSQ  returns the values  scl  and  smsq  such that */

/*  where  x( i ) = X( 1 + ( i - 1 )*INCX ). The value of  sumsq  is */
/*  assumed to be non-negative and  scl  returns the value */

/*     scl = max( scale, abs( x( i ) ) ). */

/*  scale and sumsq must be supplied in SCALE and SUMSQ and */
/*  scl and smsq are overwritten on SCALE and SUMSQ respectively. */

/*  The routine makes only one pass through the vector x. */

/*  Arguments */
/*  ========= */

/*  N       (input) INTEGER */
/*          The number of elements to be used from the vector X. */

/*  X       (input) REAL array, dimension (N) */
/*          The vector for which a scaled sum of squares is computed. */
/*             x( i )  = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n. */

/*  INCX    (input) INTEGER */
/*          The increment between successive values of the vector X. */
/*          INCX > 0. */

/*  SCALE   (input/output) REAL */
/*          On entry, the value  scale  in the equation above. */
/*          On exit, SCALE is overwritten with  scl , the scaling factor */
/*          for the sum of squares. */

/*  SUMSQ   (input/output) REAL */
/*          On entry, the value  sumsq  in the equation above. */
/*          On exit, SUMSQ is overwritten with  smsq , the basic sum of */
/*          squares from which  scl  has been factored out. */

/* ===================================================================== */

    /* Parameter adjustments */
    --x;

    /* Function Body */
    if (*n > 0) {
	i__1 = (*n - 1) * *incx + 1;
	i__2 = *incx;
	for (ix = 1; i__2 < 0 ? ix >= i__1 : ix <= i__1; ix += i__2) {
	    if (x[ix] != 0.f) {
		absxi = (r__1 = x[ix], dabs(r__1));
		if (*scale < absxi) {
/* Computing 2nd power */
		    r__1 = *scale / absxi;
		    *sumsq = *sumsq * (r__1 * r__1) + 1;
		    *scale = absxi;
		} else {
/* Computing 2nd power */
		    r__1 = absxi / *scale;
		    *sumsq += r__1 * r__1;
		}
	    }
	}
    }
    return 0;

/*     End of SLASSQ */

} /* slassq_ */

//.........这里部分代码省略.........
    b -= b_offset;
    x_dim1 = *ldx;
    x_offset = 1 + x_dim1;
    x -= x_offset;
    xact_dim1 = *ldxact;
    xact_offset = 1 + xact_dim1;
    xact -= xact_offset;
    --ferr;
    --berr;
    --reslts;

    /* Function Body */
    if (*n <= 0 || *nrhs <= 0) {
	reslts[1] = 0.f;
	reslts[2] = 0.f;
	return 0;
    }

    eps = slamch_("Epsilon");
    unfl = slamch_("Safe minimum");
    ovfl = 1.f / unfl;
    upper = lsame_(uplo, "U");

/*     Test 1:  Compute the maximum of */
/*        norm(X - XACT) / ( norm(X) * FERR ) */
/*     over all the vectors X and XACT using the infinity-norm. */

    errbnd = 0.f;
    i__1 = *nrhs;
    for (j = 1; j <= i__1; ++j) {
	imax = icamax_(n, &x[j * x_dim1 + 1], &c__1);
/* Computing MAX */
	i__2 = imax + j * x_dim1;
	r__3 = (r__1 = x[i__2].r, dabs(r__1)) + (r__2 = r_imag(&x[imax + j * 
		x_dim1]), dabs(r__2));
	xnorm = dmax(r__3,unfl);
	diff = 0.f;
	i__2 = *n;
	for (i__ = 1; i__ <= i__2; ++i__) {
	    i__3 = i__ + j * x_dim1;
	    i__4 = i__ + j * xact_dim1;
	    q__2.r = x[i__3].r - xact[i__4].r, q__2.i = x[i__3].i - xact[i__4]
		    .i;
	    q__1.r = q__2.r, q__1.i = q__2.i;
/* Computing MAX */
	    r__3 = diff, r__4 = (r__1 = q__1.r, dabs(r__1)) + (r__2 = r_imag(&
		    q__1), dabs(r__2));
	    diff = dmax(r__3,r__4);
/* L10: */
	}

	if (xnorm > 1.f) {
	    goto L20;
	} else if (diff <= ovfl * xnorm) {
	    goto L20;
	} else {
	    errbnd = 1.f / eps;
	    goto L30;
	}

L20:
	if (diff / xnorm <= ferr[j]) {
/* Computing MAX */
	    r__1 = errbnd, r__2 = diff / xnorm / ferr[j];
	    errbnd = dmax(r__1,r__2);
	} else {

//.........这里部分代码省略.........
/* ***routines called  qcheb,qk15w,qwgtc */
/* ***end prologue  qc25c */




/*           the vector x contains the values cos(k*pi/24), */
/*           k = 1, ..., 11, to be used for the chebyshev series */
/*           expansion of f */


/*           list of major variables */
/*           ---------------------- */
/*           fval   - value of the function f at the points */
/*                    cos(k*pi/24),  k = 0, ..., 24 */
/*           cheb12 - chebyshev series expansion coefficients, */
/*                    for the function f, of degree 12 */
/*           cheb24 - chebyshev series expansion coefficients, */
/*                    for the function f, of degree 24 */
/*           res12  - approximation to the integral corresponding */
/*                    to the use of cheb12 */
/*           res24  - approximation to the integral corresponding */
/*                    to the use of cheb24 */
/*           qwgtc - external function subprogram defining */
/*                    the weight function */
/*           hlgth  - half-length of the interval */
/*           centr  - mid point of the interval */


/*           check the position of c. */

/* ***first executable statement  qc25c */
    cc = ((float)2. * *c__ - *b - *a) / (*b - *a);
    if (dabs(cc) < (float)1.1) {
	goto L10;
    }

/*           apply the 15-point gauss-kronrod scheme. */

    --(*krul);
    qk15w_((E_fp)f, (E_fp)qwgtc_, c__, &p2, &p3, &p4, &kp, a, b, result, 
	    abserr, &resabs, &resasc);
    *neval = 15;
    if (resasc == *abserr) {
	++(*krul);
    }
    goto L50;

/*           use the generalized clenshaw-curtis method. */

L10:
    hlgth = (*b - *a) * (float).5;
    centr = (*b + *a) * (float).5;
    *neval = 25;
    r__1 = hlgth + centr;
    fval[0] = (*f)(&r__1) * (float).5;
    fval[12] = (*f)(&centr);
    r__1 = centr - hlgth;
    fval[24] = (*f)(&r__1) * (float).5;
    for (i__ = 2; i__ <= 12; ++i__) {
	u = hlgth * x[i__ - 2];
	isym = 26 - i__;
	r__1 = u + centr;
	fval[i__ - 1] = (*f)(&r__1);
	r__1 = centr - u;
	fval[isym - 1] = (*f)(&r__1);

//.........这里部分代码省略.........

L120:
    i__1 = l;
    for (i__ = k; i__ <= i__1; ++i__) {
	scale[i__] = 1.f;
/* L130: */
    }

    if (lsame_(job, "P", (ftnlen)1, (ftnlen)1)) {
	goto L210;
    }

/*     Balance the submatrix in rows K to L. */

/*     Iterative loop for norm reduction */

    sfmin1 = slamch_("S", (ftnlen)1) / slamch_("P", (ftnlen)1);
    sfmax1 = 1.f / sfmin1;
    sfmin2 = sfmin1 * 8.f;
    sfmax2 = 1.f / sfmin2;
L140:
    noconv = FALSE_;

    i__1 = l;
    for (i__ = k; i__ <= i__1; ++i__) {
	c__ = 0.f;
	r__ = 0.f;

	i__2 = l;
	for (j = k; j <= i__2; ++j) {
	    if (j == i__) {
		goto L150;
	    }
	    c__ += (r__1 = a[j + i__ * a_dim1], dabs(r__1));
	    r__ += (r__1 = a[i__ + j * a_dim1], dabs(r__1));
L150:
	    ;
	}
	ica = isamax_(&l, &a[i__ * a_dim1 + 1], &c__1);
	ca = (r__1 = a[ica + i__ * a_dim1], dabs(r__1));
	i__2 = *n - k + 1;
	ira = isamax_(&i__2, &a[i__ + k * a_dim1], lda);
	ra = (r__1 = a[i__ + (ira + k - 1) * a_dim1], dabs(r__1));

/*        Guard against zero C or R due to underflow. */

	if (c__ == 0.f || r__ == 0.f) {
	    goto L200;
	}
	g = r__ / 8.f;
	f = 1.f;
	s = c__ + r__;
L160:
/* Computing MAX */
	r__1 = max(f,c__);
/* Computing MIN */
	r__2 = min(r__,g);
	if (c__ >= g || dmax(r__1,ca) >= sfmax2 || dmin(r__2,ra) <= sfmin2) {
	    goto L170;
	}
	f *= 8.f;
	c__ *= 8.f;
	ca *= 8.f;
	r__ /= 8.f;
	g /= 8.f;
	ra /= 8.f;

//.........这里部分代码省略.........

	    ip = 1;
	    i__1 = *n - 1;
	    for (j = 1; j <= i__1; ++j) {
		i__2 = *n - j;
		cnorm[j] = scasum_(&i__2, &ap[ip + 1], &c__1);
		ip = ip + *n - j + 1;
/* L20: */
	    }
	    cnorm[*n] = 0.f;
	}
    }

/*     Scale the column norms by TSCAL if the maximum element in CNORM is   
       greater than BIGNUM/2. */

    imax = isamax_(n, &cnorm[1], &c__1);
    tmax = cnorm[imax];
    if (tmax <= bignum * .5f) {
	tscal = 1.f;
    } else {
	tscal = .5f / (smlnum * tmax);
	sscal_(n, &tscal, &cnorm[1], &c__1);
    }

/*     Compute a bound on the computed solution vector to see if the   
       Level 2 BLAS routine CTPSV can be used. */

    xmax = 0.f;
    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
/* Computing MAX */
	i__2 = j;
	r__3 = xmax, r__4 = (r__1 = x[i__2].r / 2.f, dabs(r__1)) + (r__2 = 
		r_imag(&x[j]) / 2.f, dabs(r__2));
	xmax = dmax(r__3,r__4);
/* L30: */
    }
    xbnd = xmax;
    if (notran) {

/*        Compute the growth in A * x = b. */

	if (upper) {
	    jfirst = *n;
	    jlast = 1;
	    jinc = -1;
	} else {
	    jfirst = 1;
	    jlast = *n;
	    jinc = 1;
	}

	if (tscal != 1.f) {
	    grow = 0.f;
	    goto L60;
	}

	if (nounit) {

/*           A is non-unit triangular.   

             Compute GROW = 1/G(j) and XBND = 1/M(j).   
             Initially, G(0) = max{x(i), i=1,...,n}. */

	    grow = .5f / dmax(xbnd,smlnum);

//.........这里部分代码省略.........
    b_offset = 1 + b_dim1 * 1;
    b -= b_offset;
    x_dim1 = *ldx;
    x_offset = 1 + x_dim1 * 1;
    x -= x_offset;
    xact_dim1 = *ldxact;
    xact_offset = 1 + xact_dim1 * 1;
    xact -= xact_offset;
    --ferr;
    --berr;
    --reslts;

    /* Function Body */
    if (*n <= 0 || *nrhs <= 0) {
	reslts[1] = 0.f;
	reslts[2] = 0.f;
	return 0;
    }

    eps = slamch_("Epsilon");
    unfl = slamch_("Safe minimum");
    ovfl = 1.f / unfl;
    notran = lsame_(trans, "N");

/*     Test 1:  Compute the maximum of   
          norm(X - XACT) / ( norm(X) * FERR )   
       over all the vectors X and XACT using the infinity-norm. */

    errbnd = 0.f;
    i__1 = *nrhs;
    for (j = 1; j <= i__1; ++j) {
	imax = isamax_(n, &x_ref(1, j), &c__1);
/* Computing MAX */
	r__2 = (r__1 = x_ref(imax, j), dabs(r__1));
	xnorm = dmax(r__2,unfl);
	diff = 0.f;
	i__2 = *n;
	for (i__ = 1; i__ <= i__2; ++i__) {
/* Computing MAX */
	    r__2 = diff, r__3 = (r__1 = x_ref(i__, j) - xact_ref(i__, j), 
		    dabs(r__1));
	    diff = dmax(r__2,r__3);
/* L10: */
	}

	if (xnorm > 1.f) {
	    goto L20;
	} else if (diff <= ovfl * xnorm) {
	    goto L20;
	} else {
	    errbnd = 1.f / eps;
	    goto L30;
	}

L20:
	if (diff / xnorm <= ferr[j]) {
/* Computing MAX */
	    r__1 = errbnd, r__2 = diff / xnorm / ferr[j];
	    errbnd = dmax(r__1,r__2);
	} else {
	    errbnd = 1.f / eps;
	}
L30:
	;
    }
    reslts[1] = errbnd;

//.........这里部分代码省略.........
    nb = ilaenv_(&c__1, "SSTEBZ", " ", n, &c_n1, &c_n1, &c_n1, 6L, 1L);
    if (nb <= 1) {
	nb = 0;
    }

/*     Special Case when N=1 */

    if (*n == 1) {
	*nsplit = 1;
	ISPLIT(1) = 1;
	if (irange == 2 && (*vl >= D(1) || *vu < D(1))) {
	    *m = 0;
	} else {
	    W(1) = D(1);
	    IBLOCK(1) = 1;
	    *m = 1;
	}
	return 0;
    }

/*     Compute Splitting Points */

    *nsplit = 1;
    WORK(*n) = 0.f;
    pivmin = 1.f;

    i__1 = *n;
    for (j = 2; j <= *n; ++j) {
/* Computing 2nd power */
	r__1 = E(j - 1);
	tmp1 = r__1 * r__1;
/* Computing 2nd power */
	r__2 = ulp;
	if ((r__1 = D(j) * D(j - 1), dabs(r__1)) * (r__2 * r__2) + safemn > 
		tmp1) {
	    ISPLIT(*nsplit) = j - 1;
	    ++(*nsplit);
	    WORK(j - 1) = 0.f;
	} else {
	    WORK(j - 1) = tmp1;
	    pivmin = dmax(pivmin,tmp1);
	}
/* L10: */
    }
    ISPLIT(*nsplit) = *n;
    pivmin *= safemn;

/*     Compute Interval and ATOLI */

    if (ir