ScaLAPACK  2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
pslaiect.c
Go to the documentation of this file.
00001 
00002 
00003 /* ---------------------------------------------------------------------
00004 *
00005 *  -- ScaLAPACK routine (version 1.5) --
00006 *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
00007 *     and University of California, Berkeley.
00008 *     May 1, 1997
00009 *
00010 *  ---------------------------------------------------------------------
00011 */
00012 
00013 /*
00014  * Include Files
00015  */
00016 #include "pxsyevx.h"
00017 #include <stdio.h>
00018 #include <math.h>
00019 #define  proto(x)       ()
00020 
00021 
00022 void pslasnbt_( int *ieflag )
00023 {
00024 /* 
00025 *
00026 *  Purpose
00027 *  ======= 
00028 *
00029 *  psalsnbt finds the position of the signbit of a single
00030 *  precision floating point number. This routine assumes IEEE
00031 *  arithmetic, and hence, tests only the 32nd bit as a possibility
00032 *  for the sign bit.
00033 *
00034 *  Note : For this release, we assume that sizeof(int) is 4 bytes.
00035 *
00036 *  Note : If a compile time flag (NO_IEEE) indicates that the
00037 *  machine does not have IEEE arithmetic, IEFLAG = 0 is returned.
00038 *
00039 *  Arguments
00040 *  =========
00041 *
00042 *  IEFLAG   (output) INTEGER
00043 *           This indicates the position of the signbit of any single
00044 *           precision floating point number.
00045 *           IEFLAG = 0 if the compile time flag, NO_IEEE, indicates
00046 *           that the machine does not have IEEE  arithmetic, or if
00047 *           sizeof(int) is different from 4 bytes.
00048 *           IEFLAG = 1 indicates that the sign bit is the 32nd bit.
00049 *
00050 *  =====================================================================
00051 *
00052 *  .. Local Scalars ..
00053 */
00054    float x;
00055    int         negone=-1, errornum;
00056    unsigned int *ix; 
00057 /* ..
00058 *  .. Executable Statements ..
00059 */
00060 
00061 #ifdef NO_IEEE
00062    *ieflag = 0;
00063 #else
00064    if(sizeof(int) != 4){
00065       *ieflag = 0;
00066       return;
00067    }
00068    x = (float) -1.0;
00069    ix = (unsigned int *) &x;
00070    if( *ix == 0xbff00000 )
00071    {
00072       *ieflag = 1;
00073    } else {
00074       *ieflag = 0; 
00075    }
00076 #endif
00077 }
00078 
00079 void pslaiect_( float *sigma, int *n, float *d, int *count )
00080 {
00081 /* 
00082 *
00083 *  Purpose
00084 *  ======= 
00085 *
00086 *  pslaiect computes the number of negative eigenvalues of (A- SIGMA I).
00087 *  This implementation of the Sturm Sequence loop exploits IEEE Arithmetic
00088 *  and has no conditionals in the innermost loop. The signbit is assumed
00089 *  to be bit 32.
00090 *
00091 *  Note that all arguments are call-by-reference so that this routine
00092 *  can be directly called from Fortran code.
00093 *
00094 *  This is a ScaLAPACK internal subroutine and arguments are not
00095 *  checked for unreasonable values.
00096 *
00097 *  Arguments
00098 *  =========
00099 *
00100 *  SIGMA    (input) REAL
00101 *           The shift. pslaiect finds the number of eigenvalues less
00102 *           than equal to SIGMA.
00103 *
00104 *  N        (input) INTEGER
00105 *           The order of the tridiagonal matrix T. N >= 1.
00106 *
00107 *  D        (input) REAL array, dimension (2*N - 1)
00108 *           Contains the diagonals and the squares of the off-diagonal
00109 *           elements of the tridiagonal matrix T. These elements are
00110 *           assumed to be interleaved in memory for better cache
00111 *           performance. The diagonal entries of T are in the entries
00112 *           D(1),D(3),...,D(2*N-1), while the squares of the off-diagonal
00113 *           entries are D(2),D(4),...,D(2*N-2). To avoid overflow, the
00114 *           matrix must be scaled so that its largest entry is no greater
00115 *           than overflow**(1/2) * underflow**(1/4) in absolute value,
00116 *           and for greatest accuracy, it should not be much smaller
00117 *           than that.
00118 *
00119 *  COUNT    (output) INTEGER
00120 *           The count of the number of eigenvalues of T less than or
00121 *           equal to SIGMA.
00122 *
00123 *  =====================================================================
00124 *
00125 *  .. Local Scalars ..
00126 */
00127    float       lsigma, tmp, *pd, *pe2;
00128    int         i;
00129 /* ..
00130 *  .. Executable Statements ..
00131 */
00132 
00133    lsigma = *sigma;
00134    pd = d; pe2 = d+1;
00135    tmp = *pd - lsigma; pd += 2;
00136    *count = (*((int *)&tmp) >> 31) & 1;
00137    for(i = 1;i < *n;i++){
00138       tmp = *pd - *pe2/tmp - lsigma;
00139       pd += 2; pe2 += 2;
00140       *count += ((*((int *)&tmp)) >> 31) & 1;
00141    }
00142 }
00143 
00144 void pslachkieee_( int *isieee, float *rmax, float *rmin )
00145 {
00146 /* 
00147 *
00148 *  Purpose
00149 *  ======= 
00150 *
00151 *  pslachkieee performs a simple check to make sure that the features
00152 *  of the IEEE standard that we rely on are implemented.  In some
00153 *  implementations, pslachkieee may not return.
00154 *
00155 *  Note that all arguments are call-by-reference so that this routine
00156 *  can be directly called from Fortran code.
00157 *
00158 *  This is a ScaLAPACK internal subroutine and arguments are not
00159 *  checked for unreasonable values.
00160 *
00161 *  Arguments
00162 *  =========
00163 *
00164 *  ISIEEE   (local output) INTEGER
00165 *           On exit, ISIEEE = 1 implies that all the features of the
00166 *           IEEE standard that we rely on are implemented.
00167 *           On exit, ISIEEE = 0 implies that some the features of the
00168 *           IEEE standard that we rely on are missing.
00169 *
00170 *  RMAX     (local input) REAL
00171 *           The overflow threshold ( = SLAMCH('O') ).
00172 *
00173 *  RMIN     (local input) REAL
00174 *           The underflow threshold ( = SLAMCH('U') ).
00175 *
00176 *  =====================================================================
00177 *
00178 *  .. Local Scalars ..
00179 */
00180    float x, pinf, pzero, ninf, nzero;
00181    int         ieflag, *ix, sbit1, sbit2, negone=-1, errornum;
00182 /* ..
00183 *  .. Executable Statements ..
00184 */
00185 
00186    pslasnbt_( &ieflag );
00187 
00188    pinf = *rmax / *rmin;
00189    pzero = 1.0 / pinf;
00190    pinf = 1.0 / pzero;
00191 
00192    if( pzero != 0.0 ){
00193       printf("pzero = %g should be zero\n",pzero);
00194       *isieee = 0; 
00195       return ;
00196    }
00197    if( ieflag == 1 ){
00198       sbit1 = (*((int *)&pzero) >> 31) & 1;
00199       sbit2 = (*((int *)&pinf) >> 31) & 1;
00200    }
00201    if( sbit1 == 1 ){
00202       printf("Sign of positive infinity is incorrect\n");
00203       *isieee = 0;
00204    }
00205    if( sbit2 == 1 ){
00206       printf("Sign of positive zero is incorrect\n");
00207       *isieee = 0;
00208    }
00209 
00210    ninf = -pinf;
00211    nzero = 1.0 / ninf;
00212    ninf = 1.0 / nzero;
00213 
00214    if( nzero != 0.0 ){
00215       printf("nzero = %g should be zero\n",nzero);
00216       *isieee = 0;
00217    }
00218    if( ieflag == 1 ){
00219       sbit1 = (*((int *)&nzero) >> 31) & 1;
00220       sbit2 = (*((int *)&ninf) >> 31) & 1;
00221    }
00222    if( sbit1 == 0 ){
00223       printf("Sign of negative infinity is incorrect\n");
00224       *isieee = 0;
00225    }
00226    if( sbit2 == 0 ){
00227       printf("Sign of negative zero is incorrect\n");
00228       *isieee = 0;
00229    }
00230 }