pdlaiectl_()

void pdlaiectl_	(	double *	sigma,
		Int *	n,
		double *	d,
		Int *	count
	)

Definition at line 151 of file pdlaiect.c.

{
/* 
*
*  Purpose
*  ======= 
*
*  pdlaiectl computes the number of negative eigenvalues of (A- SIGMA I).
*  This implementation of the Sturm Sequence loop exploits IEEE Arithmetic
*  and has no conditionals in the innermost loop. To extract the signbit,
*  this routine assumes that the double precision word is stored in
*  "Little Endian" word order, i.e, the signbit is assumed to be bit 64.
*
*  Note that all arguments are call-by-reference so that this routine
*  can be directly called from Fortran code.
*
*  This is a ScaLAPACK internal subroutine and arguments are not
*  checked for unreasonable values.
*
*  Arguments
*  =========
*
*  SIGMA    (input) DOUBLE PRECISION
*           The shift. pdlaiectl finds the number of eigenvalues
*           less than equal to SIGMA.
*
*  N        (input) INTEGER
*           The order of the tridiagonal matrix T. N >= 1.
*
*  D        (input) DOUBLE PRECISION array, dimension (2*N - 1)
*           Contains the diagonals and the squares of the off-diagonal
*           elements of the tridiagonal matrix T. These elements are
*           assumed to be interleaved in memory for better cache
*           performance. The diagonal entries of T are in the entries
*           D(1),D(3),...,D(2*N-1), while the squares of the off-diagonal
*           entries are D(2),D(4),...,D(2*N-2). To avoid overflow, the
*           matrix must be scaled so that its largest entry is no greater
*           than overflow**(1/2) * underflow**(1/4) in absolute value,
*           and for greatest accuracy, it should not be much smaller
*           than that.
*
*  COUNT    (output) INTEGER
*           The count of the number of eigenvalues of T less than or
*           equal to SIGMA.
*
*  =====================================================================
*
*  .. Local Scalars ..
*/
   double      lsigma, tmp, *pd, *pe2;
   Int         i;
/* ..
*  .. Executable Statements ..
*/
 
   lsigma = *sigma;
   pd = d; pe2 = d+1;
   tmp = *pd - lsigma; pd += 2;
   *count = (*(((Int *)&tmp)+1) >> 31) & 1;
   for(i = 1;i < *n;i++){
      tmp = *pd - *pe2/tmp - lsigma;
      pd += 2; pe2 += 2;
      *count += (*(((Int *)&tmp)+1) >> 31) & 1;
   }
}