kflops_util.c

Go to the documentation of this file.
00001 
00010 /* $Id: kflops_util.c,v 1.3 2004/09/20 15:24:36 yarkhan Exp $ */
00011 /* $UTK_Copyright: $ */
00012 
00013 /*----------------------*/ 
00014 #include "kflops.h"
00015 
00016 #ifdef __STDC__
00017 REAL kflops_epslon(REAL x)
00018 #else
00019 REAL kflops_epslon (x)
00020 REAL x;
00021 #endif
00022 /*
00023      estimate unit roundoff in quantities of size x.
00024 */
00025 
00026 {
00027   REAL a,b,c,eps;
00028 /*
00029      this program should function properly on all systems
00030      satisfying the following two assumptions,
00031         1.  the base used in representing dfloating point
00032             numbers is not a power of three.
00033         2.  the quantity  a  in statement 10 is represented to 
00034             the accuracy used in dfloating point variables
00035             that are stored in memory.
00036      the statement number 10 and the go to 10 are intended to
00037      force optimizing compilers to generate code satisfying 
00038      assumption 2.
00039      under these assumptions, it should be true that,
00040             a  is not exactly equal to four-thirds,
00041             b  has a zero for its last bit or digit,
00042             c  is not exactly equal to one,
00043             eps  measures the separation of 1.0 from
00044                  the next larger dfloating point number.
00045      the developers of eispack would appreciate being informed
00046      about any systems where these assumptions do not hold.
00047 
00048      *****************************************************************
00049      this routine is one of the auxiliary routines used by eispack iii
00050      to avoid machine dependencies.
00051      *****************************************************************
00052 
00053      this version dated 4/6/83.
00054 */
00055 
00056   a = 4.0e0/3.0e0;
00057   eps = ZERO;
00058   while (eps == ZERO) {
00059     b = a - ONE;
00060     c = b + b + b;
00061     eps = fabs((double)(c-ONE));
00062   }
00063   return(eps*fabs((double)x));
00064 }