//.........这里部分代码省略.........
		"B");
	colequ = lsame_(equed, "C") || lsame_(equed, 
		"B");
	smlnum = slamch_("Safe minimum");
	bignum = 1.f / smlnum;
    }

/*     Test the input parameters. */

    if (! nofact && ! equil && ! lsame_(fact, "F")) {
	*info = -1;
    } else if (! notran && ! lsame_(trans, "T") && ! 
	    lsame_(trans, "C")) {
	*info = -2;
    } else if (*n < 0) {
	*info = -3;
    } else if (*nrhs < 0) {
	*info = -4;
    } else if (*lda < max(1,*n)) {
	*info = -6;
    } else if (*ldaf < max(1,*n)) {
	*info = -8;
    } else if (lsame_(fact, "F") && ! (rowequ || colequ 
	    || lsame_(equed, "N"))) {
	*info = -10;
    } else {
	if (rowequ) {
	    rcmin = bignum;
	    rcmax = 0.f;
	    i__1 = *n;
	    for (j = 1; j <= i__1; ++j) {
/* Computing MIN */
		r__1 = rcmin, r__2 = r__[j];
		rcmin = dmin(r__1,r__2);
/* Computing MAX */
		r__1 = rcmax, r__2 = r__[j];
		rcmax = dmax(r__1,r__2);
/* L10: */
	    }
	    if (rcmin <= 0.f) {
		*info = -11;
	    } else if (*n > 0) {
		rowcnd = dmax(rcmin,smlnum) / dmin(rcmax,bignum);
	    } else {
		rowcnd = 1.f;
	    }
	}
	if (colequ && *info == 0) {
	    rcmin = bignum;
	    rcmax = 0.f;
	    i__1 = *n;
	    for (j = 1; j <= i__1; ++j) {
/* Computing MIN */
		r__1 = rcmin, r__2 = c__[j];
		rcmin = dmin(r__1,r__2);
/* Computing MAX */
		r__1 = rcmax, r__2 = c__[j];
		rcmax = dmax(r__1,r__2);
/* L20: */
	    }
	    if (rcmin <= 0.f) {
		*info = -12;
	    } else if (*n > 0) {
		colcnd = dmax(rcmin,smlnum) / dmin(rcmax,bignum);
	    } else {
		colcnd = 1.f;

//.........这里部分代码省略.........
            The values computed by the two tests described above.  The   
            values are currently limited to 1/ulp, to avoid overflow.   
            RESULT(1) is always modified.  RESULT(2) is modified only   
            if LDU is at least N.   
            Modified.   

    =====================================================================   


       Parameter adjustments */
    a_dim1 = *lda;
    a_offset = 1 + a_dim1 * 1;
    a -= a_offset;
    --d__;
    --e;
    u_dim1 = *ldu;
    u_offset = 1 + u_dim1 * 1;
    u -= u_offset;
    v_dim1 = *ldv;
    v_offset = 1 + v_dim1 * 1;
    v -= v_offset;
    --tau;
    --work;
    --result;

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

    unfl = slamch_("Safe minimum");
    ulp = slamch_("Precision");

/*     Do Test 1   

       Norm of A:   

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

/*     Compute error matrix:   

       ITYPE=1: error = U' A U - S */

    ssymm_("L", uplo, n, m, &c_b6, &a[a_offset], lda, &u[u_offset], ldu, &
	    c_b7, &work[1], n);
    nn = *n * *n;
    nnp1 = nn + 1;
    sgemm_("T", "N", m, m, n, &c_b6, &u[u_offset], ldu, &work[1], n, &c_b7, &
	    work[nnp1], n);
    i__1 = *m;
    for (j = 1; j <= i__1; ++j) {
	jj = nn + (j - 1) * *n + j;
	work[jj] -= d__[j];
/* L10: */
    }
    if (*kband == 1 && *n > 1) {
	i__1 = *m;
	for (j = 2; j <= i__1; ++j) {
	    jj1 = nn + (j - 1) * *n + j - 1;
	    jj2 = nn + (j - 2) * *n + j;
	    work[jj1] -= e[j - 1];
	    work[jj2] -= e[j - 1];
/* L20: */
	}
    }
    wnorm = slansy_("1", uplo, m, &work[nnp1], n, &work[1]);

    if (anorm > wnorm) {
	result[1] = wnorm / anorm / (*m * ulp);
    } else {
	if (anorm < 1.f) {
/* Computing MIN */
	    r__1 = wnorm, r__2 = *m * anorm;
	    result[1] = dmin(r__1,r__2) / anorm / (*m * ulp);
	} else {
/* Computing MIN */
	    r__1 = wnorm / anorm, r__2 = (real) (*m);
	    result[1] = dmin(r__1,r__2) / (*m * ulp);
	}
    }

/*     Do Test 2   

       Compute  U'U - I */

    if (*itype == 1) {
	i__1 = (*n << 1) * *n;
	sort01_("Columns", n, m, &u[u_offset], ldu, &work[1], &i__1, &result[
		2]);
    }

    return 0;

/*     End of SSYT22 */

} /* ssyt22_ */

//.........这里部分代码省略.........
/*     subdiagonal element has become negligible. */

    for (its = 0; its <= 30; ++its) {

/*        Look for a single small subdiagonal element. */

	i__1 = l + 1;
	for (k = i__; k >= i__1; --k) {
	    if ((r__1 = h__[k + (k - 1) * h_dim1], dabs(r__1)) <= smlnum) {
		goto L40;
	    }
	    tst = (r__1 = h__[k - 1 + (k - 1) * h_dim1], dabs(r__1)) + (r__2 =
		     h__[k + k * h_dim1], dabs(r__2));
	    if (tst == 0.f) {
		if (k - 2 >= *ilo) {
		    tst += (r__1 = h__[k - 1 + (k - 2) * h_dim1], dabs(r__1));
		}
		if (k + 1 <= *ihi) {
		    tst += (r__1 = h__[k + 1 + k * h_dim1], dabs(r__1));
		}
	    }
/*           ==== The following is a conservative small subdiagonal */
/*           .    deflation  criterion due to Ahues & Tisseur (LAWN 122, */
/*           .    1997). It has better mathematical foundation and */
/*           .    improves accuracy in some cases.  ==== */
	    if ((r__1 = h__[k + (k - 1) * h_dim1], dabs(r__1)) <= ulp * tst) {
/* Computing MAX */
		r__3 = (r__1 = h__[k + (k - 1) * h_dim1], dabs(r__1)), r__4 = 
			(r__2 = h__[k - 1 + k * h_dim1], dabs(r__2));
		ab = dmax(r__3,r__4);
/* Computing MIN */
		r__3 = (r__1 = h__[k + (k - 1) * h_dim1], dabs(r__1)), r__4 = 
			(r__2 = h__[k - 1 + k * h_dim1], dabs(r__2));
		ba = dmin(r__3,r__4);
/* Computing MAX */
		r__3 = (r__1 = h__[k + k * h_dim1], dabs(r__1)), r__4 = (r__2 
			= h__[k - 1 + (k - 1) * h_dim1] - h__[k + k * h_dim1],
			 dabs(r__2));
		aa = dmax(r__3,r__4);
/* Computing MIN */
		r__3 = (r__1 = h__[k + k * h_dim1], dabs(r__1)), r__4 = (r__2 
			= h__[k - 1 + (k - 1) * h_dim1] - h__[k + k * h_dim1],
			 dabs(r__2));
		bb = dmin(r__3,r__4);
		s = aa + ab;
/* Computing MAX */
		r__1 = smlnum, r__2 = ulp * (bb * (aa / s));
		if (ba * (ab / s) <= dmax(r__1,r__2)) {
		    goto L40;
		}
	    }
	}
L40:
	l = k;
	if (l > *ilo) {

/*           H(L,L-1) is negligible */

	    h__[l + (l - 1) * h_dim1] = 0.f;
	}

/*        Exit from loop if a submatrix of order 1 or 2 has split off. */

	if (l >= i__ - 1) {
	    goto L150;
	}

/** Calculate cell grid dimensions, cell sizes and number of cells.
 *  Calculates the cell grid, based on \ref local_box_l and \ref
 *  max_range. If the number of cells is larger than \ref
 *  max_num_cells, it increases max_range until the number of cells is
 *  smaller or equal \ref max_num_cells. It sets: \ref
 *  DomainDecomposition::cell_grid, \ref
 *  DomainDecomposition::ghost_cell_grid, \ref
 *  DomainDecomposition::cell_size, \ref
 *  DomainDecomposition::inv_cell_size, and \ref n_cells.
 */
void dd_create_cell_grid()
{
  int i,n_local_cells,new_cells,min_ind;
  double cell_range[3], min_size, scale, volume;
  CELL_TRACE(fprintf(stderr, "%d: dd_create_cell_grid: max_range %f\n",this_node,max_range));
  CELL_TRACE(fprintf(stderr, "%d: dd_create_cell_grid: local_box %f-%f, %f-%f, %f-%f,\n",this_node,my_left[0],my_right[0],my_left[1],my_right[1],my_left[2],my_right[2]));
  
  /* initialize */
  cell_range[0]=cell_range[1]=cell_range[2] = max_range;

  if (max_range < ROUND_ERROR_PREC*box_l[0]) {
    /* this is the initialization case */
    n_local_cells = dd.cell_grid[0] = dd.cell_grid[1] = dd.cell_grid[2]=1;
  }
  else {
    /* Calculate initial cell grid */
    volume = local_box_l[0];
    for(i=1;i<3;i++) volume *= local_box_l[i];
    scale = pow(max_num_cells/volume, 1./3.);
    for(i=0;i<3;i++) {
      /* this is at least 1 */
      dd.cell_grid[i] = (int)ceil(local_box_l[i]*scale);
      cell_range[i] = local_box_l[i]/dd.cell_grid[i];

      if ( cell_range[i] < max_range ) {
	/* ok, too many cells for this direction, set to minimum */
	dd.cell_grid[i] = (int)floor(local_box_l[i]/max_range);
	if ( dd.cell_grid[i] < 1 ) {
	  char *error_msg = runtime_error(ES_INTEGER_SPACE + 2*ES_DOUBLE_SPACE + 128);
	  ERROR_SPRINTF(error_msg, "{002 interaction range %g in direction %d is larger than the local box size %g} ",
			max_range, i, local_box_l[i]);
	  dd.cell_grid[i] = 1;
	}
	cell_range[i] = local_box_l[i]/dd.cell_grid[i];
      }
    }

    /* It may be necessary to asymmetrically assign the scaling to the coordinates, which the above approach will not do.
       For a symmetric box, it gives a symmetric result. Here we correct that. */
    for (;;) {
      n_local_cells = dd.cell_grid[0];
      for (i = 1; i < 3; i++)
	n_local_cells *= dd.cell_grid[i];

      /* done */
      if (n_local_cells <= max_num_cells)
	break;

      /* find coordinate with the smallest cell range */
      min_ind = 0;
      min_size = cell_range[0];
      for (i = 1; i < 3; i++)
	if (dd.cell_grid[i] > 1 && cell_range[i] < min_size) {
	  min_ind = i;
	  min_size = cell_range[i];
	}
      CELL_TRACE(fprintf(stderr, "%d: minimal coordinate %d, size %f, grid %d\n", this_node,min_ind, min_size, dd.cell_grid[min_ind]));

      dd.cell_grid[min_ind]--;
      cell_range[min_ind] = local_box_l[min_ind]/dd.cell_grid[min_ind];
    }
    CELL_TRACE(fprintf(stderr, "%d: final %d %d %d\n", this_node, dd.cell_grid[0], dd.cell_grid[1], dd.cell_grid[2]));

    /* sanity check */
    if (n_local_cells < min_num_cells) {
      char *error_msg = runtime_error(ES_INTEGER_SPACE + 2*ES_DOUBLE_SPACE + 128);
      ERROR_SPRINTF(error_msg, "{001 number of cells %d is smaller than minimum %d (interaction range too large or min_num_cells too large)} ",
		    n_local_cells, min_num_cells);
    }
  }

  /* quit program if unsuccesful */
  if(n_local_cells > max_num_cells) {
    char *error_msg = runtime_error(128);
    ERROR_SPRINTF(error_msg, "{003 no suitable cell grid found} ");
  }

  /* now set all dependent variables */
  new_cells=1;
  for(i=0;i<3;i++) {
    dd.ghost_cell_grid[i] = dd.cell_grid[i]+2;	
    new_cells              *= dd.ghost_cell_grid[i];
    dd.cell_size[i]       = local_box_l[i]/(double)dd.cell_grid[i];
    dd.inv_cell_size[i]   = 1.0 / dd.cell_size[i];
  }
  max_skin = dmin(dmin(dd.cell_size[0],dd.cell_size[1]),dd.cell_size[2]) - max_cut;

  /* allocate cell array and cell pointer arrays */
  realloc_cells(new_cells);
  realloc_cellplist(&local_cells, local_cells.n = n_local_cells);
//.........这里部分代码省略.........

/* Subroutine */ int slasq4_(integer *i0, integer *n0, real *z__, integer *pp, 
	 integer *n0in, real *dmin__, real *dmin1, real *dmin2, real *dn, 
	real *dn1, real *dn2, real *tau, integer *ttype, real *g)
{
    /* System generated locals */
    integer i__1;
    real r__1, r__2;

    /* Local variables */
    real s, a2, b1, b2;
    integer i4, nn, np;
    real gam, gap1, gap2;

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

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

/*  SLASQ4 computes an approximation TAU to the smallest eigenvalue */
/*  using values of d from the previous transform. */

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

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

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

/*  PP    (input) INTEGER */
/*        PP=0 for ping, PP=1 for pong. */

/*  NOIN  (input) INTEGER */
/*        The value of N0 at start of EIGTEST. */

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

/*  DMIN1 (input) REAL */
/*        Minimum value of d, excluding D( N0 ). */

/*  DMIN2 (input) REAL */
/*        Minimum value of d, excluding D( N0 ) and D( N0-1 ). */

/*  DN    (input) REAL */
/*        d(N) */

/*  DN1   (input) REAL */
/*        d(N-1) */

/*  DN2   (input) REAL */
/*        d(N-2) */

/*  TAU   (output) REAL */
/*        This is the shift. */

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

/*  G     (input/output) REAL */
/*        G is passed as an argument in order to save its value between */
/*        calls to SLASQ4. */

/*  Further Details */
/*  =============== */
/*  CNST1 = 9/16 */

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

/*     A negative DMIN forces the shift to take that absolute value */
/*     TTYPE records the type of shift. */

    /* Parameter adjustments */
    --z__;

    /* Function Body */
    if (*dmin__ <= 0.f) {
	*tau = -(*dmin__);
	*ttype = -1;
	return 0;
    }

    nn = (*n0 << 2) + *pp;
    if (*n0in == *n0) {

/*        No eigenvalues deflated. */

	if (*dmin__ == *dn || *dmin__ == *dn1) {

	    b1 = sqrt(z__[nn - 3]) * sqrt(z__[nn - 5]);
	    b2 = sqrt(z__[nn - 7]) * sqrt(z__[nn - 9]);
	    a2 = z__[nn - 7] + z__[nn - 5];
//.........这里部分代码省略.........

/* Subroutine */ int cbdt02_(integer *m, integer *n, complex *b, integer *ldb, 
	 complex *c__, integer *ldc, complex *u, integer *ldu, complex *work, 
	real *rwork, real *resid)
{
    /* System generated locals */
    integer b_dim1, b_offset, c_dim1, c_offset, u_dim1, u_offset, i__1;
    real r__1, r__2;

    /* Local variables */
    integer j;
    real eps;
    extern /* Subroutine */ int cgemv_(char *, integer *, integer *, complex *
, complex *, integer *, complex *, integer *, complex *, complex *
, integer *);
    real bnorm;
    extern /* Subroutine */ int ccopy_(integer *, complex *, integer *, 
	    complex *, integer *);
    extern doublereal clange_(char *, integer *, integer *, complex *, 
	    integer *, real *), slamch_(char *);
    real realmn;
    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 */
/*  ======= */

/*  CBDT02 tests the change of basis C = U' * B by computing the residual */

/*     RESID = norm( B - U * C ) / ( max(m,n) * norm(B) * EPS ), */

/*  where B and C are M by N matrices, U is an M by M orthogonal matrix, */
/*  and EPS is the machine precision. */

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

/*  M       (input) INTEGER */
/*          The number of rows of the matrices B and C and the order of */
/*          the matrix Q. */

/*  N       (input) INTEGER */
/*          The number of columns of the matrices B and C. */

/*  B       (input) COMPLEX array, dimension (LDB,N) */
/*          The m by n matrix B. */

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

/*  C       (input) COMPLEX array, dimension (LDC,N) */
/*          The m by n matrix C, assumed to contain U' * B. */

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

/*  U       (input) COMPLEX array, dimension (LDU,M) */
/*          The m by m orthogonal matrix U. */

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

/*  WORK    (workspace) COMPLEX array, dimension (M) */

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

/*  RESID   (output) REAL */
/*          RESID = norm( B - U * C ) / ( max(m,n) * norm(B) * EPS ), */

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

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

/*     Quick return if possible */

    /* Parameter adjustments */
    b_dim1 = *ldb;
    b_offset = 1 + b_dim1;
    b -= b_offset;
    c_dim1 = *ldc;
    c_offset = 1 + c_dim1;
    c__ -= c_offset;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	    ++(*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 (irange == 3) {

/*        RANGE='I': Compute the interval containing eigenvalues */
/*                   IL through IU. */

/*        Compute Gershgorin interval for entire (split) matrix */
/*        and use it as the initial interval */

	gu = d__[1];
	gl = d__[1];
	tmp1 = 0.f;

	i__1 = *n - 1;
	for (j = 1; j <= i__1; ++j) {
	    tmp2 = sqrt(work[j]);
/* Computing MAX */
	    r__1 = gu, r__2 = d__[j] + tmp1 + tmp2;
	    gu = dmax(r__1,r__2);
/* Computing MIN */
	    r__1 = gl, r__2 = d__[j] - tmp1 - tmp2;
	    gl = dmin(r__1,r__2);
	    tmp1 = tmp2;
/* L20: */
	}

/* Computing MAX */
	r__1 = gu, r__2 = d__[*n] + tmp1;
	gu = dmax(r__1,r__2);
/* Computing MIN */
	r__1 = gl, r__2 = d__[*n] - tmp1;
	gl = dmin(r__1,r__2);
/* Computing MAX */
	r__1 = dabs(gl), r__2 = dabs(gu);
	tnorm = dmax(r__1,r__2);
	gl = gl - tnorm * 2.1f * ulp * *n - pivmin * 4.2000000000000002f;
	gu = gu + tnorm * 2.1f * ulp * *n + pivmin * 2.1f;

/*        Compute Iteration parameters */

	itmax = (integer) ((log(tnorm + pivmin) - log(pivmin)) / log(2.f)) + 
		2;
	if (*abstol <= 0.f) {
	    atoli = ulp * tnorm;
	} else {
	    atoli = *abstol;
	}

	work[*n + 1] = gl;
	work[*n + 2] = gl;
	work[*n + 3] = gu;
	work[*n + 4] = gu;
	work[*n + 5] = gl;
	work[*n + 6] = gu;

//.........这里部分代码省略.........
    INFO    (output) INTEGER   
            = 0:  successful exit   
            < 0:  if INFO = -i, the i-th argument had an illegal value   
            > 0:  if INFO = i, the i-th diagonal element is nonpositive.   

    =====================================================================   


       Test the input parameters.   

       Parameter adjustments */
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2;
    real r__1, r__2;
    /* Builtin functions */
    double sqrt(doublereal);
    /* Local variables */
    static real smin;
    static integer i__;
    extern /* Subroutine */ int xerbla_(char *, integer *);
#define a_subscr(a_1,a_2) (a_2)*a_dim1 + a_1
#define a_ref(a_1,a_2) a[a_subscr(a_1,a_2)]

    a_dim1 = *lda;
    a_offset = 1 + a_dim1 * 1;
    a -= a_offset;
    --s;

    /* Function Body */
    *info = 0;
    if (*n < 0) {
	*info = -1;
    } else if (*lda < max(1,*n)) {
	*info = -3;
    }
    if (*info != 0) {
	i__1 = -(*info);
	xerbla_("CPOEQU", &i__1);
	return 0;
    }

/*     Quick return if possible */

    if (*n == 0) {
	*scond = 1.f;
	*amax = 0.f;
	return 0;
    }

/*     Find the minimum and maximum diagonal elements. */

    i__1 = a_subscr(1, 1);
    s[1] = a[i__1].r;
    smin = s[1];
    *amax = s[1];
    i__1 = *n;
    for (i__ = 2; i__ <= i__1; ++i__) {
	i__2 = a_subscr(i__, i__);
	s[i__] = a[i__2].r;
/* Computing MIN */
	r__1 = smin, r__2 = s[i__];
	smin = dmin(r__1,r__2);
/* Computing MAX */
	r__1 = *amax, r__2 = s[i__];
	*amax = dmax(r__1,r__2);
/* L10: */
    }

    if (smin <= 0.f) {

/*        Find the first non-positive diagonal element and return. */

	i__1 = *n;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    if (s[i__] <= 0.f) {
		*info = i__;
		return 0;
	    }
/* L20: */
	}
    } else {

/*        Set the scale factors to the reciprocals   
          of the diagonal elements. */

	i__1 = *n;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    s[i__] = 1.f / sqrt(s[i__]);
/* L30: */
	}

/*        Compute SCOND = min(S(I)) / max(S(I)) */

	*scond = sqrt(smin) / sqrt(*amax);
    }
    return 0;

/*     End of CPOEQU */

} /* cpoequ_ */

//.........这里部分代码省略.........
#ifndef  TWODIMS
		      hinv4 = hinv * hinv * hinv * hinv;
#else
		      hinv4 = hinv * hinv * hinv / boxSize_Z;
#endif
		      u = r * hinv;
		      if(u < 0.5)
			dwk_j = hinv4 * u * (KERNEL_COEFF_3 * u - KERNEL_COEFF_4);
		      else
			dwk_j = hinv4 * KERNEL_COEFF_6 * (1.0 - u) * (1.0 - u);
		    }
		  else
		    {
		      dwk_j = 0;
		    }

		  if(soundspeed_i + soundspeed_j > maxSignalVel)
		    maxSignalVel = soundspeed_i + soundspeed_j;

		  if(vdotr2 < 0)	/* ... artificial viscosity */
		    {
#ifndef MONAGHAN83VISC
		      mu_ij = fac_mu * vdotr2 / r;	/* note: this is negative! */
#else		      
		      h_ij = 0.5 * (h_i + h_j);
		      mu_ij = fac_mu * h_ij * vdotr2 / (r2 + 0.0001 * h_ij * h_ij);
#endif
                      vsig = soundspeed_i + soundspeed_j - 3 * mu_ij;



		      if(vsig > maxSignalVel)
			maxSignalVel = vsig;

		      rho_ij = 0.5 * (rho + SphP[j].Density);

#ifdef MORRIS97VISC		      
                      visc = 0.25 * (alpha_visc + alpha_visc_j) * vsig * (-mu_ij) / rho_ij;
#else

		      f2 =
			fabs(SphP[j].DivVel) / (fabs(SphP[j].DivVel) + SphP[j].CurlVel +
						0.0001 * soundspeed_j / fac_mu / SphP[j].Hsml);

#ifndef MONAGHAN83VISC
		      visc = 0.25 * All.ArtBulkViscConst * vsig * (-mu_ij) / rho_ij * (f1 + f2);
#else			      
		      soundspeed_ij = (soundspeed_i+soundspeed_j) * 0.5;

		      visc = ((-All.ArtBulkViscConst) * soundspeed_ij * mu_ij + All.ArtBulkViscBeta * mu_ij * mu_ij) / rho_ij;
#endif //MONAGHAN83VISC

#endif //MORRIS97VISC

		      /* .... end artificial viscosity evaluation */
#ifndef NOVISCOSITYLIMITER
		      /* make sure that viscous acceleration is not too large */
		      dt = imax(timestep, (P[j].Ti_endstep - P[j].Ti_begstep)) * All.Timebase_interval;
		      if(dt > 0 && (dwk_i + dwk_j) < 0)
			{
			  visc = dmin(visc, 0.5 * fac_vsic_fix * vdotr2 /
				      (0.5 * (mass + P[j].Mass) * (dwk_i + dwk_j) * r * dt));
			}
#endif
		    }
		  else
		    visc = 0;

		  p_over_rho2_j *= SphP[j].DhsmlDensityFactor;

		  hfc_visc = 0.5 * P[j].Mass * visc * (dwk_i + dwk_j) / r;

		  hfc = hfc_visc + P[j].Mass * (p_over_rho2_i * dwk_i + p_over_rho2_j * dwk_j) / r;

		  acc[0] -= hfc * dx;
		  acc[1] -= hfc * dy;
		  acc[2] -= hfc * dz;
		  dtEntropy += 0.5 * hfc_visc * vdotr2;
		}
	    }
	}
    }
  while(startnode >= 0);

  /* Now collect the result at the right place */
  if(mode == 0)
    {
      for(k = 0; k < 3; k++)
	SphP[target].HydroAccel[k] = acc[k];
      SphP[target].DtEntropy = dtEntropy;
      SphP[target].MaxSignalVel = maxSignalVel;
    }
  else
    {
      for(k = 0; k < 3; k++)
	HydroDataResult[target].Acc[k] = acc[k];
      HydroDataResult[target].DtEntropy = dtEntropy;
      HydroDataResult[target].MaxSignalVel = maxSignalVel;
    }
}

//.........这里部分代码省略.........
    r__1 = r1mach_(&c__4);
    tol = dmax(r__1,1e-18f);
    fid = (real) (*id);
    if (az > 1.f) {
	goto L60;
    }
/* ----------------------------------------------------------------------- */
/*     POWER SERIES FOR ABS(Z).LE.1. */
/* ----------------------------------------------------------------------- */
    s1.r = cone.r, s1.i = cone.i;
    s2.r = cone.r, s2.i = cone.i;
    if (az < tol) {
	goto L160;
    }
    aa = az * az;
    if (aa < tol / az) {
	goto L40;
    }
    trm1.r = cone.r, trm1.i = cone.i;
    trm2.r = cone.r, trm2.i = cone.i;
    atrm = 1.f;
    q__2.r = z__->r * z__->r - z__->i * z__->i, q__2.i = z__->r * z__->i + 
	    z__->i * z__->r;
    q__1.r = q__2.r * z__->r - q__2.i * z__->i, q__1.i = q__2.r * z__->i + 
	    q__2.i * z__->r;
    z3.r = q__1.r, z3.i = q__1.i;
    az3 = az * aa;
    ak = fid + 2.f;
    bk = 3.f - fid - fid;
    ck = 4.f - fid;
    dk = fid + 3.f + fid;
    d1 = ak * dk;
    d2 = bk * ck;
    ad = dmin(d1,d2);
    ak = fid * 9.f + 24.f;
    bk = 30.f - fid * 9.f;
    z3r = z3.r;
    z3i = r_imag(&z3);
    for (k = 1; k <= 25; ++k) {
	r__1 = z3r / d1;
	r__2 = z3i / d1;
	q__2.r = r__1, q__2.i = r__2;
	q__1.r = trm1.r * q__2.r - trm1.i * q__2.i, q__1.i = trm1.r * q__2.i 
		+ trm1.i * q__2.r;
	trm1.r = q__1.r, trm1.i = q__1.i;
	q__1.r = s1.r + trm1.r, q__1.i = s1.i + trm1.i;
	s1.r = q__1.r, s1.i = q__1.i;
	r__1 = z3r / d2;
	r__2 = z3i / d2;
	q__2.r = r__1, q__2.i = r__2;
	q__1.r = trm2.r * q__2.r - trm2.i * q__2.i, q__1.i = trm2.r * q__2.i 
		+ trm2.i * q__2.r;
	trm2.r = q__1.r, trm2.i = q__1.i;
	q__1.r = s2.r + trm2.r, q__1.i = s2.i + trm2.i;
	s2.r = q__1.r, s2.i = q__1.i;
	atrm = atrm * az3 / ad;
	d1 += ak;
	d2 += bk;
	ad = dmin(d1,d2);
	if (atrm < tol * ad) {
	    goto L40;
	}
	ak += 18.f;
	bk += 18.f;
/* L30: */
    }

//.........这里部分代码省略.........
/*                 or Large Orders, NPL Mathematical Tables 6, Her */
/*                 Majesty's Stationery Office, London, 1962. */
/* ***ROUTINES CALLED  ALNGAM, ASYJY, I1MACH, JAIRY, R1MACH, XERMSG */
/* ***REVISION HISTORY  (YYMMDD) */
/*   750101  DATE WRITTEN */
/*   890531  Changed all specific intrinsics to generic.  (WRB) */
/*   890531  REVISION DATE from Version 3.2 */
/*   891214  Prologue converted to Version 4.0 format.  (BAB) */
/*   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ) */
/*   900326  Removed duplicate information from DESCRIPTION section. */
/*           (WRB) */
/*   920501  Reformatted the REFERENCES section.  (WRB) */
/* ***END PROLOGUE  BESJ */
    /* Parameter adjustments */
    --y;

    /* Function Body */
/* ***FIRST EXECUTABLE STATEMENT  BESJ */
    *nz = 0;
    kt = 1;
    ns = 0;
/*     I1MACH(14) REPLACES I1MACH(11) IN A DOUBLE PRECISION CODE */
/*     I1MACH(15) REPLACES I1MACH(12) IN A DOUBLE PRECISION CODE */
    ta = r1mach_(&c__3);
    tol = dmax(ta,1e-15f);
    i1 = i1mach_(&c__11) + 1;
    i2 = i1mach_(&c__12);
    tb = r1mach_(&c__5);
    elim1 = (i2 * tb + 3.f) * -2.303f;
    rtol = 1.f / tol;
    slim = r1mach_(&c__1) * 1e3f * rtol;
/*     TOLLN = -LN(TOL) */
    tolln = tb * 2.303f * i1;
    tolln = dmin(tolln,34.5388f);
    if ((i__1 = *n - 1) < 0) {
	goto L720;
    } else if (i__1 == 0) {
	goto L10;
    } else {
	goto L20;
    }
L10:
    kt = 2;
L20:
    nn = *n;
    if (*x < 0.f) {
	goto L730;
    } else if (*x == 0) {
	goto L30;
    } else {
	goto L80;
    }
L30:
    if (*alpha < 0.f) {
	goto L710;
    } else if (*alpha == 0) {
	goto L40;
    } else {
	goto L50;
    }
L40:
    y[1] = 1.f;
    if (*n == 1) {
	return 0;
    }
    i1 = 2;

//.........这里部分代码省略.........
	    ++(*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 (irange == 3) {

/*        RANGE='I': Compute the interval containing eigenvalues   
                     IL through IU.   

          Compute Gershgorin interval for entire (split) matrix   
          and use it as the initial interval */

	gu = D(1);
	gl = D(1);
	tmp1 = 0.f;

	i__1 = *n - 1;
	for (j = 1; j <= *n-1; ++j) {
	    tmp2 = sqrt(WORK(j));
/* Computing MAX */
	    r__1 = gu, r__2 = D(j) + tmp1 + tmp2;
	    gu = dmax(r__1,r__2);
/* Computing MIN */
	    r__1 = gl, r__2 = D(j) - tmp1 - tmp2;
	    gl = dmin(r__1,r__2);
	    tmp1 = tmp2;
/* L20: */
	}

/* Computing MAX */
	r__1 = gu, r__2 = D(*n) + tmp1;
	gu = dmax(r__1,r__2);
/* Computing MIN */
	r__1 = gl, r__2 = D(*n) - tmp1;
	gl = dmin(r__1,r__2);
/* Computing MAX */
	r__1 = dabs(gl), r__2 = dabs(gu);
	tnorm = dmax(r__1,r__2);
	gl = gl - tnorm * 2.f * ulp * *n - pivmin * 4.f;
	gu = gu + tnorm * 2.f * ulp * *n + pivmin * 2.f;

/*        Compute Iteration parameters */

	itmax = (integer) ((log(tnorm + pivmin) - log(pivmin)) / log(2.f)) + 
		2;
	if (*abstol <= 0.f) {
	    atoli = ulp * tnorm;
	} else {
	    atoli = *abstol;
	}

	WORK(*n + 1) = gl;
	WORK(*n + 2) = gl;
	WORK(*n + 3) = gu;
	WORK(*n + 4) = gu;
	WORK(*n + 5) = gl;
	WORK(*n + 6) = gu;

//.........这里部分代码省略.........
	dpsi = 0.f;
	psi = 0.f;
	erretm = 0.f;
	i__1 = ii;
	for (j = 1; j <= i__1; ++j) {
	    temp = z__[j] / delta[j];
	    psi += z__[j] * temp;
	    dpsi += temp * temp;
	    erretm += psi;
/* L40: */
	}
	erretm = dabs(erretm);

/*        Evaluate PHI and the derivative DPHI */

	temp = z__[*n] / delta[*n];
	phi = z__[*n] * temp;
	dphi = temp * temp;
	erretm = (-phi - psi) * 8.f + erretm - phi + rhoinv + dabs(tau) * (
		dpsi + dphi);

	w = rhoinv + phi + psi;

/*        Test for convergence */

	if (dabs(w) <= eps * erretm) {
	    *dlam = d__[*i__] + tau;
	    goto L250;
	}

	if (w <= 0.f) {
	    dltlb = dmax(dltlb,tau);
	} else {
	    dltub = dmin(dltub,tau);
	}

/*        Calculate the new step */

	++niter;
	c__ = w - delta[*n - 1] * dpsi - delta[*n] * dphi;
	a = (delta[*n - 1] + delta[*n]) * w - delta[*n - 1] * delta[*n] * (
		dpsi + dphi);
	b = delta[*n - 1] * delta[*n] * w;
	if (c__ < 0.f) {
	    c__ = dabs(c__);
	}
	if (c__ == 0.f) {
/*          ETA = B/A   
             ETA = RHO - TAU */
	    eta = dltub - tau;
	} else if (a >= 0.f) {
	    eta = (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) / (
		    c__ * 2.f);
	} else {
	    eta = b * 2.f / (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(
		    r__1))));
	}

/*        Note, eta should be positive if w is negative, and   
          eta should be negative otherwise. However,   
          if for some reason caused by roundoff, eta*w > 0,   
          we simply use one Newton step instead. This way   
          will guarantee eta*w < 0. */

	if (w * eta > 0.f) {
	    eta = -w / (dpsi + dphi);

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

/*     Quick return if possible */

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

    if (*n == 1) {
	if (alleig || indeig) {
	    *m = 1;
	    w[1] = a[a_dim1 + 1];
	} else {
	    if (*vl < a[a_dim1 + 1] && *vu >= a[a_dim1 + 1]) {
		*m = 1;
		w[1] = a[a_dim1 + 1];
	    }
	}
	if (wantz) {
	    z__[z_dim1 + 1] = 1.f;
	}
	return 0;
    }

/*     Get machine constants. */

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

    iscale = 0;
    abstll = *abstol;
    if (valeig) {
	vll = *vl;
	vuu = *vu;
    }
    anrm = slansy_("M", uplo, n, &a[a_offset], lda, &work[1]);
    if (anrm > 0.f && anrm < rmin) {
	iscale = 1;
	sigma = rmin / anrm;
    } else if (anrm > rmax) {
	iscale = 1;
	sigma = rmax / anrm;
    }
    if (iscale == 1) {
	if (lower) {
	    i__1 = *n;
	    for (j = 1; j <= i__1; ++j) {
		i__2 = *n - j + 1;
		sscal_(&i__2, &sigma, &a[j + j * a_dim1], &c__1);
	    }
	} else {
	    i__1 = *n;
	    for (j = 1; j <= i__1; ++j) {
		sscal_(&j, &sigma, &a[j * a_dim1 + 1], &c__1);
	    }
	}
	if (*abstol > 0.f) {
	    abstll = *abstol * sigma;

//.........这里部分代码省略.........
			- dfdy[i__4].i;
		dfdy[i__1].r = q__1.r, dfdy[i__1].i = q__1.i;
	    }
	    numer = scnrm2_(n, &dfdy[dfdy_offset], matdim);
	    if (el[*nq * 13 + 1] * numer <= *uround * 100.f * y0nrm) {
		if (*rmax == 2.f) {
		    switch__ = TRUE_;
		    goto L170;
		}
	    }
	    i__1 = *n;
	    for (i__ = 1; i__ <= i__1; ++i__) {
/* L136: */
		i__3 = i__ * dfdy_dim1 + 1;
		i__4 = i__;
		dfdy[i__3].r = save1[i__4].r, dfdy[i__3].i = save1[i__4].i;
	    }
	    if (denom != 0.f) {
/* Computing MAX */
		r__1 = bnd, r__2 = numer / (denom * dabs(*h__) * el[*nq * 13 
			+ 1]);
		bnd = dmax(r__1,r__2);
	    }
	}
    }
    if (iter > 0) {
/* Computing MAX */
	r__1 = *trend * .9f, r__2 = d__ / d1;
	*trend = dmax(r__1,r__2);
    }
    d1 = d__;
/* Computing MIN */
    r__1 = *trend * 2.f;
    ctest = dmin(r__1,1.f) * d__;
    if (ctest <= *eps) {
	goto L170;
    }
    ++iter;
    if (iter < 3) {
	i__3 = *n;
	for (i__ = 1; i__ <= i__3; ++i__) {
/* L140: */
	    i__4 = i__;
	    i__1 = i__ + yh_dim1;
	    i__2 = *nq * 13 + 1;
	    i__5 = i__;
	    q__2.r = el[i__2] * save1[i__5].r, q__2.i = el[i__2] * save1[i__5]
		    .i;
	    q__1.r = yh[i__1].r + q__2.r, q__1.i = yh[i__1].i + q__2.i;
	    y[i__4].r = q__1.r, y[i__4].i = q__1.i;
	}
	(*f)(n, t, &y[1], &save2[1]);
	if (*n == 0) {
	    *jstate = 6;
	    goto L430;
	}
	++(*nfe);
	goto L130;
    }
/*                     The corrector iteration failed to converge in */
/*                     MXITER tries.  If partials are involved but are */
/*                     not up to date, they are reevaluated for the next */
/*                     try.  Otherwise the YH array is retracted to its */
/*                     values before prediction, and H is reduced, if */
/*                     possible.  If not, a no-convergence exit is taken. */
    if (*convrg) {

//.........这里部分代码省略.........
	    jlast = *n;
	    jinc = 1;
	    maind = 1;
	}

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

	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 = 1.f / dmax(xbnd,smlnum);
	    xbnd = grow;
	    i__1 = jlast;
	    i__2 = jinc;
	    for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {

/*              Exit the loop if the growth factor is too small. */

		if (grow <= smlnum) {
		    goto L50;
		}

/*              M(j) = G(j-1) / abs(A(j,j)) */

		tjj = (r__1 = ab[maind + j * ab_dim1], dabs(r__1));
/* Computing MIN */
		r__1 = xbnd, r__2 = dmin(1.f,tjj) * grow;
		xbnd = dmin(r__1,r__2);
		if (tjj + cnorm[j] >= smlnum) {

/*                 G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) ) */

		    grow *= tjj / (tjj + cnorm[j]);
		} else {

/*                 G(j) could overflow, set GROW to 0. */

		    grow = 0.f;
		}
/* L30: */
	    }
	    grow = xbnd;
	} else {

/*           A is unit triangular. */

/*           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. */

/* Computing MIN */
	    r__1 = 1.f, r__2 = 1.f / dmax(xbnd,smlnum);
	    grow = dmin(r__1,r__2);
	    i__2 = jlast;
	    i__1 = jinc;
	    for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {

/*              Exit the loop if the growth factor is too small. */

		if (grow <= smlnum) {
		    goto L50;

/* Subroutine */ int CPraxis::min_(C_INT *n, C_INT *j, C_INT *nits, C_FLOAT64 *
                                   d2, C_FLOAT64 *x1, C_FLOAT64 *f1, bool *fk, FPraxis *f, C_FLOAT64 *
                                   x, C_FLOAT64 *t, C_FLOAT64 *machep, C_FLOAT64 *h__)
{
  /* System generated locals */
  C_INT i__1;
  C_FLOAT64 d__1, d__2;

  /* Local variables */
  static C_FLOAT64 temp;
  static C_INT i__, k;
  static C_FLOAT64 s, small, d1, f0, f2, m2, m4, t2, x2, fm;
  static bool dz;
  static C_FLOAT64 xm, sf1, sx1;

  /* ...THE SUBROUTINE MIN MINIMIZES F FROM X IN THE DIRECTION V(*,J) UNLESS */
  /*   J IS LESS THAN 1, WHEN A QUADRATIC SEARCH IS MADE IN THE PLANE */
  /*   DEFINED BY Q0,Q1,X. */
  /*   D2 IS EITHER ZERO OR AN APPROXIMATION TO HALF F". */
  /*   ON ENTRY, X1 IS AN ESTIMATE OF THE DISTANCE FROM X TO THE MINIMUM */
  /*   ALONG V(*,J) (OR, IF J=0, A CURVE).  ON RETURN, X1 IS THE DISTANCE */
  /*   FOUND. */
  /*   IF FK=.TRUE., THEN F1 IS FLIN(X1).  OTHERWISE X1 AND F1 ARE IGNORED */
  /*   ON ENTRY UNLESS FINAL FX IS GREATER THAN F1. */
  /*   NITS CONTROLS THE NUMBER OF TIMES AN ATTEMPT WILL BE MADE TO HALVE */
  /*   THE INTERVAL. */
  /* Parameter adjustments */
  --x;

  /* Function Body */
  /* Computing 2nd power */
  d__1 = *machep;
  small = d__1 * d__1;
  m2 = sqrt(*machep);
  m4 = sqrt(m2);
  sf1 = *f1;
  sx1 = *x1;
  k = 0;
  xm = 0.;
  fm = global_1.fx;
  f0 = global_1.fx;
  dz = *d2 < *machep;
  /* ...FIND THE STEP SIZE... */
  s = 0.;
  i__1 = *n;

  for (i__ = 1; i__ <= i__1; ++i__)
    {
      /* L1: */
      /* Computing 2nd power */
      d__1 = x[i__];
      s += d__1 * d__1;
    }

  s = sqrt(s);
  temp = *d2;

  if (dz)
    {
      temp = global_1.dmin__;
    }

  t2 = m4 * sqrt(fabs(global_1.fx) / temp + s * global_1.ldt) + m2 *
       global_1.ldt;
  s = m4 * s + *t;

  if (dz && t2 > s)
    {
      t2 = s;
    }

  t2 = dmax(t2, small);
  /* Computing MIN */
  d__1 = t2, d__2 = *h__ * .01;
  t2 = dmin(d__1, d__2);

  if (!(*fk) || *f1 > fm)
    {
      goto L2;
    }

  xm = *x1;
  fm = *f1;
L2:

  if (*fk && fabs(*x1) >= t2)
    {
      goto L3;
    }

  temp = 1.;

  if (*x1 < 0.)
    {
      temp = -1.;
    }

  *x1 = temp * t2;
  *f1 = flin_(n, j, x1, f, &x[1], &global_1.nf);
L3:
//.........这里部分代码省略.........

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

    unfl = slamch_("Safe minimum");
    ulp = slamch_("Precision");

/*     Do Test 1   

       Copy A & Compute its 1-Norm: */

    slaset_("Full", n, n, &c_b5, &c_b5, &work[1], n);

    anorm = 0.f;
    temp1 = 0.f;

    i__1 = *n - 1;
    for (j = 1; j <= i__1; ++j) {
	work[(*n + 1) * (j - 1) + 1] = ad[j];
	work[(*n + 1) * (j - 1) + 2] = ae[j];
	temp2 = (r__1 = ae[j], dabs(r__1));
/* Computing MAX */
	r__2 = anorm, r__3 = (r__1 = ad[j], dabs(r__1)) + temp1 + temp2;
	anorm = dmax(r__2,r__3);
	temp1 = temp2;
/* L10: */
    }

/* Computing 2nd power */
    i__1 = *n;
    work[i__1 * i__1] = ad[*n];
/* Computing MAX */
    r__2 = anorm, r__3 = (r__1 = ad[*n], dabs(r__1)) + temp1, r__2 = max(r__2,
	    r__3);
    anorm = dmax(r__2,unfl);

/*     Norm of A - USU' */

    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
	r__1 = -sd[j];
	ssyr_("L", n, &r__1, &u_ref(1, j), &c__1, &work[1], n);
/* L20: */
    }

    if (*n > 1 && *kband == 1) {
	i__1 = *n - 1;
	for (j = 1; j <= i__1; ++j) {
	    r__1 = -se[j];
	    ssyr2_("L", n, &r__1, &u_ref(1, j), &c__1, &u_ref(1, j + 1), &
		    c__1, &work[1], n);
/* L30: */
	}
    }

/* Computing 2nd power */
    i__1 = *n;
    wnorm = slansy_("1", "L", n, &work[1], n, &work[i__1 * i__1 + 1]);

    if (anorm > wnorm) {
	result[1] = wnorm / anorm / (*n * ulp);
    } else {
	if (anorm < 1.f) {
/* Computing MIN */
	    r__1 = wnorm, r__2 = *n * anorm;
	    result[1] = dmin(r__1,r__2) / anorm / (*n * ulp);
	} else {
/* Computing MIN */
	    r__1 = wnorm / anorm, r__2 = (real) (*n);
	    result[1] = dmin(r__1,r__2) / (*n * ulp);
	}
    }

/*     Do Test 2   

       Compute  UU' - I */

    sgemm_("N", "C", n, n, n, &c_b19, &u[u_offset], ldu, &u[u_offset], ldu, &
	    c_b5, &work[1], n);

    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
	work[(*n + 1) * (j - 1) + 1] += -1.f;
/* L40: */
    }

/* Computing MIN   
   Computing 2nd power */
    i__1 = *n;
    r__1 = (real) (*n), r__2 = slange_("1", n, n, &work[1], n, &work[i__1 * 
	    i__1 + 1]);
    result[2] = dmin(r__1,r__2) / (*n * ulp);

    return 0;

/*     End of SSTT21 */

} /* sstt21_ */