//.........这里部分代码省略.........
/*          The leading dimension of the array AB.  LDAB >= K+1. */

/*  WORK    (workspace) REAL array, dimension (MAX(1,LWORK)), */
/*          where LWORK >= N when NORM = 'I'; otherwise, WORK is not */
/*          referenced. */

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

    /* Parameter adjustments */
    ab_dim1 = *ldab;
    ab_offset = 1 + ab_dim1;
    ab -= ab_offset;
    --work;

    /* Function Body */
    if (*n == 0) {
	value = 0.f;
    } else if (lsame_(norm, "M")) {

/*        Find max(abs(A(i,j))). */

	if (lsame_(diag, "U")) {
	    value = 1.f;
	    if (lsame_(uplo, "U")) {
		i__1 = *n;
		for (j = 1; j <= i__1; ++j) {
/* Computing MAX */
		    i__2 = *k + 2 - j;
		    i__3 = *k;
		    for (i__ = max(i__2,1); i__ <= i__3; ++i__) {
/* Computing MAX */
			r__2 = value, r__3 = (r__1 = ab[i__ + j * ab_dim1], 
				dabs(r__1));
			value = dmax(r__2,r__3);
		    }
		}
	    } else {
		i__1 = *n;
		for (j = 1; j <= i__1; ++j) {
/* Computing MIN */
		    i__2 = *n + 1 - j, i__4 = *k + 1;
		    i__3 = min(i__2,i__4);
		    for (i__ = 2; i__ <= i__3; ++i__) {
/* Computing MAX */
			r__2 = value, r__3 = (r__1 = ab[i__ + j * ab_dim1], 
				dabs(r__1));
			value = dmax(r__2,r__3);
		    }
		}
	    }
	} else {
	    value = 0.f;
	    if (lsame_(uplo, "U")) {
		i__1 = *n;
		for (j = 1; j <= i__1; ++j) {
/* Computing MAX */
		    i__3 = *k + 2 - j;
		    i__2 = *k + 1;
		    for (i__ = max(i__3,1); i__ <= i__2; ++i__) {
/* Computing MAX */
			r__2 = value, r__3 = (r__1 = ab[i__ + j * ab_dim1], 
				dabs(r__1));
			value = dmax(r__2,r__3);
		    }
		}
	    } else {

//.........这里部分代码省略.........
			nab[ji + (nab_dim1 << 1)] = itmp1;
		    } else {
			*info = *mmax + 1;
			return 0;
		    }
		} else {

/*                 IJOB=3: Binary search.  Keep only the interval */
/*                         containing  w  s.t. N(w) = NVAL */

		    if (itmp1 <= nval[ji]) {
			ab[ji + ab_dim1] = tmp1;
			nab[ji + nab_dim1] = itmp1;
		    }
		    if (itmp1 >= nval[ji]) {
			ab[ji + (ab_dim1 << 1)] = tmp1;
			nab[ji + (nab_dim1 << 1)] = itmp1;
		    }
		}
/* L100: */
	    }
	    kl = klnew;

/*           End of Serial Version of the loop */

	}

