/* --------------------------------------------------------------------- * * -- ScaLAPACK routine (version 1.5) -- * University of Tennessee, Knoxville, Oak Ridge National Laboratory, * and University of California, Berkeley. * May 1, 1997 * * --------------------------------------------------------------------- */ /* * Include Files */ #include "pxsyevx.h" #include #include #define proto(x) () void pdlasnbt_( ieflag ) /* * .. Scalar Arguments .. */ int *ieflag; { /* * * Purpose * ======= * * pdalsnbt finds the position of the signbit of a double * double precision floating point number. This routine assumes IEEE * arithmetic, and hence, tests only the 32nd and 64th bits as * possibilities for the sign bit. * * Note : For this release, we assume that sizeof(int) is 4 bytes. * * Note : If a compile time flag (NO_IEEE) indicates that the * machine does not have IEEE arithmetic, IEFLAG = 0 is returned. * * Arguments * ========= * * IEFLAG (output) INTEGER * This indicates the position of the signbit of any double * precision floating point number. * IEFLAG = 0 if the compile time flag, NO_IEEE, indicates * that the machine does not have IEEE arithmetic, or if * sizeof(int) is different from 4 bytes. * IEFLAG = 1 indicates that the sign bit is the 32nd * bit ( Big Endian ). * IEFLAG = 2 indicates that the sign bit is the 64th * bit ( Little Endian ). * * ===================================================================== * * .. Local Scalars .. */ double x; int negone=-1, errornum; unsigned int *ix; /* .. * .. Executable Statements .. */ #ifdef NO_IEEE *ieflag = 0; #else if(sizeof(int) != 4){ *ieflag = 0; return; } x = (double) -1.0; ix = (unsigned int *) &x; if(( *ix == 0xbff00000) && ( *(ix+1) == 0x0) ) { *ieflag = 1; } else if(( *(ix+1) == 0xbff00000) && ( *ix == 0x0) ) { *ieflag = 2; } else { *ieflag = 0; } #endif } void pdlaiectb_( sigma, n, d, count ) /* * .. Scalar Arguments .. */ double *sigma, *d; int *n, *count; { /* * * Purpose * ======= * * pdlaiectb 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 * "Big Endian" word order, i.e, the signbit is assumed to be bit 32. * * 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. pdlaiectb 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) >> 31) & 1; for(i = 1;i < *n;i++){ tmp = *pd - *pe2/tmp - lsigma; pd += 2; pe2 += 2; *count += ((*((int *)&tmp)) >> 31) & 1; } } void pdlaiectl_( sigma, n, d, count ) /* * .. Scalar Arguments .. */ double *sigma, *d; int *n, *count; { /* * * 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; } } pdlachkieee_( isieee, rmax, rmin ) /* * .. Scalar Arguments .. */ double *rmax, *rmin; int *isieee; { /* * * Purpose * ======= * * pdlachkieee performs a simple check to make sure that the features * of the IEEE standard that we rely on are implemented. In some * implementations, pdlachkieee may not return. * * 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 * ========= * * ISIEEE (local output) INTEGER * On exit, ISIEEE = 1 implies that all the features of the * IEEE standard that we rely on are implemented. * On exit, ISIEEE = 0 implies that some the features of the * IEEE standard that we rely on are missing. * * RMAX (local input) DOUBLE PRECISION * The overflow threshold ( = DLAMCH('O') ). * * RMIN (local input) DOUBLE PRECISION * The underflow threshold ( = DLAMCH('U') ). * * ===================================================================== * * .. Local Scalars .. */ double x, pinf, pzero, ninf, nzero; int ieflag, *ix, sbit1, sbit2, negone=-1, errornum; /* .. * .. Executable Statements .. */ pdlasnbt_( &ieflag ); pinf = *rmax / *rmin; pzero = 1.0 / pinf; pinf = 1.0 / pzero; if( pzero != 0.0 ){ printf("pzero = %g should be zero\n",pzero); *isieee = 0; return ; } if( ieflag == 1 ){ sbit1 = (*((int *)&pzero) >> 31) & 1; sbit2 = (*((int *)&pinf) >> 31) & 1; }else if(ieflag == 2){ sbit1 = (*(((int *)&pzero)+1) >> 31) & 1; sbit2 = (*(((int *)&pinf)+1) >> 31) & 1; } if( sbit1 == 1 ){ printf("Sign of positive infinity is incorrect\n"); *isieee = 0; } if( sbit2 == 1 ){ printf("Sign of positive zero is incorrect\n"); *isieee = 0; } ninf = -pinf; nzero = 1.0 / ninf; ninf = 1.0 / nzero; if( nzero != 0.0 ){ printf("nzero = %g should be zero\n",nzero); *isieee = 0; } if( ieflag == 1 ){ sbit1 = (*((int *)&nzero) >> 31) & 1; sbit2 = (*((int *)&ninf) >> 31) & 1; }else if(ieflag == 2){ sbit1 = (*(((int *)&nzero)+1) >> 31) & 1; sbit2 = (*(((int *)&ninf)+1) >> 31) & 1; } if( sbit1 == 0 ){ printf("Sign of negative infinity is incorrect\n"); *isieee = 0; } if( sbit2 == 0 ){ printf("Sign of negative zero is incorrect\n"); *isieee = 0; } }