Skip to content

GitLab

Attention une mise à jour du serveur va être effectuée le vendredi 16 avril entre 12h et 12h30. Cette mise à jour va générer une interruption du service de quelques minutes.

Commit 1032efb7 authored by berenger-bramas's avatar berenger-bramas
Browse files
Options

Add the compilation directive : -Wconversion 

Now there is possible "auto" cast.

git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/scalfmm/scalfmm/trunk@189 2616d619-271b-44dc-8df4-d4a8f33a7222
parent 524af421
Hide whitespace changes
Inline Side-by-side
Showing with 224 additions and 186 deletions
CMakeLists.txt
View file @ 1032efb7
......@@ -39,7 +39,7 @@ if( SCALFMM_USE_CBLAS )
endif()
# Compile option
ADD_DEFINITIONS(-O2 -Wall -Wshadow -Wpointer-arith -Wcast-qual)
ADD_DEFINITIONS(-O2 -Wall -Wshadow -Wpointer-arith -Wcast-qual -Wconversion)
CONFIGURE_FILE( ${CMAKE_SOURCE_DIR}/Src/ScalFmmConfig.h.cmake
${CMAKE_BINARY_DIR}/Src/ScalFmmConfig.h )
......
Src/Containers/FOctree.hpp
View file @ 1032efb7
......@@ -76,8 +76,8 @@ class FOctree {
*/
long getTreeCoordinate(const FReal inRelativePosition) const {
const FReal indexFReal = inRelativePosition / this->boxWidthAtLevel[this->leafIndex];
const long index = FMath::dfloor(indexFReal);
if( index && FMath::LookEqual(inRelativePosition, this->boxWidthAtLevel[this->leafIndex] * index ) ){
const long index = long(FMath::dfloor(indexFReal));
if( index && FMath::LookEqual(inRelativePosition, this->boxWidthAtLevel[this->leafIndex] * FReal(index) ) ){
return index - 1;
}
return index;
......@@ -99,7 +99,7 @@ public:
// pre compute box width for each level
for(int indexLevel = 0; indexLevel < this->height; ++indexLevel ){
this->boxWidthAtLevel[indexLevel] = tempWidth;
tempWidth /= 2.0;
tempWidth /= FReal(2.0);
}
}
......@@ -674,7 +674,7 @@ public:
// compute the leaf index
const MortonIndex fullIndex = inIndex >> 3 * (inLevel + 1 - (workingTree.tree->getSubOctreeHeight() + workingTree.tree->getSubOctreePosition()) );
// point to next suboctree
workingTree.tree = workingTree.middleTree->leafs(treeMiddleMask & fullIndex);
workingTree.tree = workingTree.middleTree->leafs(int(treeMiddleMask & fullIndex));
if(!workingTree.tree) return 0;
}
......@@ -774,13 +774,13 @@ public:
// compute the leaf index
const MortonIndex fullIndex = inIndex >> (3 * (leafIndex + 1 - (workingTree.tree->getSubOctreeHeight() + workingTree.tree->getSubOctreePosition()) ) );
// point to next suboctree
workingTree.tree = workingTree.middleTree->leafs(treeSubLeafMask & fullIndex);
workingTree.tree = workingTree.middleTree->leafs(int(treeSubLeafMask & fullIndex));
if(!workingTree.tree) return 0;
}
// compute correct index in the array
const MortonIndex treeLeafMask = ~(~0x00LL << (3 * (leafIndex + 1 - workingTree.tree->getSubOctreePosition()) ));
return workingTree.leafTree->getLeafSrc(treeLeafMask & inIndex);
return workingTree.leafTree->getLeafSrc(int(treeLeafMask & inIndex));
}
/** This function fill an array with the neighbors of a cell
......
Src/Containers/FSubOctree.hpp
View file @ 1032efb7
......@@ -88,9 +88,9 @@ protected:
const int realLevel = indexLevel + this->getSubOctreePosition();
const FReal widthAtLevel = inBoxWidthAtLevel[realLevel];
newNode->setPosition(F3DPosition(
treePosition.getX() * widthAtLevel + widthAtLevel/2.0,
treePosition.getY() * widthAtLevel + widthAtLevel/2.0,
treePosition.getZ() * widthAtLevel + widthAtLevel/2.0));
FReal(treePosition.getX()) * widthAtLevel + widthAtLevel/ FReal(2.0),
FReal(treePosition.getY()) * widthAtLevel + widthAtLevel/ FReal(2.0),
FReal(treePosition.getZ()) * widthAtLevel + widthAtLevel/ FReal(2.0) ));
this->cells[indexLevel][arrayIndex] = newNode;
......@@ -324,7 +324,7 @@ public:
if( !this->leafs[arrayIndex] ){
this->leafs[arrayIndex] = new LeafClass();
FAbstractSubOctree<ParticleClass,CellClass,ContainerClass,LeafClass>::newLeafInserted( arrayIndex , index, host, inBoxWidthAtLevel);
FAbstractSubOctree<ParticleClass,CellClass,ContainerClass,LeafClass>::newLeafInserted( long(arrayIndex) , index, host, inBoxWidthAtLevel);
}
// add particle to leaf list
this->leafs[arrayIndex]->push(inParticle);
......@@ -439,11 +439,11 @@ public:
// Next suboctree is a middle suboctree
if(inTreeHeight > nextSubOctreeHeight + nextSubOctreePosition){
this->subleafs[arrayIndex] = new FSubOctree(this,arrayIndex,nextSubOctreeHeight,nextSubOctreePosition);
this->subleafs[arrayIndex] = new FSubOctree(this,long(arrayIndex),nextSubOctreeHeight,nextSubOctreePosition);
}
// Or next suboctree contains the reail leaf!
else{
this->subleafs[arrayIndex] = new FSubOctreeWithLeafs<ParticleClass,CellClass,ContainerClass,LeafClass>(this,arrayIndex,nextSubOctreeHeight,nextSubOctreePosition);
this->subleafs[arrayIndex] = new FSubOctreeWithLeafs<ParticleClass,CellClass,ContainerClass,LeafClass>(this,long(arrayIndex),nextSubOctreeHeight,nextSubOctreePosition);
}
const FTreeCoordinate hostAtLevel(
......@@ -452,7 +452,7 @@ public:
host.getZ() >> (inTreeHeight - nextSubOctreePosition ));
// We need to inform parent class
FAbstractSubOctree<ParticleClass,CellClass,ContainerClass,LeafClass>::newLeafInserted( arrayIndex, index >> (3 * (inTreeHeight-nextSubOctreePosition) ), hostAtLevel, inBoxWidthAtLevel);
FAbstractSubOctree<ParticleClass,CellClass,ContainerClass,LeafClass>::newLeafInserted( long(arrayIndex), index >> (3 * (inTreeHeight-nextSubOctreePosition) ), hostAtLevel, inBoxWidthAtLevel);
}
// Ask next suboctree to insert the particle
this->subleafs[arrayIndex]->insert( index, host, inParticle, inTreeHeight, inBoxWidthAtLevel);
......
Src/Containers/FTreeCoordinate.hpp
View file @ 1032efb7
......@@ -159,11 +159,11 @@ public:
this->z = 0;
for(int indexLevel = 0; indexLevel < inLevel ; ++indexLevel){
z |= (inIndex & mask);
z |= long(inIndex & mask);
inIndex >>= 1;
y |= (inIndex & mask);
y |= long(inIndex & mask);
inIndex >>= 1;
x |= (inIndex & mask);
x |= long(inIndex & mask);
mask <<= 1;
}
......
Src/Containers/FVector.hpp
View file @ 1032efb7
......@@ -67,7 +67,7 @@ public:
if(other.getSize() > this->capacity){
delete [] reinterpret_cast< char* >(this->array);
this->capacity = other.getSize() * 1.5;
this->capacity = int(other.getSize() * 1.5);
array = reinterpret_cast< T* >( new char[SizeOfT*this->capacity]);
}
this->index = other.index;
......@@ -165,7 +165,7 @@ public:
*/
void push( const T & inValue ){
if( this->index == this->capacity ){
setCapacity(this->capacity * 1.5);
setCapacity(int(this->capacity * 1.5));
}
new((void*)&this->array[this->index]) T;
this->array[this->index] = inValue;
......
Src/Core/FFmmAlgorithmThread.hpp
View file @ 1032efb7
......@@ -323,9 +323,9 @@ public:
#pragma omp parallel
{
const float step = (float(this->leafsNumber) / omp_get_num_threads());
const int start = int(FMath::Ceil(step * omp_get_thread_num()));
const int tempEnd = int(FMath::Ceil(step * (omp_get_thread_num()+1)));
const float step = float(this->leafsNumber) / float(omp_get_num_threads());
const int start = int(FMath::Ceil(step * float(omp_get_thread_num())));
const int tempEnd = int(FMath::Ceil(step * float(omp_get_thread_num()+1)));
const int end = (tempEnd > this->leafsNumber ? this->leafsNumber : tempEnd);
typename OctreeClass::Iterator octreeIterator(tree);
......
Src/Core/FFmmAlgorithmThreadProc.hpp
View file @ 1032efb7
......@@ -355,7 +355,7 @@ public:
if( child[idxChild] && workingIntervalsPerLevel[(idxLevel+1) * nbProcess + idProcess].min <= child[idxChild]->getMortonIndex() ){
child[idxChild]->serializeUp(&sendBuffer[idxBuff]);
idxBuff += CellClass::SerializedSizeUp;
state |= (0x1 << idxChild);
state |= char(0x1 << idxChild);
}
}
sendBuffer[0] = state;
......@@ -418,7 +418,7 @@ public:
// retreive data and merge my child and the child from others
for(int idxProc = firstProcThatSend ; idxProc < endProcThatSend ; ++idxProc){
char state = recvBuffer[idxProc * recvBufferOffset];
int state = int(recvBuffer[idxProc * recvBufferOffset]);
int position = 0;
int bufferIndex = 1;
......
Src/Files/FMpiTreeBuilder.hpp
View file @ 1032efb7
......@@ -17,8 +17,8 @@ class FMpiTreeBuilder{
*/
static long getTreeCoordinate(const FReal inRelativePosition, const FReal boxWidthAtLeafLevel) {
const FReal indexFReal = inRelativePosition / boxWidthAtLeafLevel;
const long index = FMath::dfloor(indexFReal);
if( index && FMath::LookEqual(inRelativePosition, boxWidthAtLeafLevel * index ) ){
const long index = long(FMath::dfloor(indexFReal));
if( index && FMath::LookEqual(inRelativePosition, boxWidthAtLeafLevel * FReal(index) ) ){
return index - 1;
}
return index;
......@@ -48,14 +48,14 @@ class FMpiTreeBuilder{
template< class T >
static T GetLeft(const T inSize) {
const float step = (float(inSize) / MpiGetNbProcs());
return T(FMath::Ceil(step * MpiGetRank()));
const float step = (float(inSize) / float(MpiGetNbProcs()));
return T(FMath::Ceil(step * float(MpiGetRank())));
}
template< class T >
static T GetRight(const T inSize) {
const float step = (float(inSize) / MpiGetNbProcs());
const T res = T(FMath::Ceil(step * (MpiGetRank()+1)));
const float step = (float(inSize) / float(MpiGetNbProcs()));
const T res = T(FMath::Ceil(step * float(MpiGetRank()+1)));
if(res > inSize) return inSize;
else return res;
}
......@@ -132,7 +132,7 @@ public:
F3DPosition boxCorner(loader.getCenterOfBox() - (loader.getBoxWidth()/2));
FTreeCoordinate host;
const FReal boxWidthAtLeafLevel = loader.getBoxWidth() / (1 << (NbLevels - 1) );
const FReal boxWidthAtLeafLevel = loader.getBoxWidth() / FReal(1 << (NbLevels - 1) );
for(long idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
loader.fillParticle(realParticlesIndexed[idxPart].particle);
host.setX( getTreeCoordinate( realParticlesIndexed[idxPart].particle.getPosition().getX() - boxCorner.getX(), boxWidthAtLeafLevel ));
......
Src/Fmb/FFmbKernels.hpp
View file @ 1032efb7
......@@ -116,7 +116,7 @@ protected:
//std::cout << "sphereHarmoOuterCoef\n";
FReal factOuter = 1.0;
// in FMB code stoped at <= FMB_Info_M2L_P but this is not sufficient
for(int idxP = 0 ; idxP <= FMB_Info_M2L_P; factOuter *= (++idxP) ){
for(int idxP = 0 ; idxP <= FMB_Info_M2L_P; factOuter *= FReal(++idxP) ){
this->sphereHarmoOuterCoef[idxP] = factOuter;
//printf("spherical_harmonic_Outer_coefficients_array %e\n",this->sphereHarmoOuterCoef[idxP]);
//printf("fact_l %e\n",factOuter);
......@@ -128,9 +128,9 @@ protected:
FReal factInner = 1.0;
FReal powN1idxP = 1.0;
//std::cout << "sphereHarmoInnerCoef\n";
for(int idxP = 0 ; idxP <= this->FMB_Info_M2L_P ; factInner *= (++idxP), powN1idxP = -powN1idxP){
for(int idxMP = 0, fact_l_m = factInner; idxMP <= idxP ; fact_l_m *= idxP+(++idxMP), ++currentInner){
*currentInner = powN1idxP / fact_l_m;
for(int idxP = 0 ; idxP <= this->FMB_Info_M2L_P ; factInner *= FReal(++idxP), powN1idxP = -powN1idxP){
for(int idxMP = 0, fact_l_m = int(factInner); idxMP <= idxP ; fact_l_m *= idxP+(++idxMP), ++currentInner){
*currentInner = powN1idxP / FReal(fact_l_m);
//std::cout << (*currentInner) << "\n";
}
}
......@@ -261,18 +261,18 @@ protected:
for(int idxl = 2; idxl <= lmax ; ++idxl ){
for( int idxm = 0; idxm <= idxl - 2 ; ++idxm , ++idxCurrent , ++idxCurrent1m , ++idxCurrent2m ){
// Compute P_l^m, l >= m+2, using (1) and store it into results_array:
outResults[idxCurrent] = (inCosTheta * ( 2 * idxl - 1 ) * outResults[idxCurrent1m] - ( idxl + idxm - 1 )
* outResults[idxCurrent2m] ) / ( idxl - idxm );
outResults[idxCurrent] = (inCosTheta * FReal( 2 * idxl - 1 ) * outResults[idxCurrent1m] - FReal( idxl + idxm - 1 )
* outResults[idxCurrent2m] ) / FReal( idxl - idxm );
}
// p_results_array_l_minus_1_m now points on P_{l-1}^{l-1}
// Compute P_l^{l-1} using (3) and store it into ptrResults:
outResults[idxCurrent++] = inCosTheta * ( 2 * idxl - 1 ) * outResults[idxCurrent1m];
outResults[idxCurrent++] = inCosTheta * FReal( 2 * idxl - 1 ) * outResults[idxCurrent1m];
// Compute P_l^l using (2 bis) and store it into results_array:
outResults[idxCurrent++] = fact * invSinTheta * outResults[idxCurrent1m];
fact += 2.0;
fact += FReal(2.0);
++idxCurrent1m;
}
}
......@@ -328,7 +328,7 @@ protected:
FComplexe cosSin[FMB_Info_M2L_P + 1];
for(int idxl = 0, idxlMod4 = 0; idxl <= FMB_Info_M2L_P ; ++idxl, ++idxlMod4){
if(idxlMod4 == 4) idxlMod4 = 0;
const FReal angle = idxl * inSphere.phi + this->PiArrayOuter[idxlMod4];
const FReal angle = FReal(idxl) * inSphere.phi + this->PiArrayOuter[idxlMod4];
cosSin[idxl].setReal( FMath::Sin(angle + FMath::FPiDiv2) );
cosSin[idxl].setImag( FMath::Sin(angle) );
......@@ -403,7 +403,7 @@ protected:
// Initialization of precomputed_cos_and_sin_array:
for(int idxm = 0 , idxmMod4 = 0; idxm <= FMB_Info_P ; ++idxm, ++idxmMod4){
if(idxmMod4 == 4) idxmMod4 = 0;
const FReal angle = idxm*inSphere.phi + PiArrayInner[idxmMod4];
const FReal angle = FReal(idxm) *inSphere.phi + PiArrayInner[idxmMod4];
cosSin[idxm].setReal(FMath::Sin(angle + FMath::FPiDiv2));
cosSin[idxm].setImag(FMath::Sin(angle));
......@@ -452,8 +452,8 @@ protected:
);*/
// Computation of {\partial Inner_l^m(r, theta, phi)}/{\partial theta}:
magnitude = (*p_spherical_harmonic_Inner_coefficients_array) * r_l * ((l*inSphere.cosTheta*(*ptr_associated_Legendre_function_Array)
- (l+m)*(*(start_ptr_associated_Legendre_function_Array + expansion_Redirection_array_for_j[l-1] + m))) / inSphere.sinTheta);
magnitude = (*p_spherical_harmonic_Inner_coefficients_array) * r_l * ((FReal(l)*inSphere.cosTheta*(*ptr_associated_Legendre_function_Array)
- FReal(l+m)*(*(start_ptr_associated_Legendre_function_Array + expansion_Redirection_array_for_j[l-1] + m) )) / inSphere.sinTheta);
p_theta_derivated_term->setReal(magnitude * cosSin[m].getReal());
p_theta_derivated_term->setImag(magnitude * cosSin[m].getImag());
......@@ -481,7 +481,7 @@ protected:
);*/
// Computation of {\partial Inner_m^m(r, theta, phi)}/{\partial theta}:
magnitude = (*p_spherical_harmonic_Inner_coefficients_array) * r_l * (m * inSphere.cosTheta * (*ptr_associated_Legendre_function_Array) / inSphere.sinTheta);
magnitude = (*p_spherical_harmonic_Inner_coefficients_array) * FReal(r_l) * (FReal(m) * inSphere.cosTheta * (*ptr_associated_Legendre_function_Array) / inSphere.sinTheta);
p_theta_derivated_term->setReal(magnitude * cosSin[m].getReal());
p_theta_derivated_term->setImag(magnitude * cosSin[m].getImag());
......@@ -540,17 +540,17 @@ protected:
childBox.setPositionFromMorton(idxChild,1);
const F3DPosition M2MVector (
father.getX() - (treeWidthAtLevel * (1 + (childBox.getX() * 2))),
father.getY() - (treeWidthAtLevel * (1 + (childBox.getY() * 2))),
father.getZ() - (treeWidthAtLevel * (1 + (childBox.getZ() * 2)))
father.getX() - (treeWidthAtLevel * FReal(1 + (childBox.getX() * 2))),
father.getY() - (treeWidthAtLevel * FReal(1 + (childBox.getY() * 2))),
father.getZ() - (treeWidthAtLevel * FReal(1 + (childBox.getZ() * 2)))
);
harmonicInner(positionTsmphere(M2MVector),this->transitionM2M[idxLevel][idxChild]);
const F3DPosition L2LVector (
(treeWidthAtLevel * (1 + (childBox.getX() * 2))) - father.getX(),
(treeWidthAtLevel * (1 + (childBox.getY() * 2))) - father.getY(),
(treeWidthAtLevel * (1 + (childBox.getZ() * 2))) - father.getZ()
(treeWidthAtLevel * FReal(1 + (childBox.getX() * 2))) - father.getX(),
(treeWidthAtLevel * FReal(1 + (childBox.getY() * 2))) - father.getY(),
(treeWidthAtLevel * FReal(1 + (childBox.getZ() * 2))) - father.getZ()
);
harmonicInner(positionTsmphere(L2LVector),this->transitionL2L[idxLevel][idxChild]);
......@@ -587,7 +587,7 @@ protected:
//printf("x=%ld \t y=%ld \t z=%ld\n",x,y,z);
if( ( x*x + y*y + z*z ) >= ( 3*FMB_Info_ws*FMB_Info_ws + 0.1 ) ){
const F3DPosition relativePos( x*treeWidthAtLevel , y*treeWidthAtLevel , z*treeWidthAtLevel );
const F3DPosition relativePos( FReal(x) * treeWidthAtLevel , FReal(y) * treeWidthAtLevel , FReal(z) * treeWidthAtLevel );
// Not blas so
//printf("transferM2L[%d][%d][%d][%d]\n", idxLevel, idxd1, idxd2, idxd3);
......@@ -1105,7 +1105,7 @@ public:
exp_term_aux.setReal( (p_Y_term->getReal() * p_local_exp_term->getReal()) - (p_Y_term->getImag() * p_local_exp_term->getImag()) );
exp_term_aux.setImag( (p_Y_term->getReal() * p_local_exp_term->getImag()) + (p_Y_term->getImag() * p_local_exp_term->getReal()) );
force_vector_in_local_base.setX( force_vector_in_local_base.getX() + j * exp_term_aux.getReal());
force_vector_in_local_base.setX( force_vector_in_local_base.getX() + FReal(j) * exp_term_aux.getReal());
// F_phi: k=0 => nothing to do for F_phi
// F_theta:
exp_term_aux.setReal( (p_Y_theta_derivated_term->getReal() * p_local_exp_term->getReal()) - (p_Y_theta_derivated_term->getImag() * p_local_exp_term->getImag()) );
......@@ -1140,9 +1140,9 @@ public:
exp_term_aux.setReal( (p_Y_term->getReal() * p_local_exp_term->getReal()) - (p_Y_term->getImag() * p_local_exp_term->getImag()) );
exp_term_aux.setImag( (p_Y_term->getReal() * p_local_exp_term->getImag()) + (p_Y_term->getImag() * p_local_exp_term->getReal()) );
force_vector_in_local_base.setX(force_vector_in_local_base.getX() + 2 * j * exp_term_aux.getReal());
force_vector_in_local_base.setX(force_vector_in_local_base.getX() + FReal(2 * j) * exp_term_aux.getReal());
// F_phi:
force_vector_in_local_base.setZ( force_vector_in_local_base.getZ() - 2 * k * exp_term_aux.getImag());
force_vector_in_local_base.setZ( force_vector_in_local_base.getZ() - FReal(2 * k) * exp_term_aux.getImag());
// F_theta:
/*printf("\t\t k = %d \t j = %d \t exp_term_aux real = %e imag = %e \n",
......@@ -1155,7 +1155,7 @@ public:
exp_term_aux.setReal( (p_Y_theta_derivated_term->getReal() * p_local_exp_term->getReal()) - (p_Y_theta_derivated_term->getImag() * p_local_exp_term->getImag()) );
exp_term_aux.setImag( (p_Y_theta_derivated_term->getReal() * p_local_exp_term->getImag()) + (p_Y_theta_derivated_term->getImag() * p_local_exp_term->getReal()) );
force_vector_in_local_base.setY(force_vector_in_local_base.getY() + 2 * exp_term_aux.getReal());
force_vector_in_local_base.setY(force_vector_in_local_base.getY() + FReal(2.0) * exp_term_aux.getReal());
/*printf("\t\t k = %d \t j = %d \t exp_term_aux real = %lf imag = %lf \n",
k,j,
......@@ -1177,21 +1177,21 @@ public:
// We want: - gradient(POTENTIAL_SIGN potential).
// The -(- 1.0) computing is not the most efficient programming ...
//#define FMB_TMP_SIGN -(POTENTIAL_SIGN 1.0)
force_vector_in_local_base.setX( force_vector_in_local_base.getX() * (-1.0) / spherical.r);
force_vector_in_local_base.setY( force_vector_in_local_base.getY() * (-1.0) / spherical.r);
force_vector_in_local_base.setZ( force_vector_in_local_base.getZ() * (-1.0) / (spherical.r * spherical.sinTheta));
force_vector_in_local_base.setX( force_vector_in_local_base.getX() * FReal(-1.0) / spherical.r);
force_vector_in_local_base.setY( force_vector_in_local_base.getY() * FReal(-1.0) / spherical.r);
force_vector_in_local_base.setZ( force_vector_in_local_base.getZ() * FReal(-1.0) / (spherical.r * spherical.sinTheta));
//#undef FMB_TMP_SIGN
/////////////////////////////////////////////////////////////////////
//spherical_position_Set_ph
//FMB_INLINE COORDINATES_T angle_Convert_in_MinusPi_Pi(COORDINATES_T a){
FReal ph = FMath::Fmod(spherical.phi, 2*FMath::FPi);
if (ph > M_PI) ph -= 2*FMath::FPi;
if (ph < -M_PI + FMath::Epsilon) ph += 2 * FMath::Epsilon;
FReal ph = FMath::Fmod(spherical.phi, FReal(2)*FMath::FPi);
if (ph > M_PI) ph -= FReal(2) * FMath::FPi;
if (ph < -M_PI + FMath::Epsilon) ph += FReal(2) * FMath::Epsilon;
//spherical_position_Set_th
FReal th = FMath::Fmod(FMath::ACos(spherical.cosTheta), 2*FMath::FPi);
FReal th = FMath::Fmod(FMath::ACos(spherical.cosTheta), FReal(2) * FMath::FPi);
if (th < 0.0) th += 2*FMath::FPi;
if (th > FMath::FPi){
th = 2*FMath::FPi - th;
......@@ -1366,7 +1366,7 @@ public:
const FReal dy = target.getPosition().getY() - source.getPosition().getY();
const FReal dz = target.getPosition().getZ() - source.getPosition().getZ();
FReal inv_square_distance = 1.0/ (dx*dx + dy*dy + dz*dz + FMB_Info_eps_soft_square);
FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz + FReal(FMB_Info_eps_soft_square));
FReal inv_distance = FMath::Sqrt(inv_square_distance);
inv_distance *= target.getPhysicalValue() * source.getPhysicalValue();
inv_square_distance *= inv_distance;
......@@ -1442,7 +1442,7 @@ public:
const FReal dy = target.getPosition().getY() - source.getPosition().getY();
const FReal dz = target.getPosition().getZ() - source.getPosition().getZ();
FReal inv_square_distance = 1.0/ (dx*dx + dy*dy + dz*dz + FMB_Info_eps_soft_square);
FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz + FReal(FMB_Info_eps_soft_square));
FReal inv_distance = FMath::Sqrt(inv_square_distance);
inv_distance *= target.getPhysicalValue() * source.getPhysicalValue();
inv_square_distance *= inv_distance;
......
Src/Fmb/FFmbKernelsBlockBlas.hpp
View file @ 1032efb7
......@@ -143,7 +143,7 @@ protected:
//std::cout << "sphereHarmoOuterCoef\n";
FReal factOuter = 1.0;
// in FMB code stoped at <= FMB_Info_M2L_P but this is not sufficient
for(int idxP = 0 ; idxP <= FMB_Info_M2L_P; factOuter *= (++idxP) ){
for(int idxP = 0 ; idxP <= FMB_Info_M2L_P; factOuter *= FReal(++idxP) ){
this->sphereHarmoOuterCoef[idxP] = factOuter;
//printf("spherical_harmonic_Outer_coefficients_array %.15e\n",this->sphereHarmoOuterCoef[idxP]);
//printf("fact_l %.15e\n",factOuter);
......@@ -155,9 +155,9 @@ protected:
FReal factInner = 1.0;
FReal powN1idxP = 1.0;
//std::cout << "sphereHarmoInnerCoef\n";
for(int idxP = 0 ; idxP <= this->FMB_Info_M2L_P ; factInner *= (++idxP), powN1idxP = -powN1idxP){
for(int idxMP = 0, fact_l_m = factInner; idxMP <= idxP ; fact_l_m *= idxP+(++idxMP), ++currentInner){
*currentInner = powN1idxP / fact_l_m;
for(int idxP = 0 ; idxP <= this->FMB_Info_M2L_P ; factInner *= FReal(++idxP), powN1idxP = -powN1idxP){
for(int idxMP = 0, fact_l_m = int(factInner); idxMP <= idxP ; fact_l_m *= idxP+(++idxMP), ++currentInner){
*currentInner = powN1idxP / FReal(fact_l_m);
//std::cout << (*currentInner) << "\n";
}
}
......@@ -348,8 +348,8 @@ protected:
for(int idxl = 2; idxl <= lmax ; ++idxl ){
for( int idxm = 0; idxm <= idxl - 2 ; ++idxm , ++idxCurrent , ++idxCurrent1m , ++idxCurrent2m ){
// Compute P_l^m, l >= m+2, using (1) and store it into results_array:
outResults[idxCurrent] = (inCosTheta * ( 2 * idxl - 1 ) * outResults[idxCurrent1m] - ( idxl + idxm - 1 )
* outResults[idxCurrent2m] ) / ( idxl - idxm );
outResults[idxCurrent] = (inCosTheta * FReal( 2 * idxl - 1 ) * outResults[idxCurrent1m] - FReal( idxl + idxm - 1 )
* outResults[idxCurrent2m] ) / FReal( idxl - idxm );
/*printf("\tres (%.18e) = inCosTheta (%.18e) idxl (%d) outResults 1m (%.18e) idxm (%d) outResults 2m (%.18e)\n",
outResults[idxCurrent],inCosTheta,idxl,outResults[idxCurrent1m],idxm,outResults[idxCurrent2m]);
printf("\t\t temp1 = %.18e temp2 = %.18e temp3 = %.18e\n",
......@@ -360,12 +360,12 @@ protected:
// p_results_array_l_minus_1_m now points on P_{l-1}^{l-1}
// Compute P_l^{l-1} using (3) and store it into ptrResults:
outResults[idxCurrent++] = inCosTheta * ( 2 * idxl - 1 ) * outResults[idxCurrent1m];
outResults[idxCurrent++] = inCosTheta * FReal( 2 * idxl - 1 ) * outResults[idxCurrent1m];
// Compute P_l^l using (2 bis) and store it into results_array:
outResults[idxCurrent++] = fact * invSinTheta * outResults[idxCurrent1m];
fact += 2.0;
fact += FReal(2.0);
++idxCurrent1m;
}
......@@ -432,7 +432,7 @@ protected:
for(int idxl = 0, idxlMod4 = 0; idxl <= FMB_Info_M2L_P ; ++idxl, ++idxlMod4){
if(idxlMod4 == 4) idxlMod4 = 0;
const FReal angle = idxl * inSphere.phi + this->PiArrayOuter[idxlMod4];
const FReal angle = FReal(idxl) * inSphere.phi + this->PiArrayOuter[idxlMod4];
cosSin[idxl].setReal( FMath::Sin(angle + FMath::FPiDiv2) );
cosSin[idxl].setImag( FMath::Sin(angle) );
......@@ -514,7 +514,7 @@ protected:
// Initialization of precomputed_cos_and_sin_array:
for(int idxm = 0 , idxmMod4 = 0; idxm <= FMB_Info_P ; ++idxm, ++idxmMod4){
if(idxmMod4 == 4) idxmMod4 = 0;
const FReal angle = idxm*inSphere.phi + PiArrayInner[idxmMod4];
const FReal angle = FReal(idxm) *inSphere.phi + PiArrayInner[idxmMod4];
cosSin[idxm].setReal(FMath::Sin(angle + FMath::FPiDiv2));
cosSin[idxm].setImag(FMath::Sin(angle));
......@@ -561,8 +561,8 @@ protected:
);*/
// Computation of {\partial Inner_l^m(r, theta, phi)}/{\partial theta}:
magnitude = (*p_spherical_harmonic_Inner_coefficients_array) * r_l * ((l*inSphere.cosTheta*(*ptr_associated_Legendre_function_Array)
- (l+m)*(*(start_ptr_associated_Legendre_function_Array + expansion_Redirection_array_for_j[l-1] + m))) / inSphere.sinTheta);
magnitude = (*p_spherical_harmonic_Inner_coefficients_array) * FReal(r_l) * ((FReal(l)*inSphere.cosTheta*(*ptr_associated_Legendre_function_Array)
- FReal(l+m)*(*(start_ptr_associated_Legendre_function_Array + expansion_Redirection_array_for_j[l-1] + m))) / inSphere.sinTheta);
p_theta_derivated_term->setReal(magnitude * cosSin[m].getReal());
p_theta_derivated_term->setImag(magnitude * cosSin[m].getImag());
......@@ -590,7 +590,7 @@ protected:
);*/
// Computation of {\partial Inner_m^m(r, theta, phi)}/{\partial theta}:
magnitude = (*p_spherical_harmonic_Inner_coefficients_array) * r_l * (m * inSphere.cosTheta * (*ptr_associated_Legendre_function_Array) / inSphere.sinTheta);
magnitude = (*p_spherical_harmonic_Inner_coefficients_array) * FReal(r_l) * (FReal(m) * inSphere.cosTheta * (*ptr_associated_Legendre_function_Array) / inSphere.sinTheta);
p_theta_derivated_term->setReal(magnitude * cosSin[m].getReal());
p_theta_derivated_term->setImag(magnitude * cosSin[m].getImag());
......@@ -625,9 +625,9 @@ protected:
childBox.setPositionFromMorton(idxChild,1);
const F3DPosition M2MVector (
father.getX() - (treeWidthAtLevel * (1 + (childBox.getX() * 2))),
father.getY() - (treeWidthAtLevel * (1 + (childBox.getY() * 2))),
father.getZ() - (treeWidthAtLevel * (1 + (childBox.getZ() * 2)))
father.getX() - (treeWidthAtLevel * FReal(1 + (childBox.getX() * 2))),
father.getY() - (treeWidthAtLevel * FReal(1 + (childBox.getY() * 2))),
father.getZ() - (treeWidthAtLevel * FReal(1 + (childBox.getZ() * 2)))
);
harmonicInner(positionTsmphere(M2MVector),this->transitionM2M[idxLevel][idxChild]);
......@@ -641,9 +641,9 @@ protected:
}*/
const F3DPosition L2LVector (
(treeWidthAtLevel * (1 + (childBox.getX() * 2))) - father.getX(),
(treeWidthAtLevel * (1 + (childBox.getY() * 2))) - father.getY(),
(treeWidthAtLevel * (1 + (childBox.getZ() * 2))) - father.getZ()
(treeWidthAtLevel * FReal(1 + (childBox.getX() * 2))) - father.getX(),
(treeWidthAtLevel * FReal(1 + (childBox.getY() * 2))) - father.getY(),
(treeWidthAtLevel * FReal(1 + (childBox.getZ() * 2))) - father.getZ()
);
harmonicInner(positionTsmphere(L2LVector),this->transitionL2L[idxLevel][idxChild]);
......@@ -683,8 +683,8 @@ protected:
k = NN-n-m;
if (k < 0){
const int pow_of_minus_1 = ((k%2) ? -1 : 1);
data->setReal( pow_of_minus_1 * (transfer_exp + expansion_Redirection_array_for_j[MM+NN] - k)->getReal());
data->setImag((-pow_of_minus_1) * (transfer_exp + expansion_Redirection_array_for_j[MM+NN] - k)->getImag());
data->setReal( FReal(pow_of_minus_1) * (transfer_exp + expansion_Redirection_array_for_j[MM+NN] - k)->getReal());
data->setImag( FReal(-pow_of_minus_1) * (transfer_exp + expansion_Redirection_array_for_j[MM+NN] - k)->getImag());
}
else {
*data = *(transfer_exp + expansion_Redirection_array_for_j[MM+NN] + k);
......@@ -767,14 +767,14 @@ protected:
for( int idxd2 = 0; idxd2 < this->size1Dim ; ++idxd2 ){
for( int idxd3 = 0; idxd3 < this->size1Dim ; ++idxd3 ){
const long x = idxd1 - this->halphSize1Dim;
const long y = idxd2 - this->halphSize1Dim;
const long z = idxd3 - this->halphSize1Dim;
const FReal x = FReal(idxd1 - this->halphSize1Dim);
const FReal y = FReal(idxd2 - this->halphSize1Dim);
const FReal z = FReal(idxd3 - this->halphSize1Dim);
//printf("x=%ld \t y=%ld \t z=%ld\n",x,y,z);
if( ( x*x + y*y + z*z ) >= ( 3*FMB_Info_ws*FMB_Info_ws + 0.1 ) ){
const F3DPosition relativePos( x*treeWidthAtLevel , y*treeWidthAtLevel , z*treeWidthAtLevel );
const F3DPosition relativePos( x*treeWidthAtLevel , y*treeWidthAtLevel , z*treeWidthAtLevel );
//printf("blas\n");
......@@ -1035,7 +1035,7 @@ public:
* Position in 'nexp': j*(j+1) + k
*/
int jj_plus_1 = j*(j+1);
int half_jj_plus_1 = (int) jj_plus_1*0.5;
int half_jj_plus_1 = int(FReal(jj_plus_1)*0.5);
p_exp = exp + half_jj_plus_1 + j;
p_nexp = exp + jj_plus_1 + j;
......@@ -1415,7 +1415,7 @@ public:
exp_term_aux.setReal( (p_Y_term->getReal() * p_local_exp_term->getReal()) - (p_Y_term->getImag() * p_local_exp_term->getImag()) );
exp_term_aux.setImag( (p_Y_term->getReal() * p_local_exp_term->getImag()) + (p_Y_term->getImag() * p_local_exp_term->getReal()) );
force_vector_in_local_base.setX( force_vector_in_local_base.getX() + j * exp_term_aux.getReal());
force_vector_in_local_base.setX( force_vector_in_local_base.getX() + FReal(j) * exp_term_aux.getReal());
// F_phi: k=0 => nothing to do for F_phi
// F_theta:
exp_term_aux.setReal( (p_Y_theta_derivated_term->getReal() * p_local_exp_term->getReal()) - (p_Y_theta_derivated_term->getImag() * p_local_exp_term->getImag()) );
......@@ -1450,9 +1450,9 @@ public:
exp_term_aux.setReal( (p_Y_term->getReal() * p_local_exp_term->getReal()) - (p_Y_term->getImag() * p_local_exp_term->getImag()) );
exp_term_aux.setImag( (p_Y_term->getReal() * p_local_exp_term->getImag()) + (p_Y_term->getImag() * p_local_exp_term->getReal()) );
force_vector_in_local_base.setX(force_vector_in_local_base.getX() + 2 * j * exp_term_aux.getReal());
force_vector_in_local_base.setX(force_vector_in_local_base.getX() + FReal(2 * j) * exp_term_aux.getReal());
// F_phi:
force_vector_in_local_base.setZ( force_vector_in_local_base.getZ() - 2 * k * exp_term_aux.getImag());
force_vector_in_local_base.setZ( force_vector_in_local_base.getZ() - FReal(2 * k) * exp_term_aux.getImag());
// F_theta:
/*printf("\t\t k = %d \t j = %d \t exp_term_aux real = %.15e imag = %.15e \n",
......@@ -1487,9 +1487,9 @@ public:
// We want: - gradient(POTENTIAL_SIGN potential).
// The -(- 1.0) computing is not the most efficient programming ...
//#define FMB_TMP_SIGN -(POTENTIAL_SIGN 1.0)