/*        Check for convergence */

	kfnew = kf;
	i__2 = kl;
	for (ji = kf; ji <= i__2; ++ji) {
	    tmp1 = (r__1 = ab[ji + (ab_dim1 << 1)] - ab[ji + ab_dim1], dabs(
		    r__1));
/* Computing MAX */
	    r__3 = (r__1 = ab[ji + (ab_dim1 << 1)], dabs(r__1)), r__4 = (r__2 
		    = ab[ji + ab_dim1], dabs(r__2));
	    tmp2 = dmax(r__3,r__4);
/* Computing MAX */
	    r__1 = max(*abstol,*pivmin), r__2 = *reltol * tmp2;
	    if (tmp1 < dmax(r__1,r__2) || nab[ji + nab_dim1] >= nab[ji + (
		    nab_dim1 << 1)]) {

/*              Converged -- Swap with position KFNEW, */
/*                           then increment KFNEW */

		if (ji > kfnew) {
		    tmp1 = ab[ji + ab_dim1];
		    tmp2 = ab[ji + (ab_dim1 << 1)];
		    itmp1 = nab[ji + nab_dim1];
		    itmp2 = nab[ji + (nab_dim1 << 1)];
		    ab[ji + ab_dim1] = ab[kfnew + ab_dim1];
		    ab[ji + (ab_dim1 << 1)] = ab[kfnew + (ab_dim1 << 1)];
		    nab[ji + nab_dim1] = nab[kfnew + nab_dim1];
		    nab[ji + (nab_dim1 << 1)] = nab[kfnew + (nab_dim1 << 1)];
		    ab[kfnew + ab_dim1] = tmp1;
		    ab[kfnew + (ab_dim1 << 1)] = tmp2;
		    nab[kfnew + nab_dim1] = itmp1;
		    nab[kfnew + (nab_dim1 << 1)] = itmp2;
		    if (*ijob == 3) {
			itmp1 = nval[ji];
			nval[ji] = nval[kfnew];
			nval[kfnew] = itmp1;
		    }
		}
		++kfnew;
	    }
/* L110: */
	}
	kf = kfnew;

/*        Choose Midpoints */

	i__2 = kl;
	for (ji = kf; ji <= i__2; ++ji) {
	    c__[ji] = (ab[ji + ab_dim1] + ab[ji + (ab_dim1 << 1)]) * .5f;
/* L120: */
	}

/*        If no more intervals to refine, quit. */

	if (kf > kl) {
	    goto L140;
	}
/* L130: */
    }

/*     Converged */

L140:
/* Computing MAX */
    i__1 = kl + 1 - kf;
    *info = max(i__1,0);
    *mout = kl;

    return 0;

/*     End of SLAEBZ */

} /* slaebz_ */

/* Subroutine */ int sbdt01_(integer *m, integer *n, integer *kd, real *a, 
	integer *lda, real *q, integer *ldq, real *d__, real *e, real *pt, 
	integer *ldpt, real *work, real *resid)
{
    /* System generated locals */
    integer a_dim1, a_offset, pt_dim1, pt_offset, q_dim1, q_offset, i__1, 
	    i__2;
    real r__1, r__2;

    /* Local variables */
    integer i__, j;
    real eps, anorm;
    extern /* Subroutine */ int sgemv_(char *, integer *, integer *, real *, 
	    real *, integer *, real *, integer *, real *, real *, integer *);
    extern doublereal sasum_(integer *, real *, integer *);
    extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 
	    integer *);
    extern doublereal slamch_(char *), slange_(char *, integer *, 
	    integer *, real *, integer *, real *);


/*  -- LAPACK test routine (version 3.1) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/*     November 2006 */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */

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

/*  SBDT01 reconstructs a general matrix A from its bidiagonal form */
/*     A = Q * B * P' */
/*  where Q (m by min(m,n)) and P' (min(m,n) by n) are orthogonal */
/*  matrices and B is bidiagonal. */

/*  The test ratio to test the reduction is */
/*     RESID = norm( A - Q * B * PT ) / ( n * norm(A) * EPS ) */
/*  where PT = P' and EPS is the machine precision. */

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

/*  M       (input) INTEGER */
/*          The number of rows of the matrices A and Q. */

/*  N       (input) INTEGER */
/*          The number of columns of the matrices A and P'. */

/*  KD      (input) INTEGER */
/*          If KD = 0, B is diagonal and the array E is not referenced. */
/*          If KD = 1, the reduction was performed by xGEBRD; B is upper */
/*          bidiagonal if M >= N, and lower bidiagonal if M < N. */
/*          If KD = -1, the reduction was performed by xGBBRD; B is */
/*          always upper bidiagonal. */

/*  A       (input) REAL array, dimension (LDA,N) */
/*          The m by n matrix A. */

/*  LDA     (input) INTEGER */
/*          The leading dimension of the array A.  LDA >= max(1,M). */

/*  Q       (input) REAL array, dimension (LDQ,N) */
/*          The m by min(m,n) orthogonal matrix Q in the reduction */
/*          A = Q * B * P'. */

