Small cleanup

CMakeLists.txt

@@ -76,7 +76,7 @@ ELSE()
   ENDIF()
 ENDIF()
 
-SET(CMAKE_CXX_FLAGS  ${CMAKE_CXX_FLAGS} " -fPIC -O1 --param=ssp-buffer-size=4 -Wformat -Wformat-security -Werror=format-security -Wall -std=c++0x")
+SET(CMAKE_CXX_FLAGS  ${CMAKE_CXX_FLAGS} " -fPIC -O0 --param=ssp-buffer-size=4 -Wformat -Wformat-security -Werror=format-security -Wall -std=c++0x")
 
 IF (DEBUG)
   SET(CMAKE_VERBOSE_MAKEFILE 1)

examples/qm/NDimTableTest.cpp

@@ -1,3 +1,4 @@
+#define DEBUG_NDIMTABLE
 
 #define YADE_FFTW3
 #include <iostream>
@@ -84,19 +85,30 @@ int main(void){
 	parts.push_back(&T2);
 	NDimTable<float> T3(parts);
 	std::cout << "constructed, print now ---- checking TENSOR PRODUCT\n";
+	double minRealT3(T3.at(0,0,0)),maxRealT3(T3.at(0,0,0)),sumAllT3(0);
 	for(int i=0 ; i<3; i++) {
 		for(int j=0 ; j<4 ; j++) {
 			for(int k=0 ; k<5 ; k++) {
 				std::cout << T3.at(i,j,k) << " (" << (100*i+10*j +1)*(1000*k+1000) << ")   " ;
 				assert  (    T3.at(i,j,k) ==         (100*i+10*j +1)*(1000*k+1000) );
+				minRealT3 = std::min(minRealT3,  (double)    (100*i+10*j +1)*(1000*k+1000)    );
+				maxRealT3 = std::max(maxRealT3,  (double)    (100*i+10*j +1)*(1000*k+1000)    );
+				sumAllT3  +=                                (100*i+10*j +1)*(1000*k+1000)     ;
 				//std::cout << T1.at(i,j) << " ";
 			};
 			std::cout << "\n";
 		};
 		std::cout << "\n";
 	};
+	std::cerr.precision(20);
+	std::cerr << "checking minReal.T3: " << minRealT3  << "  " <<  T3.minReal() << "\n";
+	std::cerr << "checking maxReal.T3: " << maxRealT3  << "  " <<  T3.maxReal() << "\n";
+	std::cerr << "checking sumAll.T3 : " << sumAllT3   << "  " <<  T3.sumAll()  << "\n";
+	assert(minRealT3 == T3.minReal());
+	assert(maxRealT3 == T3.maxReal());
+	assert(sumAllT3  == T3.sumAll());
 
-	std::cout << "checking print() call\n\n";
+	std::cout << "checking print() call on T3\n\n";
 	T3.print();
 
 	std::cout << "T4\n";
@@ -246,8 +258,8 @@ int main(void){
 	in.becomesIFFT(out);
 	in.print();
 	out.print();
-	std::cout << "\n" << out.max();
-	std::cout << "\n" << out.min();
+	std::cout << "last printed tensor, out.maxReal() is " << out.maxReal() << "\n";
+	std::cout << "last printed tensor, out.minReal() is " << out.minReal() << "\n";
 // debug output:
 // FIXME, FIXME - add this to yade --check or test.
 // FIXME, FIXME - this is important, because various FFT libraries divide by sqrt(N) or some other numbers.

examples/qm/NDimTableTest_correct_output.txt

@@ -47,7 +47,10 @@ constructed, print now ---- checking TENSOR PRODUCT
 221000 (221000)   442000 (442000)   663000 (663000)   884000 (884000)   1.105e+06 (1105000)   
 231000 (231000)   462000 (462000)   693000 (693000)   924000 (924000)   1.155e+06 (1155000)   
 
-checking print() call
+checking minReal.T3: 1000  1000
+checking maxReal.T3: 1155000  1155000
+checking sumAll.T3 : 20880000  20880000
+checking print() call on T3
 
 1000 2000 3000 4000 5000 
 11000 22000 33000 44000 55000 
@@ -329,6 +332,22 @@ conversion! rank:2
 101 111 121 131 
 201 211 221 231 
 
+new tensor rank=1
+new tensor dim_n: 8, 
+new tensor total=8
+new tensor rank=1
+new tensor dim_n: 8, 
+new tensor total=8
+(10,0) (20,0) (1,0) (-4,0) (5,0) (6,0) (0,0) (1,0) 
+(0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) 
+move failed! rank:1
+(10,0) (20,0) (1,0) (-4,0) (5,0) (6,0) (0,0) (1,0) 
+(39,0) (18.435,-7.36396) (14,-29) (-8.43503,-5.36396) (-7,0) (-8.43503,5.36396) (14,29) (18.435,7.36396) 
+move failed! rank:1
+(10,0) (20,0) (1,0) (-4,0) (5,0) (6,0) (0,0) (1,0) 
+(39,0) (18.435,-7.36396) (14,-29) (-8.43503,-5.36396) (-7,0) (-8.43503,5.36396) (14,29) (18.435,7.36396) 
+last printed tensor, out.maxReal() is 39
+last printed tensor, out.minReal() is -8.43503
 ----- fftwf_FREE
 ----- fftwf_FREE
 ----- fftwf_FREE

lib/base/NDimTable.hpp

@@ -53,6 +53,10 @@ class NDimTable : private std::vector<K
 		typedef std::size_t        size_type;
 		typedef std::ptrdiff_t     difference_type;
 
+		template<typename Num> struct real_type                     {typedef Num type;};
+		template<typename Num> struct real_type<std::complex<Num> > {typedef Num type;};
+		typedef typename real_type<K>::type not_complex;
+
 		template<typename L> friend class NDimTable;
 	public:
 		typedef std::vector<std::size_t>  DimN;
@@ -80,16 +84,18 @@ class NDimTable : private std::vector<K
 			for(std::size_t i=0 ; i<others.size() ; i++) { dim_n.insert( dim_n.end(), others[i]->dim_n.begin(), others[i]->dim_n.end() ); };
 			total  = 1;            // total numer of elements (really explodes)
 			BOOST_FOREACH( std::size_t d_n, dim_n ) { total *= d_n; };
-
-//			printDebugInfo();
+			#ifdef DEBUG_NDIMTABLE
+			printDebugInfo();
+			#endif
 		}
 		inline void calcDimRankTotal(const std::vector<std::size_t>& d) {
 			dim_n  = d;
 			rank_d = dim_n.size();
 			total  = 1;
-			BOOST_FOREACH( std::size_t d_n, dim_n ) { total *= d_n; };
-
-//			printDebugInfo();
+			BOOST_FOREACH( std::size_t d_n, dim_n ) { total *= d_n; };	
+			#ifdef DEBUG_NDIMTABLE
+			printDebugInfo();
+			#endif
 		}
 		// element access: http://www.fftw.org/fftw3_doc/Row_002dmajor-Format.html#Row_002dmajor-Format
 		// element is located at the position iᵈ⁻¹ + nᵈ⁻¹ *(iᵈ⁻² + nᵈ⁻² *( ... n³*( i² + n²*(i¹ + n¹ * i⁰))))
@@ -126,18 +132,24 @@ class NDimTable : private std::vector<K
 		NDimTable(const NDimTable& other) 
 			: parent(static_cast<const parent&>(other)), rank_d(other.rank_d), dim_n(other.dim_n), total(other.total) 
 		{
-//			std::cerr << "move failed! rank:" << rank_d << "\n";
+			#ifdef DEBUG_NDIMTABLE
+			std::cerr << "move failed! rank:" << rank_d << "\n";
+			#endif
 		};
 		template<typename L> NDimTable(const NDimTable<L>& other) 
 			: parent(other.begin(),other.end()), rank_d(other.rank_d), dim_n(other.dim_n), total(other.total) 
 		{
+			#ifdef DEBUG_NDIMTABLE
 			std::cerr << "conversion! rank:" << rank_d << "\n";
+			#endif
 		};
 		// move constructor
 		NDimTable(NDimTable&& other)
 			: parent(static_cast<parent&&>(other)), rank_d(std::move(other.rank_d)), dim_n(std::move(other.dim_n)), total(std::move(other.total))
 		{
-//			std::cerr << "moved! rank:" << rank_d << "\n";
+			#ifdef DEBUG_NDIMTABLE
+			std::cerr << "moved! rank:" << rank_d << "\n";
+			#endif
 			other={};
 		};
 		// tensor product constructor
@@ -203,8 +215,9 @@ class NDimTable : private std::vector<K
 		// FIXME: should be 'K'-type not 'double'-type. But min(), max() works only with real numbers. 
 		// FIXME: this is because potential should be real valued (but isn't)
 		// FIXME: if it's 'double' here then better it should be 'Real' so that changing precision works correctly
-		double min() const {double ret(std::real(this->front())); for(K v : (*this)){ret = std::min(std::real(v),ret);}; return ret;};
-		double max() const {double ret(std::real(this->front())); for(K v : (*this)){ret = std::max(std::real(v),ret);}; return ret;};
+		not_complex minReal() const {not_complex ret(std::real(this->front())); for(K v : (*this)){ret  = std::min(std::real(v),ret);}; return ret;};
+		not_complex maxReal() const {not_complex ret(std::real(this->front())); for(K v : (*this)){ret  = std::max(std::real(v),ret);}; return ret;};
+		K           sumAll()  const {K           ret(          0             ); for(K v : (*this)){ret += v                         ;}; return ret;};
 		// !!!!!!!!!!!
 		// !IMPORTANT! for effciency, these do not copy construct new data, they modify in-place!
 		NDimTable& abs()           {std::transform(this->begin(),this->end(),this->begin(),[ ](K& v){return std::abs(v    );}); return *this;};

pkg/qm/QMInteraction.cpp

@@ -8,7 +8,6 @@
 
 YADE_PLUGIN(
 	(Ip2_QMParameters_QMParameters_QMPotPhysics)
-	(Ip2_Material_QMParameters_QMPotPhysics)
 	(Law2_QMPotGeometry_QMPotPhysics_QMPotPhysics)
 	);
 
@@ -34,10 +33,6 @@ void Ip2_QMParameters_QMParameters_QMPotPhysics::go(
 	std::cerr <<"####### iphys created in QMInteraction\n";
 }
 
-CREATE_LOGGER(Ip2_Material_QMParameters_QMPotPhysics);
-void Ip2_Material_QMParameters_QMPotPhysics::go(const shared_ptr<Material>& pp1, const shared_ptr<Material>& pp2, const shared_ptr<Interaction>& interaction){
-}
-
 /*********************************************************************************
 *
 * L A W 2   In DEM it was used to calculate Fn and Fs between two interacting bodies

pkg/qm/QMInteraction.hpp

@@ -35,40 +35,14 @@ class Ip2_QMParameters_QMParameters_QMPotPhysics: public IPhysFunctor
 			, // base class
 			IPhysFunctor
 			, // class description
-			"Currently does nothing"
+			"Convert :yref:`QMParameters` instance and :yref:`QMParameters` instance to \
+			:yref:`QMPotPhysics` with corresponding parameters."
 			, // attributes, public variables
 //			((long,nothing,10,,"placeholder"))
 		);
 };
 REGISTER_SERIALIZABLE(Ip2_QMParameters_QMParameters_QMPotPhysics);
 
-
-/*! @brief When any Material meets the QMParameters a QMPotPhysics is made.
- *
- * This can be used for infinite potential wells, where the DEM boxes serve as walls.
- *
- */
-
-class Ip2_Material_QMParameters_QMPotPhysics: public IPhysFunctor
-{
-	public:
-		virtual void go(const shared_ptr<Material>& pp1, const shared_ptr<Material>& pp2, const shared_ptr<Interaction>& interaction);
-		FUNCTOR2D(Material,QMParameters);
-		DECLARE_LOGGER;
-		YADE_CLASS_BASE_DOC_ATTRS(
-			  // class name
-			Ip2_Material_QMParameters_QMPotPhysics
-			, // base class
-			IPhysFunctor
-			, // class description
-			"Convert :yref:`QMParameters` instance and :yref:`Material` instance to \
-			:yref:`QMPotPhysics` with corresponding parameters."
-			, // attributes, public variables
-		);
-};
-REGISTER_SERIALIZABLE(Ip2_Material_QMParameters_QMPotPhysics);
-
-
 /*********************************************************************************
 *
 * L A W 2   In DEM it was used to calculate Fn and Fs between two interacting bodies

pkg/qm/QMStateDiscrete.hpp

@@ -40,7 +40,7 @@ or directly by filling in the discrete values in the table. It is used for numer
 			((bool      ,firstRun,true,Attr::readonly,"It is used to mark that postLoad() already generated the wavefunction from its creator analytic function."))
 			((boost::shared_ptr<QMStateAnalytic>,creator,,Attr::hidden,"Analytic wavepacket used to create the discretized version for calculations. The analytic shape can be anything: square packets, triangle, Gaussian - as long as it is normalized."))
 			((vector<size_t>,gridSize,vector<size_t>({}),,"Lattice grid size used to describe the wave function. For FFT purposes that should be a power of 2."))
-			((vector<Real>,size,vector<Real>({}),,"Wavepacket size in position representation space."))
+			((vector<Real>,size,vector<Real>({}),,"Wavepacket size in position representation space, for each DOF. Can be higher than 4D due to tensor products between wavefunctions."))
 			, // constructor
 			createIndex();
 			, // python bindings

pkg/qm/SchrodingerPropagator.cpp

@@ -65,7 +65,7 @@ Real SchrodingerKosloffPropagator::eMin()
 		}
 	};
 	// FIXME end
-	return ((Vpsi.rank()!=0) ? (Vpsi.min()) : (0));
+	return ((Vpsi.rank()!=0) ? (Vpsi.minReal()) : (0));
 };
 
 Real SchrodingerKosloffPropagator::eMax()
@@ -92,7 +92,7 @@ Real SchrodingerKosloffPropagator::eMax()
 			else               Vpsi+=igeom->potentialValues;  // ψᵥ: V = ∑Vᵢ // FIXME i używając jej rozmiar bym tworzył potencjał?
 		}
 	};
-	ret += ((Vpsi.rank()!=0) ? (Vpsi.max()) : (0));
+	ret += ((Vpsi.rank()!=0) ? (Vpsi.maxReal()) : (0));
 	// FIXME end
 
 	return ret;