Remove direct computation in utest

Examples/DirectComputation.cpp

+// ===================================================================================
+// ===================================================================================
+// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas, Matthias Messner
+// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
+// This software is a computer program whose purpose is to compute the FMM.
+//
+// This software is governed by the CeCILL-C and LGPL licenses and
+// abiding by the rules of distribution of free software.  
+// 
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public and CeCILL-C Licenses for more details.
+// "http://www.cecill.info". 
+// "http://www.gnu.org/licenses".
+// ===================================================================================
+
+#include <iostream>
+#include <iomanip>
+
+#include <cstdio>
+#include <cstdlib>
+#include <cmath>
+#include <algorithm>
+
+#include  "ScalFmmConfig.h"
+#include "Utils/FTic.hpp"
+//#include "Utils/FMath.hpp"
+#include "Utils/FParameters.hpp"
+//#include "Utils/FIOVtk.hpp"
+
+//#include "Containers/FVector.hpp"
+#include "Files/FFmaGenericLoader.hpp"
+#include "Kernels/P2P/FP2P.hpp"
+
+
+
+
+// Simply create particles and try the kernels
+int main(int argc, char ** argv){
+	//
+	///////////////////////What we do/////////////////////////////
+	if( FParameters::existParameter(argc, argv, "-help" ) || argc < 4){
+		std::cout << ">> This executable has to be used to compute  interaction either for periodic or non periodic system.\n";
+		std::cout << ">> Example -fin filenameIN.{fma or bfma)     -fout filenameOUT{fma or bfma) \n";
+		std::cout << ">> Default input file : ../Data/unitCubeXYZQ20k.fma\n";
+		std::cout << " Options " << std::endl;
+		std::cout << "     -verbose : print index x y z Q V fx fy fz " << std::endl;
+		std::cout << "     -fin filename. Extension specifies if the file is binary or not. " << std::endl;
+		std::cout << "                            Only our FMA (.bma, .bfma) is allowed " << std::endl;
+		std::cout << "     -fout filenameOUT   output file  with extension (default output.bfma)" << std::endl;
+		exit(-1);
+
+	}
+
+	//////////////////////////////////////////////////////////////
+
+	const std::string filenameIn(FParameters::getStr(argc,argv,"-fin", "../Data/unitCubeXYZQ20k.fma"));
+	const std::string filenameOut(FParameters::getStr(argc,argv,"-fout", "output.bfma"));
+	//
+	FTic counter;
+
+	// -----------------------------------------------------
+	//  LOADER
+	//  -----------------------------------------------------
+	// ---------------------------------------------------------------------------------
+	// Read  particles in the Octree
+	// ---------------------------------------------------------------------------------
+	std::cout << "Opening : " << filenameIn << "\n";
+	//
+	FFmaGenericLoader loader(filenameIn);
+	//
+	int nbParticles = static_cast<int>(loader.getNumberOfParticles());
+	std::cout << "Read " << nbParticles << " particles ..." << std::endl;
+	double BoxWith=loader.getBoxWidth();
+	FPoint Centre(loader.getCenterOfBox().getX(), loader.getCenterOfBox().getY() , loader.getCenterOfBox().getZ());
+	std::cout << "\tWidth : " <<BoxWith << " \t center x : " << loader.getCenterOfBox().getX()
+	    																	<< " y : " << loader.getCenterOfBox().getY() << " z : " << loader.getCenterOfBox().getZ() << std::endl;
+
+	counter.tic();
+	//
+	FmaRParticle *  particles = new FmaRParticle[nbParticles];
+	memset(particles, 0, sizeof(FmaRParticle) * nbParticles) ;
+	//
+	double totalCharge = 0.0;
+	//
+	int nbDataToRead = particles[0].getReadDataNumber();
+	for(int idx = 0 ; idx<nbParticles ; ++idx){
+		//
+		loader.fillParticle(&particles[idx].position, &particles[idx].physicalValue);
+		//	loader.fillParticle(particles[idx].getPtrFirstData(), nbDataToRead);    // OK
+		//  loader.fillParticle(particles[idx]); // OK
+		std::cout << idx <<"  "<<  particles[idx].position << " "<<particles[idx].physicalValue << " "<<particles[idx].potential
+				<<"  " << particles[idx].forces[0]<<"  " <<particles[idx].forces[1]<<"  " <<particles[idx].forces[2]<<"  " <<std::endl;
+		//
+		totalCharge += particles[idx].physicalValue ;
+	}
+
+	counter.tac();
+
+	std::cout << std::endl;
+	std::cout << "Total Charge         = "<< totalCharge <<std::endl;
+	std::cout << std::endl;
+
+	std::cout << "Done  " << "(@ reading Particles  " << counter.elapsed() << " s)." << std::endl;
+	//
+	// ----------------------------------------------------------------------------------------------------------
+	//                                   COMPUTATION
+	// ----------------------------------------------------------------------------------------------------------
+	FReal denergy = 0.0;
+	//
+	//  computation
+	//
+	{
+		printf("Compute :\n");
+		counter.tic();
+#pragma omp parallel shared(nbParticles, particles,denergy)
+		{
+#pragma omp for
+			for(int idxTarget = 0 ; idxTarget < nbParticles ; ++idxTarget){
+				//
+				// compute with all other except itself
+				//
+				// Compute force and potential between  particles[idxTarget] and particles inside the box
+				//
+				for(int idxOther = 0; idxOther < nbParticles ; ++idxOther){
+					if( idxOther != idxTarget ){
+						FP2P::NonMutualParticles(
+								particles[idxOther].position.getX(), particles[idxOther].position.getY(),
+								particles[idxOther].position.getZ(),particles[idxOther].physicalValue,
+								particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
+								particles[idxTarget].position.getZ(),particles[idxTarget].physicalValue,
+								&particles[idxTarget].forces[0],&particles[idxTarget].forces[1],
+								&particles[idxTarget].forces[2],&particles[idxTarget].potential);
+					}
+				}
+			} // end for
+			// Compute the energy
+#pragma omp  for reduction(+:denergy)
+			for(int idx = 0 ; idx < nbParticles ; ++idx){
+				denergy +=  particles[idx].potential*particles[idx].physicalValue ;
+			}
+		} // end pragma parallel
+		//
+		denergy *= 0.5 ;
+		counter.tac();
+		//
+		printf("Energy =   %.14e\n",denergy);
+		//
+		std::cout << "Done  " << "(@ Direct computation done = " << counter.elapsed() << " s)." << std::endl;
+		std::cout << "\n"<< "END  "
+				<< "-------------------------------------------------------------------------"
+				<< std::endl << std::endl ;
+	} // END
+
+	//
+	// ----------------------------------------------------------------
+	//  Save  computation in binary format
+	//
+	//
+
+	std::cout << "Generate " << filenameOut <<"  for output file" << std::endl;
+	//
+	std::cout << " nbParticles: " << nbParticles <<"  " << sizeof(nbParticles) <<std::endl;
+	std::cout << " denergy: " << denergy <<"  " << sizeof(denergy) <<std::endl;
+	std::cout << " Box size: " << loader.getBoxWidth() << "  " << sizeof(loader.getBoxWidth())<<std::endl;
+	//
+	FFmaGenericWriter writer(filenameOut) ;
+	writer.writeHeader(Centre,BoxWith, nbParticles,*particles) ;
+	writer.writeArrayOfParticles(particles, nbParticles);
+	//
+	// end generate
+	// -----------------------------------------------------
+	//
+	if(FParameters::existParameter(argc, argv, "-verbose")){
+		denergy = 0 ;
+		for(int idx = 0 ; idx < nbParticles ; ++idx){
+			std::cout << ">> index " << idx << std::endl;
+			std::cout << " x   " << particles[idx].position.getX() << " y  " << particles[idx].position.getY() << " z  " << particles[idx].position.getZ() << std::endl;
+			std::cout << " Q   " << particles[idx].physicalValue   << " V  " << particles[idx].potential << std::endl;
+			std::cout << " fx  " << particles[idx].forces[0]       << " fy " << particles[idx].forces[1]       << " fz " << particles[idx].forces[2] << std::endl;
+			std::cout << "\n";
+			denergy +=  particles[idx].potential*particles[idx].physicalValue ;
+		}
+	}
+	std::cout << " ENERGY " << denergy << std::endl;
+	//
+	delete[] particles;
+	return 0;
+}
+
+
+

Src/Files/FFmaGenericLoader.hpp

@@ -42,7 +42,7 @@
 class FmaBasicParticle {
 public:
 	FPoint position;            ///< position of the particle
-//	FReal  x,y,z;                 ///< position of the particle
+	//	FReal  x,y,z;                 ///< position of the particle
 	FReal  physicalValue;    ///<  its physical value
 	/**
 	 *  return a pointer on the first value of the structure
@@ -128,6 +128,10 @@ public:
 	int getReadDataNumber()
 	{ return 8;}
 };
+
+typedef FmaBasicParticle FmaR4W4Particle ;
+typedef FmaRParticle       FmaR4W8Particle ;
+typedef FmaParticle         FmaR8W8Particle ;
 //
 //! \class  FFmaGenericLoader
 //!
@@ -345,11 +349,11 @@ public:
 	 */
 	template <class dataPart>
 	void fillParticle(dataPart &dataToRead){
-		int otherDataToRead = typeData[1] - dataToRead.getReadDataNumber() ;
-	if(binaryFile){
+		int otherDataRead = typeData[1] - dataToRead.getReadDataNumber() ;
+		if(binaryFile){
 			file->read((char*)(dataToRead.getPtrFirstData()), sizeof(FReal)*(dataToRead.getReadDataNumber()));
-			if( otherDataToRead > 0){
-				file->read((char*)(this->tmpVal), sizeof(FReal)*(otherDataToRead));
+			if( otherDataRead > 0){
+				file->read((char*)(this->tmpVal), sizeof(FReal)*(otherDataRead));
 			}
 		}
 		else{
@@ -358,9 +362,9 @@ public:
 				(*this->file)  >>*val;
 				++val;
 			}
-			if( otherDataToRead > 0){
+			if( otherDataRead > 0){
 				FReal x;
-				for (int i = 0 ; i <otherDataToRead ;++i){
+				for (int i = 0 ; i <otherDataRead ;++i){
 					(*this->file)  >>x;
 				}
 			}
@@ -374,8 +378,8 @@ public:
 	 */
 	template <class dataPart>
 	void fillParticle(dataPart *dataToRead, const int N){
-		int otherDataToRead = typeData[1] - (*dataToRead).getReadDataNumber() ;
-		if(binaryFile && otherDataToRead == 0 ){
+		int otherDataRead = typeData[1] - (*dataToRead).getReadDataNumber() ;
+		if(binaryFile && otherDataRead == 0 ){
 			file->read((char*)((*dataToRead).getPtrFirstData()), sizeof(FReal)*(N*(*dataToRead).getReadDataNumber()));
 		}
 		else {
@@ -638,7 +642,8 @@ public:
 		}
 		else{ // ASCII part
 			const int ndata = dataToWrite[0].getWriteDataNumber();
-			std::cout << "typeData "<< sizeof(FReal) << " "<<ndata << std::endl;
+//			std::cout << "typeData "<< sizeof(FReal) << " "<<ndata << std::endl;
+			this->file->precision(10);
 
 			for (int i = 0 ; i <N ; ++i){
 				const FReal * val = dataToWrite[i].getPtrFirstData() ;
@@ -673,6 +678,7 @@ public:
 			file->write((char*)(dataToWrite), N*nbData*sizeof(FReal));
 		}
 		else{
+			this->file->precision(10);
 			//			std::cout << "N "<< N << " nbData "<< nbData<<std::endl;
 			//			exit(-1);
 			int k = 0;
@@ -691,6 +697,7 @@ public:
 private:
 	void writerAscciHeader( const FPoint &centerOfBox,const FReal &boxWidth,
 			const FSize &nbParticles, const unsigned int *typeFReal) {
+		this->file->precision(10);
 		(*this->file) << typeFReal[0] <<"   "<<typeFReal[1]<<std::endl;
 		(*this->file) << nbParticles << "   "<<  boxWidth << "   "
 				<<  centerOfBox.getX()  << "  " << centerOfBox.getY() << " "<<centerOfBox.getZ()

UTests/utestChebyshev.cpp

+// ===================================================================================
+// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas, Matthias Messner
+// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
+// This software is a computer program whose purpose is to compute the FMM.
+//
+// This software is governed by the CeCILL-C and LGPL licenses and
+// abiding by the rules of distribution of free software.  
+// 
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public and CeCILL-C Licenses for more details.
+// "http://www.cecill.info". 
+// "http://www.gnu.org/licenses".
+// ===================================================================================
+
+// ==== CMAKE =====
+// @FUSE_BLAS
+// ================
+
+#include "Utils/FGlobal.hpp"
+
+#include "Containers/FOctree.hpp"
+
+#include "Files/FFmaGenericLoader.hpp"
+
+#include "Core/FFmmAlgorithmThread.hpp"
+#include "Core/FFmmAlgorithm.hpp"
+
+#include "FUTester.hpp"
+
+#include "Components/FSimpleLeaf.hpp"
+
+
+#include "Kernels/Chebyshev/FChebCell.hpp"
+#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
+#include "Kernels/Chebyshev/FChebKernel.hpp"
+#include "Kernels/Chebyshev/FChebSymKernel.hpp"
+
+#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
+/*
+  In this test we compare the Chebyschev fmm results and the direct results.
+ */
+
+
+/** the test class
+ *
+ */
+class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
+
+	///////////////////////////////////////////////////////////
+	// The tests!
+	///////////////////////////////////////////////////////////
+
+	template <class CellClass, class ContainerClass, class KernelClass, class MatrixKernelClass,
+	class LeafClass, class OctreeClass, class FmmClass>
+	void RunTest()	{
+		//
+		// Load particles
+		//
+		if(sizeof(FReal) == sizeof(float) ) {
+			std::cerr << "No input data available for Float "<< std::endl;
+			exit(-1);
+		}
+		const char* const filename = (sizeof(FReal) == sizeof(float))?
+				"../Data/UTest/DirectFloatbfma":
+				"../Data/UTest/DirectDouble.bfma";
+		//checkDirect.bfma
+		FFmaGenericLoader loader(filename);
+		Print("Number of particles:");
+		Print(loader.getNumberOfParticles());
+
+		const int NbLevels        = 4;
+		const int SizeSubLevels = 2;
+
+		// Create octree
+
+		FSize nbParticles = loader.getNumberOfParticles() ;
+		FmaR8W8Particle* const particles = new FmaR8W8Particle[nbParticles];
+
+		loader.fillParticle(particles,nbParticles);
+         //
+		//   Insert particle in the tree
+		//
+		OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+		for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
+			tree.insert(particles[idxPart].position , idxPart, particles[idxPart].physicalValue );
+		}
+		//
+		/////////////////////////////////////////////////////////////////////////////////////////////////
+		// Run FMM computation
+		/////////////////////////////////////////////////////////////////////////////////////////////////
+		Print("Fmm...");
+		KernelClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+		FmmClass algo(&tree,&kernels);
+		algo.execute();
+		//0
+		FReal energy= 0.0 , energyD = 0.0 ;
+		/////////////////////////////////////////////////////////////////////////////////////////////////
+		// Compute direct energy
+		/////////////////////////////////////////////////////////////////////////////////////////////////
+
+		for(int idx = 0 ; idx <  loader.getNumberOfParticles()  ; ++idx){
+			energyD +=  particles[idx].potential*particles[idx].physicalValue ;
+		}
+		/////////////////////////////////////////////////////////////////////////////////////////////////
+		// Compare
+		/////////////////////////////////////////////////////////////////////////////////////////////////
+		Print("Compute Diff...");
+		FMath::FAccurater potentialDiff;
+		FMath::FAccurater fx, fy, fz;
+		{ // Check that each particle has been summed with all other
+
+			tree.forEachLeaf([&](LeafClass* leaf){
+				const FReal*const potentials        = leaf->getTargets()->getPotentials();
+				const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
+				const FReal*const forcesX            = leaf->getTargets()->getForcesX();
+				const FReal*const forcesY            = leaf->getTargets()->getForcesY();
+				const FReal*const forcesZ            = leaf->getTargets()->getForcesZ();
+				const int nbParticlesInLeaf           = leaf->getTargets()->getNbParticles();
+				const FVector<int>& indexes      = leaf->getTargets()->getIndexes();
+
+				for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
+					const int indexPartOrig = indexes[idxPart];
+					potentialDiff.add(particles[indexPartOrig].potential,potentials[idxPart]);
+					fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]);
+					fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
+					fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
+					energy   += potentials[idxPart]*physicalValues[idxPart];
+				}
+			});
+		}
+
+		delete[] particles;
+
+		// Print for information
+
+		Print("Potential diff is = ");
+		printf("         Pot L2Norm     %e\n",potentialDiff.getL2Norm());
+		printf("         Pot RL2Norm   %e\n",potentialDiff.getRelativeL2Norm());
+		printf("         Pot RMSError   %e\n",potentialDiff.getRMSError());
+		Print("Fx diff is = ");
+		printf("         Fx L2Norm     %e\n",fx.getL2Norm());
+		printf("         Fx RL2Norm   %e\n",fx.getRelativeL2Norm());
+		printf("         Fx RMSError   %e\n",fx.getRMSError());
+		Print("Fy diff is = ");
+		printf("        Fy L2Norm     %e\n",fy.getL2Norm());
+		printf("        Fy RL2Norm   %e\n",fy.getRelativeL2Norm());
+		printf("        Fy RMSError   %e\n",fy.getRMSError());
+		Print("Fz diff is = ");
+		printf("        Fz L2Norm     %e\n",fz.getL2Norm());
+		printf("        Fz RL2Norm   %e\n",fz.getRelativeL2Norm());
+		printf("        Fz RMSError   %e\n",fz.getRMSError());
+		FReal L2error = (fx.getRelativeL2Norm()*fx.getRelativeL2Norm() + fy.getRelativeL2Norm()*fy.getRelativeL2Norm()  + fz.getRelativeL2Norm() *fz.getRelativeL2Norm()  );
+		printf(" Total L2 Force Error= %e\n",FMath::Sqrt(L2error)) ;
+		printf("  Energy Error  =   %.12e\n",FMath::Abs(energy-energyD));
+		printf("  Energy FMM    =   %.12e\n",FMath::Abs(energy));
+		printf("  Energy DIRECT =   %.12e\n",FMath::Abs(energyD));
+
+		// Assert
+		const FReal MaximumDiffPotential = FReal(9e-3);
+		const FReal MaximumDiffForces     = FReal(9e-2);
+
+		Print("Test1 - Error Relative L2 norm Potential ");
+		uassert(potentialDiff.getRelativeL2Norm() < MaximumDiffPotential);    //1
+		Print("Test2 - Error RMS L2 norm Potential ");
+		uassert(potentialDiff.getRMSError() < MaximumDiffPotential);  //2
+		Print("Test3 - Error Relative L2 norm FX ");
+		uassert(fx.getRelativeL2Norm()  < MaximumDiffForces);                       //3
+		Print("Test4 - Error RMS L2 norm FX ");
+		uassert(fx.getRMSError() < MaximumDiffForces);                      //4
+		Print("Test5 - Error Relative L2 norm FY ");
+		uassert(fy.getRelativeL2Norm()  < MaximumDiffForces);                       //5
+		Print("Test6 - Error RMS L2 norm FY ");
+		uassert(fy.getRMSError() < MaximumDiffForces);                      //6
+		Print("Test7 - Error Relative L2 norm FZ ");
+		uassert(fz.getRelativeL2Norm()  < MaximumDiffForces);                      //8
+		Print("Test8 - Error RMS L2 norm FZ ");
+		uassert(fz.getRMSError() < MaximumDiffForces);                                           //8
+		Print("Test9 - Error Relative L2 norm F ");
+		uassert(L2error              < MaximumDiffForces);                                            //9   Total Force
+		Print("Test10 - Relative error Energy ");
+		uassert(FMath::Abs(energy-energyD) /energyD< MaximumDiffPotential);                     //10  Total Energy
+
+
+	}
+
+	/** If memstas is running print the memory used */
+	void PostTest() {
+		if( FMemStats::controler.isUsed() ){
+			std::cout << "Memory used at the end " << FMemStats::controler.getCurrentAllocated()
+										<< " Bytes (" << FMemStats::controler.getCurrentAllocatedMB() << "MB)\n";
+			std::cout << "Max memory used " << FMemStats::controler.getMaxAllocated()
+										<< " Bytes (" << FMemStats::controler.getMaxAllocatedMB() << "MB)\n";
+			std::cout << "Total memory used " << FMemStats::controler.getTotalAllocated()
+										<< " Bytes (" << FMemStats::controler.getTotalAllocatedMB() << "MB)\n";
+		}
+	}
+
+
+	///////////////////////////////////////////////////////////
+	// Set the tests!
+	///////////////////////////////////////////////////////////
+
+
+	/** TestChebKernel */
+	void TestChebKernel(){
+		const unsigned int ORDER = 6;
+		typedef FP2PParticleContainerIndexed<> ContainerClass;
+		typedef FSimpleLeaf<ContainerClass> LeafClass;
+		typedef FInterpMatrixKernelR MatrixKernelClass;
+		typedef FChebCell<ORDER> CellClass;
+		typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
+		typedef FChebKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
+		typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
+		// run test
+		RunTest<CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass>();
+	}
+
+	/** TestChebSymKernel */
+	void TestChebSymKernel(){
+		const unsigned int ORDER = 6;
+		typedef FP2PParticleContainerIndexed<> ContainerClass;
+		typedef FSimpleLeaf<ContainerClass> LeafClass;
+		typedef FInterpMatrixKernelR MatrixKernelClass;
+		typedef FChebCell<ORDER> CellClass;
+		typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
+		typedef FChebSymKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
+		typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
+		// run test
+		RunTest<CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass>();
+	}
+
+
+
+	///////////////////////////////////////////////////////////
+	// Set the tests!
+	///////////////////////////////////////////////////////////
+
+	/** set test */
+	void SetTests(){
+		AddTest(&TestChebyshevDirect::TestChebKernel,"Test Chebyshev Kernel with one big SVD");
+		AddTest(&TestChebyshevDirect::TestChebSymKernel,"Test Chebyshev Kernel with 16 small SVDs and symmetries");
+	}
+};
+
+
+// You must do this
+TestClass(TestChebyshevDirect)
+
+
+
+