/*  LDQ     (input) INTEGER */
/*          The leading dimension of the array Q.  LDQ >= max(1,M). */

/*  D       (input) REAL array, dimension (min(M,N)) */
/*          The diagonal elements of the bidiagonal matrix B. */

/*  E       (input) REAL array, dimension (min(M,N)-1) */
/*          The superdiagonal elements of the bidiagonal matrix B if */
/*          m >= n, or the subdiagonal elements of B if m < n. */

/*  PT      (input) REAL array, dimension (LDPT,N) */
/*          The min(m,n) by n orthogonal matrix P' in the reduction */
/*          A = Q * B * P'. */

/*  LDPT    (input) INTEGER */
/*          The leading dimension of the array PT. */
/*          LDPT >= max(1,min(M,N)). */

/*  WORK    (workspace) REAL array, dimension (M+N) */

/*  RESID   (output) REAL */
/*          The test ratio:  norm(A - Q * B * P') / ( n * norm(A) * EPS ) */

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

/*     .. Parameters .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. External Functions .. */
/*     .. */
/*     .. External Subroutines .. */
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    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:
	    ;
	}
    }
    reslts[1] = errbnd;

//.........这里部分代码省略.........
			if (iju == 3 && ijvt == 3 || iju == 1 && ijvt == 1) {
			    goto L90;
			}
			*(unsigned char *)jobu = *(unsigned char *)&cjob[iju];
			*(unsigned char *)jobvt = *(unsigned char *)&cjob[
				ijvt];
			clacpy_("F", &m, &n, &asav[asav_offset], lda, &a[
				a_offset], lda);
			cgesvd_(jobu, jobvt, &m, &n, &a[a_offset], lda, &s[1],
				 &u[u_offset], ldu, &vt[vt_offset], ldvt, &
				work[1], &lswork, &rwork[1], &iinfo);

/*                    Compare U */

			dif = 0.f;
			if (m > 0 && n > 0) {
			    if (iju == 1) {
				cunt03_("C", &m, &mnmin, &m, &mnmin, &usav[
					usav_offset], ldu, &a[a_offset], lda, 
					&work[1], lwork, &rwork[1], &dif, &
					iinfo);
			    } else if (iju == 2) {
				cunt03_("C", &m, &mnmin, &m, &mnmin, &usav[
					usav_offset], ldu, &u[u_offset], ldu, 
					&work[1], lwork, &rwork[1], &dif, &
					iinfo);
			    } else if (iju == 3) {
				cunt03_("C", &m, &m, &m, &mnmin, &usav[
					usav_offset], ldu, &u[u_offset], ldu, 
					&work[1], lwork, &rwork[1], &dif, &
					iinfo);
			    }
			}
			result[4] = dmax(result[4],dif);

