Commit 54bc0e27 authored by miller86's avatar miller86

FPE fixups

git-svn-id: http://visit.ilight.com/svn/visit/trunk/src@24627 18c085ea-50e0-402c-830e-de6fd14e8384
parent 822dd8f1
......@@ -163,7 +163,9 @@ avtMagnitudeExpression::DeriveVariable(vtkDataSet *in_ds, int currentDomainsInde
dtype *r = (dtype*)results->GetVoidPointer(0); \
for (int i = 0, idx = 0 ; i < ntuples ; i++, idx += 3) \
{ \
r[i] = sqrt(x[idx+0]*x[idx+0]+x[idx+1]*x[idx+1]+x[idx+2]*x[idx+2]); \
r[i] = sqrt((double)x[idx+0]*(double)x[idx+0]+\
(double)x[idx+1]*(double)x[idx+1]+\
(double)x[idx+2]*(double)x[idx+2]); \
} \
}
......
......@@ -130,6 +130,8 @@
// Hank Childs, Fri Jan 9 14:10:25 PST 2009
// Initialize jitter.
//
// Mark C. Miller, Thu Oct 2 09:41:37 PDT 2014
// Initialize lightDirection.
// ****************************************************************************
avtSamplePointExtractor::avtSamplePointExtractor(int w, int h, int d)
......@@ -175,6 +177,7 @@ avtSamplePointExtractor::avtSamplePointExtractor(int w, int h, int d)
rayCastingSLIVR = false;
lighting = false;
lightPosition[0] = lightPosition[1] = lightPosition[2] = 0.0; lightPosition[3] = 1.0;
lightDirection[0] = 0; lightDirection[1] = 0; lightDirection[2] = -1;
materialProperties[0] = 0.4; materialProperties[1] = 0.75; materialProperties[3] = 0.0; materialProperties[3] = 15.0;
}
......
......@@ -59,8 +59,8 @@ public:
union {
struct {
float lm;
float rm;
double lm;
double rm;
};
struct {
......@@ -69,7 +69,7 @@ public:
};
};
void make_node( unsigned int left, unsigned int d, float b[2] )
void make_node( unsigned int left, unsigned int d, double b[2] )
{
index = (d & 3) | (left << 2);
lm = b[0];
......@@ -101,12 +101,12 @@ public:
return index & 3;
}
const float& lmax() const
const double& lmax() const
{
return lm;
}
const float& rmin() const
const double& rmin() const
{
return rm;
}
......@@ -145,9 +145,9 @@ public:
const celltree& m_ct;
unsigned int m_stack[32];
unsigned int* m_sp;
const float* m_pos;
const double* m_pos;
point_traversal( const celltree& ct, const float* pos ) :
point_traversal( const celltree& ct, const double* pos ) :
m_ct(ct), m_pos(pos)
{
m_stack[0] = 0;
......@@ -166,7 +166,7 @@ public:
if( n->is_leaf() )
return n;
const float p = m_pos[n->dim()];
const double p = m_pos[n->dim()];
const unsigned int left = n->left();
bool l = p <= n->lmax();
......@@ -202,18 +202,18 @@ private:
struct bucket
{
float min;
float max;
double min;
double max;
unsigned int cnt;
bucket()
{
cnt = 0;
min = std::numeric_limits<float>::max();
max = -std::numeric_limits<float>::max();
min = std::numeric_limits<double>::max();
max = -std::numeric_limits<double>::max();
}
void add( const float _min, const float _max )
void add( const double _min, const double _max )
{
++cnt;
......@@ -227,8 +227,8 @@ private:
struct per_cell
{
float min[3];
float max[3];
double min[3];
double max[3];
unsigned int ind;
};
......@@ -250,9 +250,9 @@ private:
struct left_predicate
{
unsigned int d;
float p;
double p;
left_predicate( unsigned int _d, float _p ) :
left_predicate( unsigned int _d, double _p ) :
d(_d), p(2.0f*_p)
{
}
......@@ -267,7 +267,7 @@ private:
// -------------------------------------------------------------------------
void find_min_max( const per_cell* begin, const per_cell* end,
float* min, float* max )
double* min, double* max )
{
if( begin == end )
return;
......@@ -291,7 +291,7 @@ private:
// -------------------------------------------------------------------------
void find_min_d( const per_cell* begin, const per_cell* end,
unsigned int d, float& min )
unsigned int d, double& min )
{
min = begin->min[d];
......@@ -301,7 +301,7 @@ private:
}
void find_max_d( const per_cell* begin, const per_cell* end,
unsigned int d, float& max )
unsigned int d, double& max )
{
max = begin->max[d];
......@@ -312,7 +312,7 @@ private:
// -------------------------------------------------------------------------
void split( unsigned int index, float min[3], float max[3] )
void split( unsigned int index, double min[3], double max[3] )
{
unsigned int start = m_nodes[index].start();
unsigned int size = m_nodes[index].size();
......@@ -326,16 +326,17 @@ private:
const int nbuckets = 6;
const float ext[3] = { max[0]-min[0], max[1]-min[1], max[2]-min[2] };
const float iext[3] = { nbuckets/ext[0], nbuckets/ext[1], nbuckets/ext[2] };
const double ext[3] = { max[0]-min[0], max[1]-min[1], max[2]-min[2] };
const double iext[3] = { ext[0]!=0 ? nbuckets/ext[0] : -1,
ext[1]!=0 ? nbuckets/ext[1] : -1,
ext[2]!=0 ? nbuckets/ext[2] : -1};
bucket b[3][nbuckets];
for( const per_cell* pc=begin; pc!=end; ++pc )
{
for( unsigned int d=0; d<3; ++d )
{
float cen = (pc->min[d] + pc->max[d])/2.0f;
double cen = (pc->min[d] + pc->max[d])/2.0f;
int ind = (int)( (cen-min[d])*iext[d] );
if( ind<0 )
......@@ -348,8 +349,8 @@ private:
}
}
float cost = std::numeric_limits<float>::max();
float plane = 0;
double cost = std::numeric_limits<double>::max();
double plane = 0;
unsigned int dim = 0;
for( unsigned int d=0; d<3; ++d )
......@@ -358,8 +359,8 @@ private:
for( unsigned int n=0; n<nbuckets-1; ++n )
{
float lmax = -std::numeric_limits<float>::max();
float rmin = std::numeric_limits<float>::max();
double lmax = -std::numeric_limits<double>::max();
double rmin = std::numeric_limits<double>::max();
for( unsigned int m=0; m<=n; ++m )
if( b[d][m].max > lmax )
......@@ -371,21 +372,27 @@ private:
sum += b[d][n].cnt;
float lvol = (lmax-min[d])/ext[d];
float rvol = (max[d]-rmin)/ext[d];
if (ext[d] != 0 &&
lmax != -std::numeric_limits<double>::max() &&
rmin != std::numeric_limits<double>::max())
{
double lvol = (lmax-min[d])/ext[d];
double rvol = (max[d]-rmin)/ext[d];
float c = lvol*sum + rvol*(size-sum);
double c = lvol*sum + rvol*(size-sum);
if( sum > 0 && sum < size && c < cost )
{
cost = c;
dim = d;
plane = min[d] + (n+1)/iext[d];
if( sum > 0 && sum < size && c < cost )
{
cost = c;
dim = d;
if (iext[d] > 0)
plane = min[d] + (n+1)/iext[d];
}
}
}
}
if( cost != std::numeric_limits<float>::max() )
if( cost != std::numeric_limits<double>::max() )
mid = std::partition( begin, end, left_predicate( dim, plane ) );
// fallback
......@@ -397,12 +404,12 @@ private:
std::nth_element( begin, mid, end, center_order( dim ) );
}
float lmin[3], lmax[3], rmin[3], rmax[3];
double lmin[3], lmax[3], rmin[3], rmax[3];
find_min_max( begin, mid, lmin, lmax );
find_min_max( mid, end, rmin, rmax );
float clip[2] = { lmax[dim], rmin[dim] };
double clip[2] = { lmax[dim], rmin[dim] };
celltree::node child[2];
child[0].make_leaf( begin - m_pc, mid-begin );
......@@ -430,16 +437,16 @@ public:
m_pc = new per_cell[size];
float min[3] = {
std::numeric_limits<float>::max(),
std::numeric_limits<float>::max(),
std::numeric_limits<float>::max()
double min[3] = {
std::numeric_limits<double>::max(),
std::numeric_limits<double>::max(),
std::numeric_limits<double>::max()
};
float max[3] = {
-std::numeric_limits<float>::max(),
-std::numeric_limits<float>::max(),
-std::numeric_limits<float>::max(),
double max[3] = {
-std::numeric_limits<double>::max(),
-std::numeric_limits<double>::max(),
-std::numeric_limits<double>::max(),
};
for( unsigned int i=0; i<size; ++i )
......@@ -608,7 +615,7 @@ void avtCellLocatorBIH::FindCellRecursive( const double pos[3],
else
{
// else descend the tree
const float p = pos[n.dim()];
const double p = pos[n.dim()];
const unsigned int left = n.left();
bool l = p <= n.lmax();
......
......@@ -662,19 +662,28 @@ avtICAlgorithm::ComputeStatistic(ICStatistics &stats)
nVals++;
}
}
stats.mean = stats.total / (float)nVals;
if (nVals != 0)
stats.mean = stats.total / (float)nVals;
else
stats.mean = stats.value;
float sum = 0.0;
double sum = 0.0;
for (int i = 0; i < nProcs; i++)
{
if (output[i] >= 0.0)
{
float x = output[i] - stats.mean;
double x = output[i] - stats.mean;
sum += (x*x);
}
}
sum /= (float)nVals;
stats.sigma = sqrt(sum);
if (nVals != 0)
{
sum /= (float)nVals;
if (sum > 0)
stats.sigma = sqrt(sum);
else
stats.sigma = 0.0;
}
stats.histogram.resize(nVals);
int i, j;
......@@ -1002,6 +1011,8 @@ avtICAlgorithm::ReportCounters(ostream &os, bool totals)
// Dave Pugmire, Thu Feb 12 08:47:29 EST 2009
// Better formatting for stats output.
//
// Mark C. Miller, Thu Oct 2 09:23:56 PDT 2014
// Defend against FPE div by zero
// ****************************************************************************
void
avtICAlgorithm::PrintTiming(ostream &os,
......@@ -1018,7 +1029,10 @@ avtICAlgorithm::PrintTiming(ostream &os,
if (total)
{
os<<s.total;
os<<" ["<<100.0*(s.total/t.total)<<"%] ";
if (t.total != 0)
os<<" ["<<100.0*(s.total/t.total)<<"%] ";
else
os<<" ["<<0<<"%] ";
os<<" ["<<s.min<<", "<<s.max<<", "<<s.mean<<" : "<<s.sigma<<"]";
if (s.mean != 0.0)
......@@ -1035,7 +1049,10 @@ avtICAlgorithm::PrintTiming(ostream &os,
v = 0.0;
os<<v;
os<<" ["<<100.0*(v/t.value)<<"%] ";
if (t.value != 0)
os<<" ["<<100.0*(v/t.value)<<"%] ";
else
os<<" ["<<0<<"%] ";
os<<endl;
}
}
......
......@@ -55,8 +55,10 @@
//
// ****************************************************************************
avtIVPSolver::avtIVPSolver() : convertToCartesian(0), convertToCylindrical(0),
order(1), h(1e-5), tol(1e-8), t(0.0), period(0)
avtIVPSolver::avtIVPSolver() : convertToCartesian(false), convertToCylindrical(false),
order(1), yCur(avtVector()), h(1e-5), h_max(1e-5),
tol(1e-8), t(0.0), period(0), baseTime(0), maxTime(1),
direction(DIRECTION_BACKWARD)
{
}
......
......@@ -163,6 +163,7 @@ avtPICSFilter::avtPICSFilter()
pathlineOverrideTime = false;
seedTimeStep0 = 0;
seedTime0 = 0.0;
baseTime = 0.0;
period = 0;
rollover = false;
......
......@@ -1718,6 +1718,8 @@ GetDataAllComponentsRange(vtkDataSet *ds, double *exts, const char *vname,
// Added 'GetNodalMagnitudeRangeViaCells' to be used when data is nodal
// and 'onlyConnectedNodes' is requested.
//
// Mark C. Miller, Wed Oct 1 19:41:32 PDT 2014
// Add some casts to double precision to avoid FPE issues.
// ****************************************************************************
template <class T> static void
......@@ -1732,7 +1734,8 @@ GetMagnitudeRange(T *buf, int n, int ncomps, double *exts,
double mag = 0.0;
for (int j = 0; j < ncomps; j++, buf++)
mag += *buf * *buf;
mag += (double) *buf * (double) *buf;
if (checkFinite)
if (! visitIsFinite(mag))
......@@ -1756,8 +1759,8 @@ GetMagnitudeRange(T *buf, int n, int ncomps, double *exts,
return GetMagnitudeRange(buf_orig, n, ncomps, exts, ghosts, true);
}
exts[0] = sqrt(exts[0]);
exts[1] = sqrt(exts[1]);
exts[0] = exts[0]>0?sqrt(exts[0]):0;
exts[1] = exts[1]>0?sqrt(exts[1]):0;
}
template <class T> static void
......
......@@ -143,6 +143,9 @@ Digits(double min, double max)
// Programmer: Eric Brugger
// Creation: December 9, 2008
//
// Modifications
// Mark C. Miller, Thu Oct 2 09:07:24 PDT 2014
// Defend against negative logs.
// ****************************************************************************
int
......@@ -153,6 +156,9 @@ LabelExponent(double min, double max)
//
double range = (fabs(min) > fabs(max) ? fabs(min) : fabs(max));
if (range <= 0)
return 0;
double pow10 = log10(range);
//
......
......@@ -229,51 +229,6 @@ InlineExtract(char * const &target, const char *&src, const int &amount)
src += amount;
}
// ****************************************************************************
// Function: VisItInitExtentsToLimits
//
// Purpose: Initialize an extents array to numerical limits of platform
//
// Programmer: Mark C. Miller
// Creation: November 15, 2007
//
// ****************************************************************************
inline bool
VisItInitExtentsToLimits(double *exts, int n)
{
if (!(n == 2 || n == 4 || n == 6))
return false;
for (int i = 0; i < n; i++)
exts[i] = n % 2 ? -DBL_MAX : DBL_MAX;
return true;
}
inline bool
VisItInitExtentsToLimits(float *exts, int n)
{
if (!(n == 2 || n == 4 || n == 6))
return false;
for (int i = 0; i < n; i++)
exts[i] = n % 2 ? -FLT_MAX : FLT_MAX;
return true;
}
inline bool
VisItInitExtentsToLimits(int *exts, int n)
{
if (!(n == 2 || n == 4 || n == 6))
return false;
for (int i = 0; i < n; i++)
exts[i] = n % 2 ? -INT_MAX : INT_MAX;
return true;
}
#endif
......@@ -76,10 +76,15 @@ using std::string;
//
// Mark C. Miller, Tue Jan 4 10:23:19 PST 2005
// Added window id
//
// Mark C. Miller, Wed Oct 1 19:55:02 PDT 2014
// Ensure extents are initialized.
// ****************************************************************************
SetWinAnnotAttsRPC::SetWinAnnotAttsRPC() : BlockingRPC("aaasaIDsi")
{
const double init_extents[6] = {0,0,0,0,0,0};
SetViewExtents(init_extents);
}
// ****************************************************************************
......
......@@ -6395,14 +6395,8 @@ ViewerWindow::SendWindowEnvironmentToEngine(const EngineKey &ek)
int fns[7];
visWindow->GetFrameAndState(fns[0], fns[1], fns[2], fns[3],
fns[4], fns[5], fns[6]);
double vexts[6];
VisItInitExtentsToLimits(vexts, 6);
double vexts[6]={0,0,0,0,0,0};
GetExtents(visWindow->GetWindowMode()==WINMODE_3D?3:2, vexts);
if (visWindow->GetWindowMode()!=WINMODE_3D)
{
vexts[4] = 0.0;
vexts[5] = 0.0;
}
return GetViewerEngineManager()->SetWinAnnotAtts(ek,
&winAtts, &annotAtts, &annotObjs, extStr, &visCues,
......
This diff is collapsed.
......@@ -117,6 +117,9 @@ vtkVectorReduceFilter::SetLimitToOriginal(bool lto)
// Eric Brugger, Thu Jan 10 10:13:10 PST 2013
// Modified to inherit from vtkPolyDataAlgorithm.
//
// Mark C. Miller, Wed Oct 1 19:42:14 PDT 2014
// Add logic to selectively InsertNextValue(double|float[3]) as well as
// safely down-cast from double to float to avoid FPE issues.
// ****************************************************************************
int
......@@ -244,7 +247,18 @@ vtkVectorReduceFilter::RequestData(
input->GetPoint(i, pt);
outpts->InsertNextPoint(pt);
outVecs->InsertNextTuple(v);
if (inctype == VTK_DOUBLE || inptype == VTK_DOUBLE)
{
outVecs->InsertNextTuple(v);
}
else
{
float fv[3];
fv[0] = vtkVisItUtility::SafeDoubleToFloat(v[0]);
fv[1] = vtkVisItUtility::SafeDoubleToFloat(v[1]);
fv[2] = vtkVisItUtility::SafeDoubleToFloat(v[2]);
outVecs->InsertNextTuple(fv);
}
outPd->CopyData(inPd, i, count++);
}
}
......
......@@ -41,7 +41,9 @@
// ************************************************************************* //
#include <algorithm>
#include <cmath>
#include <float.h>
#include <limits>
#include <list>
#include <vtkVisItUtility.h>
......@@ -1314,3 +1316,23 @@ void vtkVisItUtility::CleanupStaticVTKObjects()
std::for_each(vtkobjects.begin(), vtkobjects.end(), DeleteVTK<vtkObject>);
vtkobjects.clear();
}
// ****************************************************************************
// Function: SafeDoubleToFloat
//
// Purpose: Safely convert a double precision value to float precision
// without triggering an FPE issue.
//
// Programmer: Mark C. Miller
// Creation: October 1, 2014
// ****************************************************************************
float vtkVisItUtility::SafeDoubleToFloat(double v)
{
if (std::abs(v) < (double) std::numeric_limits<float>::min()) return 0;
if (std::abs(v) > (double) std::numeric_limits<float>::max())
{
if (v < 0) return -std::numeric_limits<float>::max();
if (v > 0) return std::numeric_limits<float>::min();
}
return (float) v;
}
......@@ -74,6 +74,8 @@ class vtkRectilinearGrid;
// Added option to allow negaitve inddices (false by default) to
// GetLogicalIndices.
//
// Mark C. Miller, Wed Oct 1 19:44:34 PDT 2014
// Add SafeDoubleToFloat
// ****************************************************************************
namespace vtkVisItUtility
......@@ -121,5 +123,6 @@ namespace vtkVisItUtility
const double *_eps = 0);
VISIT_VTK_LIGHT_API void RegisterStaticVTKObject(vtkObject*);
VISIT_VTK_LIGHT_API void CleanupStaticVTKObjects();
VISIT_VTK_LIGHT_API float SafeDoubleToFloat(double);
}
#endif
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment