Commit 9f9e134f authored by Pavel Burovskiy's avatar Pavel Burovskiy
Browse files

Changing double to NekDouble in LibUtilitites (except LinearAlgebra)

parent 601d404d
......@@ -37,8 +37,8 @@
#define NEKTAR_LIBUTILITIES_EQUATION_HPP
#include <string>
#include <map> // for map
#include <stdexcept> // for runtime_error
#include <map>
#include <stdexcept>
#include <LibUtilities/Interpreter/AnalyticExpressionEvaluator.hpp>
#include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
......@@ -191,7 +191,7 @@ namespace Nektar
}
}
LIB_UTILITIES_EXPORT void SetParameter(const std::string& name, double value)
LIB_UTILITIES_EXPORT void SetParameter(const std::string& name, NekDouble value)
{
m_evaluator.SetParameter(name, value);
}
......@@ -208,7 +208,7 @@ namespace Nektar
/// Returns time spend on expression evaluation at
/// points (it does not include parse/pre-processing time).
LIB_UTILITIES_EXPORT double GetTime() const
LIB_UTILITIES_EXPORT NekDouble GetTime() const
{
return m_evaluator.GetTime();
}
......
......@@ -71,9 +71,9 @@ namespace Nektar
struct MeshVertex
{
int id;
double x;
double y;
double z;
NekDouble x;
NekDouble y;
NekDouble z;
};
struct MeshEdge
......
///////////////////////////////////////////////////////////////////////////////
//
// File: Timer.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Scientific Computing and Imaging Institute,
// University of Utah (USA) and Department of Aeronautics, Imperial
// College London (UK).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Time getting class
//
///////////////////////////////////////////////////////////////////////////////
#include <LibUtilities/BasicUtils/Timer.h>
......@@ -71,21 +105,21 @@ namespace Nektar
#endif
}
double Timer::TimePerTest(unsigned int n)
NekDouble Timer::TimePerTest(unsigned int n)
{
#ifdef _WIN32
CounterType frequency;
QueryPerformanceFrequency(&frequency);
return Elapsed().QuadPart/static_cast<double>(n) * 1.0/frequency.QuadPart;
return Elapsed().QuadPart/static_cast<NekDouble>(n) * 1.0/frequency.QuadPart;
#elif defined(__APPLE__)
CounterType elapsed = Elapsed();
double result = elapsed.tv_sec/static_cast<double>(n) +
( elapsed.tv_usec/static_cast<double>(n) * 1.0e-6);
NekDouble result = elapsed.tv_sec/static_cast<NekDouble>(n) +
( elapsed.tv_usec/static_cast<NekDouble>(n) * 1.0e-6);
return result;
#else
CounterType elapsed = Elapsed();
double result = elapsed.tv_sec/static_cast<double>(n) +
( elapsed.tv_nsec/static_cast<double>(n) * 1.0e-9);
NekDouble result = elapsed.tv_sec/static_cast<NekDouble>(n) +
( elapsed.tv_nsec/static_cast<NekDouble>(n) * 1.0e-9);
return result;
#endif
}
......
///////////////////////////////////////////////////////////////////////////////
//
// File: Timer.h
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
......@@ -27,7 +29,7 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description:
// Description: Time getting class
//
///////////////////////////////////////////////////////////////////////////////
......@@ -43,6 +45,7 @@
#endif
#include <LibUtilities/LibUtilitiesDeclspec.h>
#include <LibUtilities/BasicConst/NektarUnivConsts.hpp>
namespace Nektar
{
......@@ -64,7 +67,10 @@ namespace Nektar
LIB_UTILITIES_EXPORT void Start();
LIB_UTILITIES_EXPORT void Stop();
LIB_UTILITIES_EXPORT CounterType Elapsed();
LIB_UTILITIES_EXPORT double TimePerTest(unsigned int n);
/// \brief Returns amount of seconds per iteration in
/// a test with n iterations.
LIB_UTILITIES_EXPORT NekDouble TimePerTest(unsigned int n);
private:
Timer(const Timer& rhs);
......
......@@ -945,7 +945,7 @@ namespace Vmath
const Nektar::NekDouble *x, const int incx,
const int *y, const int incy);
/*
// \brief copy one int vector to another
void Vcopy(int n, const int *x, const int incx, int *y,
const int incy)
......@@ -1002,6 +1002,32 @@ namespace Vmath
}
}
}
*/
// \brief copy one vector to another
template<typename T>
void Vcopy(int n, const T *x, const int incx,
T *y, const int incy)
{
if( incx ==1 && incy == 1)
{
memcpy(y,x,n*sizeof(T));
}
else
{
while( n-- )
{
*y = *x;
x += incx;
y += incy;
}
}
}
template LIB_UTILITIES_EXPORT void Vcopy( int n, const int *x, const int incx, int *y, const int incy);
template LIB_UTILITIES_EXPORT void Vcopy( int n, const unsigned int *x, const int incx, unsigned int *y, const int incy);
template LIB_UTILITIES_EXPORT void Vcopy( int n, const Nektar::NekDouble *x, const int incx, Nektar::NekDouble *y, const int incy);
// \brief reverse the ordering of vector to another
template<class T> void Reverse( int n, const T *x, const int incx, T *y, const int incy)
......
......@@ -273,31 +273,14 @@ namespace Vmath
/********** Memory routines ***********************/
#if 0
// \brief copy one double vector to another - This is just a wrapper
// around Blas
static inline void Vcopy(int n, const double *x, int incx, double *y,
const int incy)
{
Blas::Dcopy(n,x,incx,y,incy);
}
#else
// \brief copy one int vector to another
LIB_UTILITIES_EXPORT void Vcopy(int n, const int *x, const int incx, int *y,
const int incy);
// \brief copy one int vector to another
LIB_UTILITIES_EXPORT void Vcopy(int n, const unsigned int *x, const int incx, unsigned int *y,
const int incy);
// \brief copy one double vector to another
LIB_UTILITIES_EXPORT void Vcopy(int n, const double *x, const int incx, double *y,
const int incy);
// \brief copy one vector to another
template<class T>
LIB_UTILITIES_EXPORT void Vcopy(int n, const T *x, const int incx,
T *y, const int incy);
// \brief reverse the ordering of vector to another
template<class T> LIB_UTILITIES_EXPORT void Reverse( int n, const T *x, const int incx, T *y, const int incy);
#endif
}
#endif //VECTORMATH_HPP
......
......@@ -56,7 +56,7 @@ namespace Nektar
m_FFTW_w = Array<OneD,NekDouble>(m_N);
m_FFTW_w_inv = Array<OneD,NekDouble>(m_N);
m_FFTW_w[0] = 1.0/(double)m_N;
m_FFTW_w[0] = 1.0/(NekDouble)m_N;
m_FFTW_w[1] = m_FFTW_w[0];
m_FFTW_w_inv[0] = m_N;
......
......@@ -45,7 +45,7 @@ namespace Nektar
namespace LibUtilities
{
// Default value.
double BLPoints::delta_star = 1.0;
NekDouble BLPoints::delta_star = 1.0;
void BLPoints::CalculatePoints()
{
......@@ -54,13 +54,13 @@ namespace Nektar
unsigned int npts = m_pointsKey.GetNumPoints();
// Derived power coefficient.
double rn = powf(2.0/BLPoints::delta_star,1.0/(npts-2.0));
NekDouble rn = powf(2.0/BLPoints::delta_star,1.0/(npts-2.0));
m_points[0][0] = -1.0;
for (unsigned int i = 1; i < npts; ++i)
{
m_points[0][i] = -1.0 + delta_star*pow(rn,(double)i-1.0);
m_points[0][i] = -1.0 + delta_star*pow(rn,(NekDouble)i-1.0);
}
}
......
......@@ -113,7 +113,7 @@ namespace Nektar
boost::bind(&BLPoints::CreateMatrix, this, _1));
}
LIB_UTILITIES_EXPORT static double delta_star;
LIB_UTILITIES_EXPORT static NekDouble delta_star;
private:
/// Default constructor should not be called except by Create
......
......@@ -60,7 +60,7 @@ namespace Nektar
}
else
{
double dx = 2.0/(double)(npts);
NekDouble dx = 2.0/(NekDouble)(npts);
for(unsigned int i=0;i<npts;++i)
{
m_points[0][i] = -1.0 + i*dx;
......@@ -82,7 +82,7 @@ namespace Nektar
{
for(unsigned int i=0; i<npts; ++i)
{
m_weights[i] = 1.0/(double)npts;
m_weights[i] = 1.0/(NekDouble)npts;
}
}
}
......@@ -101,7 +101,7 @@ namespace Nektar
for(unsigned int i=1;i<npts;++i)
{
m_derivmatrix[0]->SetValue(0,i, -0.5*M_PI*pow(-1.0,double(i))*cos(M_PI*i/npts)/sin(M_PI*i/npts));
m_derivmatrix[0]->SetValue(0,i, -0.5*M_PI*pow(-1.0,NekDouble(i))*cos(M_PI*i/npts)/sin(M_PI*i/npts));
}
for(unsigned int i=1;i<npts;++i)
......@@ -174,13 +174,13 @@ namespace Nektar
}
}
double FourierPoints::PeriodicSincFunction(const NekDouble x, const NekDouble h)
NekDouble FourierPoints::PeriodicSincFunction(const NekDouble x, const NekDouble h)
{
// This formula is based upon a mapped version of
// the periodic sinc presented in Trefethen's "Spectral Methods
// in Matlab"
double y = 1.0;
NekDouble y = 1.0;
if(fabs(x)>1.0e-12)
{
......
......@@ -107,7 +107,7 @@ namespace Nektar
void CalculateInterpMatrix(unsigned int npts, const Array<OneD, const NekDouble>& xpoints,
Array<OneD, NekDouble>& interp);
double PeriodicSincFunction(const NekDouble x, const NekDouble h);
NekDouble PeriodicSincFunction(const NekDouble x, const NekDouble h);
}; // class FourierPoints
} // end of namespace
} // end of namespace
......
......@@ -110,7 +110,7 @@ namespace Nektar
void CalculateInterpMatrix(unsigned int npts, const Array<OneD, const NekDouble>& xpoints,
Array<OneD, NekDouble>& interp);
double PeriodicSincFunction(const NekDouble x, const NekDouble h);
NekDouble PeriodicSincFunction(const NekDouble x, const NekDouble h);
}; // class FourierPoints
} // end of namespace
} // end of namespace
......
......@@ -51,7 +51,7 @@ namespace Nektar
void NodalTetElec::CalculatePoints()
{
// Allocate the storage for points
Points<double>::CalculatePoints();
Points<NekDouble>::CalculatePoints();
int index=0,isum=0;
const int offset = 5; //offset to match Datafile
......
......@@ -22,7 +22,7 @@ namespace Nektar
const unsigned int NodalTetElecAvailable = 10;
static const unsigned int NodalTetElecNPTS[NodalTetElecAvailable] = {1,2,3,5,6,9,11,15,18,23};
static const double NodalTetElecData[][9] = {
static const NekDouble NodalTetElecData[][9] = {
// %%% n_1 n_4 n_6 n_12 n_24 l_1 l_2 l_3 l_4
// 1 1 %%% Order / Number of Points
{0, 1, 0, 0, 0, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000},
......
......@@ -9,7 +9,7 @@ namespace Nektar
{2,1,0},{0,2,1},{1,2,0}};
const unsigned int NodalTriElecAvailable = 16;
static const unsigned int NodalTriElecNPTS[NodalTriElecAvailable] = {1,2,3,4,5,7,8,10,12,14,16,19,21,24,27,30};
static const double NodalTriElecData[][6] = {
static const NekDouble NodalTriElecData[][6] = {
// %%% n_1 n_3 n_6 l_1 l_2 l_3
// 1 1 %%% Order / Number of Points
......
......@@ -9,7 +9,7 @@ namespace Nektar
{2,1,0},{0,2,1},{1,2,0}}; // Works for abc
const unsigned int NodalTriFeketeAvailable = 16;
static const unsigned int NodalTriFeketeNPTS[NodalTriFeketeAvailable] = {1,2,3,4,5,7,8,10,12,14,16,19,21,24,27,30};
static const double NodalTriFeketeData[][6] = {
static const NekDouble NodalTriFeketeData[][6] = {
// %%% n_1 n_3 n_6 l_1 l_2 l_3
// 1 1 %%% Order / Number of Points
{0, 2, 0, 1.0000000000, 0.0000000000, 0.0000000000},
......
......@@ -154,7 +154,7 @@ namespace Nektar
}
// Formatted version of matrix ostream output
std::string MatrixToString( const NekMatrix<NekDouble> & A, int precision, double expSigFigs )
std::string MatrixToString( const NekMatrix<NekDouble> & A, int precision, NekDouble expSigFigs )
{
stringstream s;
s << setprecision(precision);
......@@ -164,7 +164,7 @@ namespace Nektar
{
for(int j=0; j<N; ++j)
{
double a = MakeRound(expSigFigs * A(i, j)) / expSigFigs;
NekDouble a = MakeRound(expSigFigs * A(i, j)) / expSigFigs;
s << setw(7) << right << a;
if( j < N-1 )
{
......@@ -180,14 +180,14 @@ namespace Nektar
}
// Formatted version of vector ostream output
std::string VectorToString( const NekVector<NekDouble> & v, int precision, double expSigFigs )
std::string VectorToString( const NekVector<NekDouble> & v, int precision, NekDouble expSigFigs )
{
stringstream s;
s << setprecision(precision) << "[ ";
int N = int(v.GetRows());
for(int j=0; j<N; ++j )
{
double x = MakeRound(expSigFigs * v(j)) / expSigFigs;
NekDouble x = MakeRound(expSigFigs * v(j)) / expSigFigs;
s << setw(7) << right << x;
if( j < N-1 )
{
......@@ -221,7 +221,7 @@ namespace Nektar
// Get Tetrahedral number, where Tn = (d+1)(d+2)(d+3)/6
int GetTetDegree(int nBasisFunc)
{
double eq = pow( 81.0 * nBasisFunc + 3.0 * sqrt(-3.0 + 729.0 * nBasisFunc * nBasisFunc), 1.0/3.0);
NekDouble eq = pow( 81.0 * nBasisFunc + 3.0 * sqrt(-3.0 + 729.0 * nBasisFunc * nBasisFunc), 1.0/3.0);
int degree = int(MakeRound(eq/3.0 + 1.0/eq - 1.0)) - 1;
ASSERTL1( GetTetNumPoints(degree) == nBasisFunc, "The number of points defines an expansion of fractional degree, which is not supported." );
......@@ -422,7 +422,7 @@ namespace Nektar
// Initialize the horizontal coordinate of the Tetrahedral
for(int el=0; el<size; ++el)
{
if( z[el] < -y[el] - numeric_limits<double>::epsilon() )
if( z[el] < -y[el] - numeric_limits<NekDouble>::epsilon() )
{
eta_1[el] = 2.0*(1.0 + x[el])/(-y[el]-z[el]) - 1.0;
}
......@@ -435,7 +435,7 @@ namespace Nektar
// Initialize the coordinate of the Tetrahedral
for(int el=0; el<size; ++el)
{
if( z[el] < 1.0 - numeric_limits<double>::epsilon())
if( z[el] < 1.0 - numeric_limits<NekDouble>::epsilon())
{
eta_2[el] = 2.0*(1.0 + y[el]) / (1.0 - z[el]) - 1.0;
}
......@@ -719,7 +719,7 @@ namespace Nektar
// Initialize the horizontal coordinate of the Tetrahedral (beta in Barycentric coordinate)
for(int el=0; el<size; ++el)
{
if( y[el] < -z[el] - numeric_limits<double>::epsilon())
if( y[el] < -z[el] - numeric_limits<NekDouble>::epsilon())
{
eta_1[el] = 2.0*(1.0 + x[el])/(-y[el]-z[el]) - 1.0;
} else
......@@ -731,7 +731,7 @@ namespace Nektar
// Initialize the coordinate of the Tetrahedral (gamma in Barycentric coordinate)
for(int el=0; el<size; ++el)
{
if( z[el] < 1.0 - numeric_limits<double>::epsilon())
if( z[el] < 1.0 - numeric_limits<NekDouble>::epsilon())
{
eta_2[el] = 2.0*(1.0 + y[el]) / (1.0 - z[el]) - 1.0;
}
......@@ -894,7 +894,7 @@ namespace Nektar
// Initialize the collapsed horizontal coordinate of the Tetrahedral
for(int el=0; el<size; ++el)
{
if( y[el] < -z[el] - numeric_limits<double>::epsilon())
if( y[el] < -z[el] - numeric_limits<NekDouble>::epsilon())
{
eta_1[el] = 2.0*(1.0 + x[el])/(-y[el]-z[el]) - 1.0;
eta_1_dy[el] = 2.0*(1.0 + x[el])/((y[el]+z[el])*(y[el]+z[el]));
......@@ -909,7 +909,7 @@ namespace Nektar
// Initialize the collapsed depth coordinate of the Tetrahedral
for(int el=0; el<size; ++el)
{
if( z[el] < 1.0 - numeric_limits<double>::epsilon())
if( z[el] < 1.0 - numeric_limits<NekDouble>::epsilon())
{
eta_2[el] = 2.0*(1.0 + y[el]) / (1.0 - z[el]) - 1.0;
eta_2_dy[el] = 2.0/(1.0 - z[el]);
......@@ -974,7 +974,7 @@ namespace Nektar
// Fix singularity at z=1
for(int k=0; k<size; ++k)
{
if( z[k] >= 1.0 - numeric_limits<double>::epsilon() )
if( z[k] >= 1.0 - numeric_limits<NekDouble>::epsilon() )
{
if( p + q > 0 )
{
......@@ -1071,7 +1071,7 @@ namespace Nektar
// Initialize the collapsed horizontal coordinate of the Tetrahedral
for(int el=0; el<size; ++el)
{
if( y[el] < -z[el] - numeric_limits<double>::epsilon())
if( y[el] < -z[el] - numeric_limits<NekDouble>::epsilon())
{
eta_1[el] = 2.0*(1.0 + x[el]) / (-y[el]-z[el]) - 1.0;
eta_1_dz[el] = 2.0*(1.0 + x[el]) / ((y[el]+z[el])*(y[el]+z[el]));
......@@ -1086,7 +1086,7 @@ namespace Nektar
// Initialize the collapsed depth coordinate of the Tetrahedral
for(int el=0; el<size; ++el)
{
if( z[el] < 1.0 - numeric_limits<double>::epsilon())
if( z[el] < 1.0 - numeric_limits<NekDouble>::epsilon())
{
eta_2[el] = 2.0*(1.0 + y[el]) / (1.0 - z[el]) - 1.0;
eta_2_dz[el] = 2.0*(1.0 + y[el]) / ((1.0 - z[el])*(1.0 - z[el]));
......
......@@ -68,8 +68,8 @@ namespace Nektar
LIB_UTILITIES_EXPORT Array<OneD, NekDouble> ToArray( const NekVector<NekDouble> & x );
LIB_UTILITIES_EXPORT NekVector<NekDouble> Hadamard( const NekVector<NekDouble> & x, const NekVector<NekDouble> & y );
LIB_UTILITIES_EXPORT NekVector<NekDouble> VectorPower( const NekVector<NekDouble> & x, NekDouble p );
LIB_UTILITIES_EXPORT std::string MatrixToString( const NekMatrix<NekDouble> & A, int precision = 2, double threshold = 1e12 );
LIB_UTILITIES_EXPORT std::string VectorToString( const NekVector<NekDouble> & v, int precision = 2, double threshold = 1e12 );
LIB_UTILITIES_EXPORT std::string MatrixToString( const NekMatrix<NekDouble> & A, int precision = 2, NekDouble threshold = 1e12 );
LIB_UTILITIES_EXPORT std::string VectorToString( const NekVector<NekDouble> & v, int precision = 2, NekDouble threshold = 1e12 );
// ////////////////////////////////////////////////////////////////
// Polynomials(Jacobi, Legendre, and Dubiner) and its derivations
......
......@@ -71,11 +71,11 @@ namespace Nektar
// AnalyticExpression definitions for Spirit Parser
// =========================================================================
typedef double (*PFD)();
typedef double (*PFD1)(double);
typedef double (*PFD2)(double, double);
typedef double (*PFD3)(double, double, double);
typedef double (*PFD4)(double, double, double, double);
typedef NekDouble (*PFD)();
typedef NekDouble (*PFD1)(NekDouble);
typedef NekDouble (*PFD2)(NekDouble, NekDouble);
typedef NekDouble (*PFD3)(NekDouble, NekDouble, NekDouble);
typedef NekDouble (*PFD4)(NekDouble, NekDouble, NekDouble, NekDouble);
struct func
{
func(PFD1 p) : func1(p), size(1) {};
......@@ -94,7 +94,7 @@ namespace Nektar
};
// signum function
double sign(double arg)
NekDouble sign(NekDouble arg)
{
return (arg > 0.0) - (arg < 0.0);
}
......@@ -103,7 +103,7 @@ namespace Nektar
// Arg: sigma of the zero-mean gaussian distribution
// Attention: this function is not actually used for
// evaluation purposes.
double awgn(double sigma)
NekDouble awgn(NekDouble sigma)
{
AnalyticExpressionEvaluator::RandomGeneratorType rng;
boost::variate_generator<
......@@ -116,7 +116,7 @@ namespace Nektar
/** This struct creates a parser that matches the function
definitions from math.h. All of the functions accept one
of more doubles as arguments and returns a double. **/
of more NekDoubles as arguments and returns a NekDouble. **/
static struct functions : symbols<func>
{
functions()
......@@ -140,7 +140,7 @@ namespace Nektar
("sqrt", sqrt)
("tan", tan)
("tanh", tanh)
// and one more
// and few more custom functions
("sign", sign)
("awgn", awgn)
;
......@@ -302,15 +302,15 @@ namespace Nektar
}
void AnalyticExpressionEvaluator::AddConstants(std::map<std::string, double> const& constants)
void AnalyticExpressionEvaluator::AddConstants(std::map<std::string, NekDouble> const& constants)
{
for (std::map<std::string, double>::const_iterator it = constants.begin(); it != constants.end(); ++it)
for (std::map<std::string, NekDouble>::const_iterator it = constants.begin(); it != constants.end(); ++it)
{
AddConstant(it->first, it->second);
}
}
int AnalyticExpressionEvaluator::AddConstant(std::string const& name, double value)
int AnalyticExpressionEvaluator::AddConstant(std::string const& name, NekDouble value)
{
ConstantMap::const_iterator it = m_constantMapNameToId.find(name);
if (it == m_constantMapNameToId.end())
......@@ -335,24 +335,24 @@ namespace Nektar
return it->second;
}
double AnalyticExpressionEvaluator::GetConstant(std::string const& name)
NekDouble AnalyticExpressionEvaluator::GetConstant(std::string const& name)
{
double* value = find(m_constantsParser, name.c_str());
NekDouble* value = find(m_constantsParser, name.c_str());
ASSERTL1(value != NULL, "Constant variable not found: " + name);