From 655fc19203fb8da008d3a1017087f22b01124bbd Mon Sep 17 00:00:00 2001
From: messner <messner@2616d619-271b-44dc-8df4-d4a8f33a7222>
Date: Mon, 27 Feb 2012 08:14:39 +0000
Subject: [PATCH] update of weighting functionality of M2L operators (WIP)

git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/scalfmm/scalfmm/trunk@445 2616d619-271b-44dc-8df4-d4a8f33a7222
---
 Src/Chebyshev/FChebKernels.hpp    |  6 ++---
 Src/Chebyshev/FChebM2LHandler.hpp | 42 +++++++++++++++++++++++--------
 Src/Chebyshev/FChebTensor.hpp     | 10 ++++++--
 3 files changed, 41 insertions(+), 17 deletions(-)

diff --git a/Src/Chebyshev/FChebKernels.hpp b/Src/Chebyshev/FChebKernels.hpp
index 13641ae7..aa65a221 100644
--- a/Src/Chebyshev/FChebKernels.hpp
+++ b/Src/Chebyshev/FChebKernels.hpp
@@ -85,12 +85,10 @@ public:
 			Epsilon(inEpsilon)
 	{
 		// read precomputed compressed m2l operators from binary file
-		M2LHandler->ReadFromBinaryFileAndSet();
+		//M2LHandler->ReadFromBinaryFileAndSet();
+		M2LHandler->ComputeAndCompressAndSet();
 	}
 
-//	/** Default destructor */
-//	~FChebKernels()	{}
-
 
 	void P2M(CellClass* const LeafCell,
 					 const ContainerClass* const SourceParticles)
diff --git a/Src/Chebyshev/FChebM2LHandler.hpp b/Src/Chebyshev/FChebM2LHandler.hpp
index abfe6a22..7f17a414 100644
--- a/Src/Chebyshev/FChebM2LHandler.hpp
+++ b/Src/Chebyshev/FChebM2LHandler.hpp
@@ -86,7 +86,8 @@ public:
 		if (U||C||B) throw std::runtime_error("Compressed M2L operator already set");
 		rank = ComputeAndCompress(epsilon, U, C, B);
 		// write info
-		std::cout << "Compressed and set M2L operators in "	<< time.tacAndElapsed() << "sec."	<< std::endl;
+		std::cout << "Compressed and set M2L operators (" << rank << ") in "
+							<< time.tacAndElapsed() << "sec."	<< std::endl;
 	}
 
 	/**
@@ -205,6 +206,9 @@ FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompress(const FReal epsilo
 																															FReal* &C,
 																															FReal* &B)
 {
+	// allocate memory and store compressed M2L operators
+	if (U||C||B) throw std::runtime_error("Compressed M2L operators are already set");
+
 	// interpolation points of source (Y) and target (X) cell
 	F3DPosition X[nnodes], Y[nnodes];
 	// set roots of target cell (X)
@@ -237,13 +241,23 @@ FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompress(const FReal epsilo
 	}
 	if (counter != ninteractions)
 		throw std::runtime_error("Number of interactions must correspond to 316");
-		
+
+
+	//////////////////////////////////////////////////////////		
+	FReal weights[nnodes];
+	FChebTensor<order>::setRootOfWeights(weights);
+	for (unsigned int i=0; i<316; ++i)
+		for (unsigned int n=0; n<nnodes; ++n) {
+			FBlas::scal(nnodes, weights[n], _C+i*nnodes*nnodes + n,  nnodes); // scale rows
+			FBlas::scal(nnodes, weights[n], _C+i*nnodes*nnodes + n * nnodes); // scale cols
+		}
+	//////////////////////////////////////////////////////////		
+
 	// svd compression of M2L
 	const unsigned int rank	= Compress<ORDER>(epsilon, ninteractions, _U, _C, _B);
 	if (!(rank>0)) throw std::runtime_error("Low rank must be larger then 0!");
 
-	// allocate memory and store compressed M2L operators
-	if (U||C||B) throw std::runtime_error("Compressed M2L operator already set");
+
 	// store U
 	U = new FReal [nnodes * rank];
 	FBlas::copy(rank*nnodes, _U, U);
@@ -262,12 +276,21 @@ FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompress(const FReal epsilo
 				if (abs(i)>1 || abs(j)>1 || abs(k)>1) {
 					FBlas::copy(rank*rank, _C + counter*rank*rank, C + idx*rank*rank);
 					counter++;
-				} else FBlas::scal(rank*rank, FReal(0.), C + idx*rank*rank);
+				} else FBlas::setzero(rank*rank, C + idx*rank*rank);
 			}
 	if (counter != ninteractions)
 		throw std::runtime_error("Number of interactions must correspond to 316");
 	delete [] _C;
 		
+
+	//////////////////////////////////////////////////////////		
+	for (unsigned int n=0; n<nnodes; ++n) {
+		FBlas::scal(rank, FReal(1.) / weights[n], U+n, nnodes); // scale rows
+		FBlas::scal(rank, FReal(1.) / weights[n], B+n, nnodes); // scale rows
+	}
+	//////////////////////////////////////////////////////////		
+
+
 	// return low rank
 	return rank;
 }
@@ -312,8 +335,7 @@ FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompressAndStoreInBinaryFil
 	if (B != NULL) delete [] B;
 	if (C != NULL) delete [] C;
 	// write info
-	std::cout << "Compressed M2L operators (rank="<< rank
-						<< ") stored in binary file "	<< filename
+	std::cout << "Compressed M2L operators ("<< rank << ") stored in binary file "	<< filename
 						<< " in " << time.tacAndElapsed() << "sec."	<< std::endl;
 }
 
@@ -356,9 +378,8 @@ FChebM2LHandler<ORDER, MatrixKernelClass>::ReadFromBinaryFileAndSet()
 	}	else throw std::runtime_error("File could not be opened to read");
 	stream.close();
 	// write info
-	std::cout << "Compressed M2L operators read from binary file "
-						<< filename << " in " << time.tacAndElapsed() << "sec."
-						<< std::endl;
+	std::cout << "Compressed M2L operators (" << rank << ") read from binary file "
+						<< filename << " in " << time.tacAndElapsed() << "sec."	<< std::endl;
 }
 
 /*
@@ -532,7 +553,6 @@ unsigned int Compress(const FReal epsilon, const unsigned int ninteractions,
 
 
 
-
 #endif // FCHEBM2LHANDLER_HPP
 
 // [--END--]
diff --git a/Src/Chebyshev/FChebTensor.hpp b/Src/Chebyshev/FChebTensor.hpp
index 89fb7ed8..274ff2ab 100644
--- a/Src/Chebyshev/FChebTensor.hpp
+++ b/Src/Chebyshev/FChebTensor.hpp
@@ -76,8 +76,14 @@ public:
 		// weights in 3d (tensor structure)
 		unsigned int node_ids[nnodes][3];
 		setNodeIds(node_ids);
-		for (unsigned int n=0; n<nnodes; ++n) 
-			weights[n] = sqrt(weights_1d[node_ids[n][0]]*weights_1d[node_ids[n][1]]*weights_1d[node_ids[n][2]]);
+		for (unsigned int n=0; n<nnodes; ++n) {
+			FReal weight = FReal(1.); // no weighting
+			weight *= sqrt(weights_1d[node_ids[n][0]]*weights_1d[node_ids[n][1]]*weights_1d[node_ids[n][2]]); // 1/2
+			//weight *= weight; // 1
+			//weight *= weight; // 2
+			//weight *= weight; // 4
+			weights[n] = weight;
+		}
 	}
 
 
-- 
GitLab