/*                    Compare VT */

			dif = 0.f;
			if (m > 0 && n > 0) {
			    if (ijvt == 1) {
				cunt03_("R", &n, &mnmin, &n, &mnmin, &vtsav[
					vtsav_offset], ldvt, &a[a_offset], 
					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;

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

		if (wantq) {
		    crot_(n, &q[(*n - *l + j) * q_dim1 + 1], &c__1, &q[(*n - *
			    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 */

//.........这里部分代码省略.........
					rwork[*nrhs + 1], &work[1], &rwork[(*
					nrhs << 1) + 1], &info);

/*                          Check the error code from CGBSVX. */

				if (info != izero) {
/* Writing concatenation */
				    i__11[0] = 1, a__1[0] = fact;
				    i__11[1] = 1, a__1[1] = trans;
				    s_cat(ch__1, a__1, i__11, &c__2, (ftnlen)
					    2);
				    alaerh_(path, "CGBSVX", &info, &izero, 
					    ch__1, &n, &n, &kl, &ku, nrhs, &
					    imat, &nfail, &nerrs, nout);
				}
/*                          Compare RWORK(2*NRHS+1) from CGBSVX with the */
/*                          computed reciprocal pivot growth RPVGRW */

				if (info != 0) {
				    anrmpv = 0.f;
				    i__8 = info;
				    for (j = 1; j <= i__8; ++j) {
/* Computing MAX */
					i__6 = ku + 2 - j;
/* Computing MIN */
					i__9 = n + ku + 1 - j, i__10 = kl + 
						ku + 1;
					i__7 = min(i__9,i__10);
					for (i__ = max(i__6,1); i__ <= i__7; 
						++i__) {
/* Computing MAX */
					    r__1 = anrmpv, r__2 = c_abs(&a[
						    i__ + (j - 1) * lda]);
					    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 << 

/* Subroutine */ int cgtt02_(char *trans, integer *n, integer *nrhs, complex *
	dl, complex *d__, complex *du, complex *x, integer *ldx, complex *b, 
	integer *ldb, real *rwork, real *resid)
{
    /* System generated locals */
    integer b_dim1, b_offset, x_dim1, x_offset, i__1;
    real r__1, r__2;

    /* Local variables */
    integer j;
    real eps;
    extern logical lsame_(char *, char *);
    real anorm, bnorm, xnorm;
    extern doublereal slamch_(char *), clangt_(char *, integer *, 
	    complex *, complex *, complex *);
    extern /* Subroutine */ int clagtm_(char *, integer *, integer *, real *, 
	    complex *, complex *, complex *, complex *, integer *, real *, 
	    complex *, integer *);
    extern doublereal scasum_(integer *, complex *, integer *);


/*  -- LAPACK test routine (version 3.1) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/*     November 2006 */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */

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

/*  CGTT02 computes the residual for the solution to a tridiagonal */
/*  system of equations: */
/*     RESID = norm(B - op(A)*X) / (norm(A) * norm(X) * EPS), */
/*  where EPS is the machine epsilon. */

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

/*  TRANS   (input) CHARACTER */
/*          Specifies the form of the residual. */
/*          = 'N':  B - A * X     (No transpose) */
/*          = 'T':  B - A**T * X  (Transpose) */
/*          = 'C':  B - A**H * X  (Conjugate transpose) */

/*  N       (input) INTEGTER */
/*          The order of the matrix A.  N >= 0. */

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

/*  DL      (input) COMPLEX array, dimension (N-1) */
/*          The (n-1) sub-diagonal elements of A. */

/*  D       (input) COMPLEX array, dimension (N) */
/*          The diagonal elements of A. */

/*  DU      (input) COMPLEX array, dimension (N-1) */
/*          The (n-1) super-diagonal elements of A. */

/*  X       (input) COMPLEX array, dimension (LDX,NRHS) */
/*          The computed solution vectors X. */

/*  LDX     (input) INTEGER */
/*          The leading dimension of the array X.  LDX >= max(1,N). */

/*  B       (input/output) COMPLEX array, dimension (LDB,NRHS) */
/*          On entry, the right hand side vectors for the system of */
/*          linear equations. */
/*          On exit, B is overwritten with the difference B - op(A)*X. */

/*  LDB     (input) INTEGER */
/*          The leading dimension of the array B.  LDB >= max(1,N). */

/*  RWORK   (workspace) REAL array, dimension (N) */

/*  RESID   (output) REAL */
/*          norm(B - op(A)*X) / (norm(A) * norm(X) * EPS) */

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

/*     .. Parameters .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. External Functions .. */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */

/*     Quick exit if N = 0 or NRHS = 0 */

    /* Parameter adjustments */
    --dl;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
		s = 0.f;
		i__5 = k + j * x_dim1;
		xk = (r__1 = X(k,j).r, dabs(r__1)) + (r__2 = r_imag(&X(k,j)), dabs(r__2));
		i__5 = k * ab_dim1 + 1;
		RWORK(k) += (r__1 = AB(1,k).r, dabs(r__1)) * xk;
		l = 1 - k;
/* Computing MIN */
		i__3 = *n, i__4 = k + *kd;
		i__5 = min(i__3,i__4);
		for (i = k + 1; i <= min(*n,k+*kd); ++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)));
/* L60: */
		}
		RWORK(k) += s;
/* L70: */
	    }
	}
	s = 0.f;
	i__2 = *n;
	for (i = 1; i <= *n; ++i) {
	    if (RWORK(i) > safe2) {
/* Computing MAX */
		i__5 = i;
		r__3 = s, r__4 = ((r__1 = WORK(i).r, dabs(r__1)) + (r__2 = 
			r_imag(&WORK(i)), dabs(r__2))) / RWORK(i);
		s = dmax(r__3,r__4);
	    } else {
/* Computing MAX */
		i__5 = i;
		r__3 = s, r__4 = ((r__1 = WORK(i).r, dabs(r__1)) + (r__2 = 
			r_imag(&WORK(i)), dabs(r__2)) + safe1) / (RWORK(i) + 
			safe1);
		s = dmax(r__3,r__4);
	    }
/* L80: */
	}
	BERR(j) = s;

/*        Test stopping criterion. Continue iterating if   
             1) The residual BERR(J) is larger than machine epsilon, a
nd   
             2) BERR(J) decreased by at least a factor of 2 during the
   
                last iteration, and   
             3) At most ITMAX iterations tried. */

	if (BERR(j) > eps && BERR(j) * 2.f <= lstres && count <= 5) {

/*           Update solution and try again. */

	    cpbtrs_(uplo, n, kd, &c__1, &AFB(1,1), ldafb, &WORK(1), n, 
		    info);
	    caxpy_(n, &c_b1, &WORK(1), &c__1, &X(1,j), &c__1);
	    lstres = BERR(j);
	    ++count;
	    goto L20;
	}

//.........这里部分代码省略.........
    /* Function Body */
    result[1] = 0.f;
    result[2] = 0.f;
    if (*n <= 0) {
	return 0;
    }

/* Computing MAX */
/* Computing MIN */
    i__3 = *n - 1;
    i__1 = 0, i__2 = min(i__3,*ka);
    ika = max(i__1,i__2);
    lw = *n * (*n + 1) / 2;

    if (lsame_(uplo, "U")) {
	lower = FALSE_;
	*(unsigned char *)cuplo = 'U';
    } else {
	lower = TRUE_;
	*(unsigned char *)cuplo = 'L';
    }

    unfl = slamch_("Safe minimum");
    ulp = slamch_("Epsilon") * slamch_("Base");

/*     Some Error Checks */

/*     Do Test 1 */

/*     Norm of A: */

/* Computing MAX */
    r__1 = slansb_("1", cuplo, n, &ika, &a[a_offset], lda, &work[1]);
    anorm = dmax(r__1,unfl);

/*     Compute error matrix:    Error = A - U S U' */

/*     Copy A from SB to SP storage format. */

    j = 0;
    i__1 = *n;
    for (jc = 1; jc <= i__1; ++jc) {
	if (lower) {
/* Computing MIN */
	    i__3 = ika + 1, i__4 = *n + 1 - jc;
	    i__2 = min(i__3,i__4);
	    for (jr = 1; jr <= i__2; ++jr) {
		++j;
		work[j] = a[jr + jc * a_dim1];
/* L10: */
	    }
	    i__2 = *n + 1 - jc;
	    for (jr = ika + 2; jr <= i__2; ++jr) {
		++j;
		work[j] = 0.f;
/* L20: */
	    }
	} else {
	    i__2 = jc;
	    for (jr = ika + 2; jr <= i__2; ++jr) {
		++j;
		work[j] = 0.f;
/* L30: */
	    }
/* Computing MIN */
	    i__2 = ika, i__3 = jc - 1;

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

	    cgemm_("C", "N", &n, &n, &n, &c_b2, &u[u_offset], ldu, &uz[
		    uz_offset], ldu, &c_b1, &z__[z_offset], ldu);
	    ntest = 8;

/*           Do Tests 3: | H - Z T Z' | / ( |H| n ulp ) */
/*                and 4: | I - Z Z' | / ( n ulp ) */

	    chst01_(&n, &ilo, &ihi, &h__[h_offset], lda, &t1[t1_offset], lda, 
		    &z__[z_offset], ldu, &work[1], nwork, &rwork[1], &result[
		    3]);

/*           Do Tests 5: | A - UZ T (UZ)' | / ( |A| n ulp ) */
/*                and 6: | I - UZ (UZ)' | / ( n ulp ) */

	    chst01_(&n, &ilo, &ihi, &a[a_offset], lda, &t1[t1_offset], lda, &
		    uz[uz_offset], ldu, &work[1], nwork, &rwork[1], &result[5]
);

/*           Do Test 7: | T2 - T1 | / ( |T| n ulp ) */

	    cget10_(&n, &n, &t2[t2_offset], lda, &t1[t1_offset], lda, &work[1]
, &rwork[1], &result[7]);

/*           Do Test 8: | W3 - W1 | / ( max(|W1|,|W3|) ulp ) */

	    temp1 = 0.f;
	    temp2 = 0.f;
	    i__3 = n;
	    for (j = 1; j <= i__3; ++j) {
/* Computing MAX */
		r__1 = temp1, r__2 = c_abs(&w1[j]), r__1 = max(r__1,r__2), 
			r__2 = c_abs(&w3[j]);
		temp1 = dmax(r__1,r__2);
/* Computing MAX */
		i__4 = j;
		i__5 = j;
		q__1.r = w1[i__4].r - w3[i__5].r, q__1.i = w1[i__4].i - w3[
			i__5].i;
		r__1 = temp2, r__2 = c_abs(&q__1);
		temp2 = dmax(r__1,r__2);
/* L130: */
	    }

/* Computing MAX */
	    r__1 = unfl, r__2 = ulp * dmax(temp1,temp2);
	    result[8] = temp2 / dmax(r__1,r__2);

/*           Compute the Left and Right Eigenvectors of T */

/*           Compute the Right eigenvector Matrix: */

	    ntest = 9;
	    result[9] = ulpinv;

/*           Select every other eigenvector */

	    i__3 = n;
	    for (j = 1; j <= i__3; ++j) {
		select[j] = FALSE_;
/* L140: */
	    }
	    i__3 = n;
	    for (j = 1; j <= i__3; j += 2) {
		select[j] = TRUE_;
/* L150: */

/* Subroutine */ int slasq3_(integer *i0, integer *n0, real *z__, integer *pp, 
	 real *dmin__, real *sigma, real *desig, real *qmax, integer *nfail, 
	integer *iter, integer *ndiv, logical *ieee, integer *ttype, real *
	dmin1, real *dmin2, real *dn, real *dn1, real *dn2, real *g, real *
	tau)
{
    /* System generated locals */
    integer i__1;
    real r__1, r__2;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    real s, t;
    integer j4, nn;
    real eps, tol;
    integer n0in, ipn4;
    real tol2, temp;
    extern /* Subroutine */ int slasq4_(integer *, integer *, real *, integer 
	    *, integer *, real *, real *, real *, real *, real *, real *, 
	    real *, integer *, real *), slasq5_(integer *, integer *, real *, 
	    integer *, real *, real *, real *, real *, real *, real *, real *, 
	     logical *), slasq6_(integer *, integer *, real *, integer *, 
	    real *, real *, real *, real *, real *, real *);
    extern doublereal slamch_(char *);
    extern logical sisnan_(real *);


/*  -- LAPACK routine (version 3.2)                                    -- */

/*  -- Contributed by Osni Marques of the Lawrence Berkeley National   -- */
/*  -- Laboratory and Beresford Parlett of the Univ. of California at  -- */
/*  -- Berkeley                                                        -- */
/*  -- November 2008                                                   -- */

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

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */

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

/*  SLASQ3 checks for deflation, computes a shift (TAU) and calls dqds. */
/*  In case of failure it changes shifts, and tries again until output */
/*  is positive. */

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

/*  I0     (input) INTEGER */
/*         First index. */

/*  N0     (input) INTEGER */
/*         Last index. */

/*  Z      (input) REAL array, dimension ( 4*N ) */
/*         Z holds the qd array. */

/*  PP     (input/output) INTEGER */
/*         PP=0 for ping, PP=1 for pong. */
/*         PP=2 indicates that flipping was applied to the Z array */
/*         and that the initial tests for deflation should not be */
/*         performed. */

/*  DMIN   (output) REAL */
/*         Minimum value of d. */

/*  SIGMA  (output) REAL */
/*         Sum of shifts used in current segment. */

/*  DESIG  (input/output) REAL */
/*         Lower order part of SIGMA */

/*  QMAX   (input) REAL */
/*         Maximum value of q. */

/*  NFAIL  (output) INTEGER */
/*         Number of times shift was too big. */

/*  ITER   (output) INTEGER */
/*         Number of iterations. */

/*  NDIV   (output) INTEGER */
/*         Number of divisions. */

/*  IEEE   (input) LOGICAL */
/*         Flag for IEEE or non IEEE arithmetic (passed to SLASQ5). */

/*  TTYPE  (input/output) INTEGER */
/*         Shift type. */

/*  DMIN1, DMIN2, DN, DN1, DN2, G, TAU (input/output) REAL */
/*         These are passed as arguments in order to save their values */
/*         between calls to SLASQ3. */

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

/* 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__ */

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

    /* Function Body */
    *resid = 0.f;
    if (*n <= 0) {
	return 0;
    }

/*     Compute B - U * S * V' one column at a time. */

    bnorm = 0.f;
    if (*kd >= 1) {

/*        B is bidiagonal. */

	if (lsame_(uplo, "U")) {

/*           B is upper bidiagonal. */

	    i__1 = *n;
	    for (j = 1; j <= i__1; ++j) {
		i__2 = *n;
		for (i__ = 1; i__ <= i__2; ++i__) {
		    work[*n + i__] = s[i__] * vt[i__ + j * vt_dim1];
/* L10: */
		}
		sgemv_("No transpose", n, n, &c_b6, &u[u_offset], ldu, &work[*
			n + 1], &c__1, &c_b8, &work[1], &c__1);
		work[j] += d__[j];
		if (j > 1) {
		    work[j - 1] += e[j - 1];
/* Computing MAX */
		    r__3 = bnorm, r__4 = (r__1 = d__[j], dabs(r__1)) + (r__2 =
			     e[j - 1], dabs(r__2));
		    bnorm = dmax(r__3,r__4);
		} else {
/* Computing MAX */
		    r__2 = bnorm, r__3 = (r__1 = d__[j], dabs(r__1));
		    bnorm = dmax(r__2,r__3);
		}
/* Computing MAX */
		r__1 = *resid, r__2 = sasum_(n, &work[1], &c__1);
		*resid = dmax(r__1,r__2);
/* L20: */
	    }
	} else {

/*           B is lower bidiagonal. */

	    i__1 = *n;
	    for (j = 1; j <= i__1; ++j) {
		i__2 = *n;
		for (i__ = 1; i__ <= i__2; ++i__) {
		    work[*n + i__] = s[i__] * vt[i__ + j * vt_dim1];
/* L30: */
		}
		sgemv_("No transpose", n, n, &c_b6, &u[u_offset], ldu, &work[*
			n + 1], &c__1, &c_b8, &work[1], &c__1);
		work[j] += d__[j];
		if (j < *n) {
		    work[j + 1] += e[j];
/* Computing MAX */
		    r__3 = bnorm, r__4 = (r__1 = d__[j], dabs(r__1)) + (r__2 =
			     e[j], dabs(r__2));
		    bnorm = dmax(r__3,r__4);
		} else {
/* Computing MAX */

/* the function below detects particles that have a number of neighbours 
 * outside the allowed tolerance range. For these, particles the smoothing
 * length is adjusted accordingly. Note that the smoothing length is
 */
void sidm_ensure_neighbours(int mode)
{
#define MAXITER 30

  int    i, ntot, last=0;
  float  *r2list;
  int    *ngblist, count, candidates;
  int    iter=0;
  double save;
  
  /*
  int    IndFirstUpdateBackup, NumForceUpdateBackup;
  
  IndFirstUpdateBackup= IndFirstUpdate;
  NumForceUpdateBackup= NumForceUpdate;
  for(i=IndFirstUpdate, count=0; count<NumForceUpdate; 
      i=P[i].ForceFlag, count++){
    P[i].ForceFlagBackup= P[i].ForceFlag;
  }
  */

  for(i=IndFirstUpdate, count=0, candidates=0; count<NumForceUpdate; 
      i=P[i].ForceFlag, count++) {
    if(P[i].Type>0) {
      if(P[i].NgbVelDisp < (All.DesNumNgb-All.MaxNumNgbDeviation) || 
	 (P[i].NgbVelDisp > (All.DesNumNgb+All.MaxNumNgbDeviation)))
        candidates++;
    }
  }

  MPI_Reduce(&candidates, &ntot, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Bcast(&ntot, 1, MPI_INT, 0, MPI_COMM_WORLD);
  
  if(ntot > 0) {
    //#ifdef FINDNBRLOG
      //      if(ThisTask==0)
      //{
      //printf("\n%d particles have too few/too many neighbours in sidm calculation!\n", ntot);
      // printf("Now fixing that...\n"); 
      // }
    //#endif
      
    for(i=IndFirstUpdate, count=0; count<NumForceUpdate; 
	i=P[i].ForceFlag, count++)   
      P[i].Left= P[i].Right= 0;
      
    do {
      for(i=1+N_gas, NumForceUpdate= 0, NumSphUpdate=0; i<=NumPart; i++) { 
	if( P[i].NgbVelDisp < (All.DesNumNgb-All.MaxNumNgbDeviation) || 
	    P[i].NgbVelDisp > (All.DesNumNgb+All.MaxNumNgbDeviation)) {
          if(P[i].Left>0 && P[i].Right>0)
            if((P[i].Right-P[i].Left) < 1.0e-3 * P[i].Left)
              continue;
	  
          if(NumForceUpdate==0)
            IndFirstUpdate= i;
          else
            P[last].ForceFlag= i;

          NumForceUpdate++;
          last=i;
          
          if(P[i].NgbVelDisp < (All.DesNumNgb-All.MaxNumNgbDeviation))
            P[i].Left= dmax(P[i].HsmlVelDisp, P[i].Left);
          else
            if(P[i].Right!=0)
              {
            if(P[i].HsmlVelDisp<P[i].Right)
              P[i].Right= P[i].HsmlVelDisp;
              }
            else
              P[i].Right= P[i].HsmlVelDisp;
        }
      }
      
      MPI_Allreduce(&NumForceUpdate, &ntot, 1, MPI_INT, MPI_SUM, 
		    MPI_COMM_WORLD);
      
      if(ntot > 0) {
	//#ifdef FINDNBRLOG
	//  if(ThisTask==0)
	//printf("ngb iteration %d.  still %d particles\n", iter, ntot);
	//#endif
	for(i=IndFirstUpdate, count=0; count<NumForceUpdate; 
	    i=P[i].ForceFlag, count++) {
          if(iter >= 20) {
	    printf("i=%d ID=%d Hsml=%g Left=%g Right=%g Ngbs=%d Right-Left=%g\n   pos=(%g|%g|%g)\n",
		   i, P[i].ID, P[i].HsmlVelDisp, P[i].Left, P[i].Right, 
		   P[i].NgbVelDisp, P[i].Right-P[i].Left,
		   P[i].PosPred[0], P[i].PosPred[1], P[i].PosPred[2]);
	  }
          
          if(iter == MAXITER) {
	    printf("ThisTask=%d Mi=(%g|%g|%g) Ma=(%g|%g|%g)\n", ThisTask,
		   DomainMin[P[i].Type][0], DomainMin[P[i].Type][1], 
		   DomainMin[P[i].Type][2], DomainMax[P[i].Type][0], 
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	    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);
	    xbnd = grow;
	    ip = jfirst * (jfirst + 1) / 2;

//.........这里部分代码省略.........
    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;

//.........这里部分代码省略.........
    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);
	}
/* L