Une MAJ de sécurité est nécessaire sur notre version actuelle. Elle sera effectuée lundi 02/08 entre 12h30 et 13h. L'interruption de service devrait durer quelques minutes (probablement moins de 5 minutes).

FChebM2LHandler.hpp 19.3 KB
Newer Older
BRAMAS Berenger's avatar
BRAMAS Berenger committed
1
// ===================================================================================
2 3 4 5 6 7 8 9 10 11 12 13 14
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger 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".
BRAMAS Berenger's avatar
BRAMAS Berenger committed
15
// ===================================================================================
16 17 18 19 20 21 22 23
#ifndef FCHEBM2LHANDLER_HPP
#define FCHEBM2LHANDLER_HPP

#include <numeric>
#include <stdexcept>
#include <string>
#include <sstream>
#include <fstream>
24
#include <typeinfo>
25

26 27
#include "../../Utils/FBlas.hpp"
#include "../../Utils/FTic.hpp"
28 29 30 31 32

#include "./FChebTensor.hpp"


template <int ORDER>
33 34
unsigned int Compress(const FReal epsilon, const unsigned int ninteractions,
											FReal* &U,	FReal* &C, FReal* &B);
35 36 37 38 39 40 41


/**
 * @author Matthias Messner (matthias.messner@inria.fr)
 * @class FChebM2LHandler
 * Please read the license
 *
messner's avatar
messner committed
42 43 44 45 46 47 48 49 50 51 52
 * This class precomputes and compresses the M2L operators
 * \f$[K_1,\dots,K_{316}]\f$ for all (\f$7^3-3^3 = 316\f$ possible interacting
 * cells in the far-field) interactions for the Chebyshev interpolation
 * approach. The class uses the compression via a truncated SVD and represents
 * the compressed M2L operator as \f$K_t \sim U C_t B^\top\f$ with
 * \f$t=1,\dots,316\f$. The truncation rank is denoted by \f$r\f$ and is
 * determined by the prescribed accuracy \f$\varepsilon\f$. Hence, the
 * originally \f$K_t\f$ of size \f$\ell^3\times\ell^3\f$ times \f$316\f$ for
 * all interactions is reduced to only one \f$U\f$ and one \f$B\f$, each of
 * size \f$\ell^3\times r\f$, and \f$316\f$ \f$C_t\f$, each of size \f$r\times
 * r\f$.
53
 *
54 55 56 57 58 59 60 61 62 63
 * PB: FChebM2LHandler does not seem to support non_homogeneous kernels!
 * In fact nothing appears to handle this here (i.e. adapt scaling and storage 
 * to MatrixKernelClass::Type). Given the relatively important cost of the 
 * Chebyshev variant, it is probably a choice not to have implemented this 
 * feature here but instead in the ChebyshevSym variant. But what if the 
 * kernel is non homogeneous and non symmetric (e.g. Dislocations)... 
 * 
 * TODO Specialize class (see UnifM2LHandler) OR prevent from using this 
 * class with non homogeneous kernels ?!
 *
64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
 * @tparam ORDER interpolation order \f$\ell\f$
 */
template <int ORDER, class MatrixKernelClass>
class FChebM2LHandler : FNoCopyable
{
	enum {order = ORDER,
				nnodes = TensorTraits<ORDER>::nnodes,
				ninteractions = 316}; // 7^3 - 3^3 (max num cells in far-field)

	const MatrixKernelClass MatrixKernel;

	FReal *U, *C, *B;
	const FReal epsilon; //<! accuracy which determines trucation of SVD
	unsigned int rank;   //<! truncation rank, satisfies @p epsilon


	static const std::string getFileName(FReal epsilon)
	{
82
		const char precision_type = (typeid(FReal)==typeid(double) ? 'd' : 'f');
83
		std::stringstream stream;
84
		stream << "m2l_k"<< MatrixKernelClass::Identifier << "_" << precision_type
85 86 87 88 89 90 91
					 << "_o" << order << "_e" << epsilon << ".bin";
		return stream.str();
	}

	
public:
	FChebM2LHandler(const FReal _epsilon)
messner's avatar
messner committed
92
		: MatrixKernel(), U(NULL), C(NULL), B(NULL), epsilon(_epsilon), rank(0)
93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111
	{}

	~FChebM2LHandler()
	{
		if (U != NULL) delete [] U;
		if (B != NULL) delete [] B;
		if (C != NULL) delete [] C;
	}

	/**
	 * Computes, compresses and sets the matrices \f$Y, C_t, B\f$
	 */
	void ComputeAndCompressAndSet()
	{
		// measure time
		FTic time; time.tic();
		// check if aready set
		if (U||C||B) throw std::runtime_error("Compressed M2L operator already set");
		rank = ComputeAndCompress(epsilon, U, C, B);
112 113 114 115

    unsigned long sizeM2L = 343*rank*rank*sizeof(FReal);


116
		// write info
117
		std::cout << "Compressed and set M2L operators (" << long(sizeM2L) << " B) in "
118
							<< time.tacAndElapsed() << "sec."	<< std::endl;
119 120 121 122 123 124 125 126 127
	}

	/**
	 * Computes, compresses, writes to binary file, reads it and sets the matrices \f$Y, C_t, B\f$
	 */
	void ComputeAndCompressAndStoreInBinaryFileAndReadFromFileAndSet()
	{
		FChebM2LHandler<ORDER,MatrixKernelClass>::ComputeAndCompressAndStoreInBinaryFile(epsilon);
		this->ReadFromBinaryFileAndSet();
128 129 130 131 132 133 134 135 136 137
	}

	/**
	 * Computes and compressed all \f$K_t\f$.
	 *
	 * @param[in] epsilon accuracy
	 * @param[out] U matrix of size \f$\ell^3\times r\f$
	 * @param[out] C matrix of size \f$r\times 316 r\f$ storing \f$[C_1,\dots,C_{316}]\f$
	 * @param[out] B matrix of size \f$\ell^3\times r\f$
	 */
messner's avatar
messner committed
138
	static unsigned int ComputeAndCompress(const FReal epsilon, FReal* &U, FReal* &C, FReal* &B);
139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165

	/**
	 * Computes, compresses and stores the matrices \f$Y, C_t, B\f$ in a binary
	 * file
	 */
	static void ComputeAndCompressAndStoreInBinaryFile(const FReal epsilon);

	/**
	 * Reads the matrices \f$Y, C_t, B\f$ from the respective binary file
	 */
	void ReadFromBinaryFileAndSet();
		

	/**
	 * @return rank of the SVD compressed M2L operators
	 */
	unsigned int getRank() const {return rank;}

  /**
	 * Expands potentials \f$x+=UX\f$ of a target cell. This operation can be
	 * seen as part of the L2L operation.
	 *
	 * @param[in] X compressed local expansion of size \f$r\f$
	 * @param[out] x local expansion of size \f$\ell^3\f$
	 */
  void applyU(const FReal *const X, FReal *const x) const
  {
messner's avatar
messner committed
166
    FBlas::gemva(nnodes, rank, 1., U, const_cast<FReal*>(X), x);
167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185
  }

  /**
	 * Compressed M2L operation \f$X+=C_tY\f$, where \f$Y\f$ is the compressed
	 * multipole expansion and \f$X\f$ is the compressed local expansion, both
	 * of size \f$r\f$. The index \f$t\f$ denotes the transfer vector of the
	 * target cell to the source cell.
	 *
	 * @param[in] transfer transfer vector
	 * @param[in] Y compressed multipole expansion
	 * @param[out] X compressed local expansion
	 * @param[in] CellWidth needed for the scaling of the compressed M2L operators which are based on a homogeneous matrix kernel computed for the reference cell width \f$w=2\f$, ie in \f$[-1,1]^3\f$.
	 */
  void applyC(const int transfer[3], FReal CellWidth,
							const FReal *const Y, FReal *const X) const
  {
		const unsigned int idx
			= (transfer[0]+3)*7*7 + (transfer[1]+3)*7 + (transfer[2]+3);
		const FReal scale(MatrixKernel.getScaleFactor(CellWidth));
messner's avatar
messner committed
186
    FBlas::gemva(rank, rank, scale, C + idx*rank*rank, const_cast<FReal*>(Y), X);
187 188 189 190 191
  }
  void applyC(const unsigned int idx, FReal CellWidth,
							const FReal *const Y, FReal *const X) const
  {
		const FReal scale(MatrixKernel.getScaleFactor(CellWidth));
messner's avatar
messner committed
192
    FBlas::gemva(rank, rank, scale, C + idx*rank*rank, const_cast<FReal*>(Y), X);
193 194 195 196 197
  }
  void applyC(FReal CellWidth,
							const FReal *const Y, FReal *const X) const
  {
		const FReal scale(MatrixKernel.getScaleFactor(CellWidth));
messner's avatar
messner committed
198
    FBlas::gemva(rank, rank * 343, scale, C, const_cast<FReal*>(Y), X);
199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230
  }

  /**
	 * Compresses densities \f$Y=B^\top y\f$ of a source cell. This operation
	 * can be seen as part of the M2M operation.
	 *
	 * @param[in] y multipole expansion of size \f$\ell^3\f$
	 * @param[out] Y compressed multipole expansion of size \f$r\f$
	 */
  void applyB(FReal *const y, FReal *const Y) const
  {
    FBlas::gemtv(nnodes, rank, 1., B, y, Y);
  }


};






//////////////////////////////////////////////////////////////////////
// definition ////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////






template <int ORDER, class MatrixKernelClass>
messner's avatar
messner committed
231
unsigned int
232 233 234 235 236
FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompress(const FReal epsilon,
																															FReal* &U,
																															FReal* &C,
																															FReal* &B)
{
237 238 239
	// allocate memory and store compressed M2L operators
	if (U||C||B) throw std::runtime_error("Compressed M2L operators are already set");

240
	// interpolation points of source (Y) and target (X) cell
BRAMAS Berenger's avatar
BRAMAS Berenger committed
241
	FPoint X[nnodes], Y[nnodes];
242
	// set roots of target cell (X)
BRAMAS Berenger's avatar
BRAMAS Berenger committed
243
	FChebTensor<order>::setRoots(FPoint(0.,0.,0.), FReal(2.), X);
244 245 246 247 248 249 250 251 252 253 254 255 256
	// init matrix kernel
	const MatrixKernelClass MatrixKernel;

	// allocate memory and compute 316 m2l operators
	FReal *_U, *_C, *_B;
	_U = _B = NULL;
	_C = new FReal [nnodes*nnodes * ninteractions];
	unsigned int counter = 0;
	for (int i=-3; i<=3; ++i) {
		for (int j=-3; j<=3; ++j) {
			for (int k=-3; k<=3; ++k) {
				if (abs(i)>1 || abs(j)>1 || abs(k)>1) {
					// set roots of source cell (Y)
BRAMAS Berenger's avatar
BRAMAS Berenger committed
257
					const FPoint cy(FReal(2.*i), FReal(2.*j), FReal(2.*k));
258 259 260 261 262 263 264 265 266 267 268 269 270 271
					FChebTensor<order>::setRoots(cy, FReal(2.), Y);
					// evaluate m2l operator
					for (unsigned int n=0; n<nnodes; ++n)
						for (unsigned int m=0; m<nnodes; ++m)
							_C[counter*nnodes*nnodes + n*nnodes + m]
								= MatrixKernel.evaluate(X[m], Y[n]);
					// increment interaction counter
					counter++;
				}
			}
		}
	}
	if (counter != ninteractions)
		throw std::runtime_error("Number of interactions must correspond to 316");
272 273 274 275 276 277 278 279 280 281 282 283


	//////////////////////////////////////////////////////////		
	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
		}
	//////////////////////////////////////////////////////////		

284
	// svd compression of M2L
285
	const unsigned int rank	= Compress<ORDER>(epsilon, ninteractions, _U, _C, _B);
286 287
	if (!(rank>0)) throw std::runtime_error("Low rank must be larger then 0!");

288

289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306
	// store U
	U = new FReal [nnodes * rank];
	FBlas::copy(rank*nnodes, _U, U);
	delete [] _U;
	// store B
	B = new FReal [nnodes * rank];
	FBlas::copy(rank*nnodes, _B, B);
	delete [] _B;
	// store C
	counter = 0;
	C = new FReal [343 * rank*rank];
	for (int i=-3; i<=3; ++i)
		for (int j=-3; j<=3; ++j)
			for (int k=-3; k<=3; ++k) {
				const unsigned int idx = (i+3)*7*7 + (j+3)*7 + (k+3);
				if (abs(i)>1 || abs(j)>1 || abs(k)>1) {
					FBlas::copy(rank*rank, _C + counter*rank*rank, C + idx*rank*rank);
					counter++;
307
				} else FBlas::setzero(rank*rank, C + idx*rank*rank);
308 309 310 311 312
			}
	if (counter != ninteractions)
		throw std::runtime_error("Number of interactions must correspond to 316");
	delete [] _C;
		
313 314 315 316 317 318 319 320 321

	//////////////////////////////////////////////////////////		
	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
	}
	//////////////////////////////////////////////////////////		


322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365
	// return low rank
	return rank;
}






template <int ORDER, class MatrixKernelClass>
void
FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompressAndStoreInBinaryFile(const FReal epsilon)
{
	// measure time
	FTic time; time.tic();
	// start computing process
	FReal *U, *C, *B;
	U = C = B = NULL;
	const unsigned int rank = ComputeAndCompress(epsilon, U, C, B);
	// store into binary file
	const std::string filename(getFileName(epsilon));
	std::ofstream stream(filename.c_str(),
											 std::ios::out | std::ios::binary | std::ios::trunc);
	if (stream.good()) {
		stream.seekp(0);
		// 1) write number of interpolation points (int)
		int _nnodes = nnodes;
		stream.write(reinterpret_cast<char*>(&_nnodes), sizeof(int));
		// 2) write low rank (int)
		int _rank = rank;
		stream.write(reinterpret_cast<char*>(&_rank), sizeof(int));
		// 3) write U (rank*nnodes * FReal)
		stream.write(reinterpret_cast<char*>(U), sizeof(FReal)*rank*nnodes);
		// 4) write B (rank*nnodes * FReal)
		stream.write(reinterpret_cast<char*>(B), sizeof(FReal)*rank*nnodes);
		// 5) write 343 C (343 * rank*rank * FReal)
		stream.write(reinterpret_cast<char*>(C), sizeof(FReal)*rank*rank*343);
	} 	else throw std::runtime_error("File could not be opened to write");
	stream.close();
	// free memory
	if (U != NULL) delete [] U;
	if (B != NULL) delete [] B;
	if (C != NULL) delete [] C;
	// write info
366
	std::cout << "Compressed M2L operators ("<< rank << ") stored in binary file "	<< filename
367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382
						<< " in " << time.tacAndElapsed() << "sec."	<< std::endl;
}


template <int ORDER, class MatrixKernelClass>
void
FChebM2LHandler<ORDER, MatrixKernelClass>::ReadFromBinaryFileAndSet()
{
	// measure time
	FTic time; time.tic();
	// start reading process
	if (U||C||B) throw std::runtime_error("Compressed M2L operator already set");
	const std::string filename(getFileName(epsilon));
	std::ifstream stream(filename.c_str(),
											 std::ios::in | std::ios::binary | std::ios::ate);
	const std::ifstream::pos_type size = stream.tellg();
383 384 385 386 387 388
	if (size<=0) {
		std::cout << "Info: The requested binary file " << filename
							<< " does not yet exist. Compute it now ... " << std::endl;
		this->ComputeAndCompressAndStoreInBinaryFileAndReadFromFileAndSet();
		return;
	} 
389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408
	if (stream.good()) {
		stream.seekg(0);
		// 1) read number of interpolation points (int)
		int npts;
		stream.read(reinterpret_cast<char*>(&npts), sizeof(int));
		if (npts!=nnodes) throw std::runtime_error("nnodes and npts do not correspond");
		// 2) read low rank (int)
		stream.read(reinterpret_cast<char*>(&rank), sizeof(int));
		// 3) write U (rank*nnodes * FReal)
		U = new FReal [rank*nnodes];
		stream.read(reinterpret_cast<char*>(U), sizeof(FReal)*rank*nnodes);
		// 4) write B (rank*nnodes * FReal)
		B = new FReal [rank*nnodes];
		stream.read(reinterpret_cast<char*>(B), sizeof(FReal)*rank*nnodes);
		// 5) write 343 C (343 * rank*rank * FReal)
		C = new FReal [343 * rank*rank];
		stream.read(reinterpret_cast<char*>(C), sizeof(FReal)*rank*rank*343);
	}	else throw std::runtime_error("File could not be opened to read");
	stream.close();
	// write info
409 410
	std::cout << "Compressed M2L operators (" << rank << ") read from binary file "
						<< filename << " in " << time.tacAndElapsed() << "sec."	<< std::endl;
411 412
}

messner's avatar
messner committed
413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433
/*
unsigned int ReadRankFromBinaryFile(const std::string& filename)
{
	// start reading process
	std::ifstream stream(filename.c_str(), std::ios::in | std::ios::binary | std::ios::ate);
	const std::ifstream::pos_type size = stream.tellg();
	if (size<=0) throw std::runtime_error("The requested binary file does not exist.");
	unsigned int rank = -1;
	if (stream.good()) {
		stream.seekg(0);
		// 1) read number of interpolation points (int)
		int npts;
		stream.read(reinterpret_cast<char*>(&npts), sizeof(int));
		// 2) read low rank (int)
		stream.read(reinterpret_cast<char*>(&rank), sizeof(int));
		return rank;
	}	else throw std::runtime_error("File could not be opened to read");
	stream.close();
	return rank;
}
*/
434 435 436 437 438 439


//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////

440
/**
messner's avatar
messner committed
441
 * Computes the low-rank \f$k\f$ based on \f$\|K-U\Sigma_kV^\top\|_F \le
442 443 444
 * \epsilon \|K\|_F\f$, ie., the truncation rank of the singular value
 * decomposition. With the definition of the Frobenius norm \f$\|K\|_F =
 * \left(\sum_{i=1}^N \sigma_i^2\right)^{\frac{1}{2}}\f$ the determination of
messner's avatar
messner committed
445
 * the low-rank follows as \f$\|K-U\Sigma_kV^\top\|_F^2 = \sum_{i=k+1}^N
446 447 448 449 450 451 452
 * \sigma_i^2 \le \epsilon^2 \sum_{i=1}^N \sigma_i^2 = \epsilon^2
 * \|K\|_F^2\f$.
 *
 * @param[in] singular_values array of singular values ordered as \f$\sigma_1
 * \ge \sigma_2 \ge \dots \ge \sigma_N\f$
 * @param[in] eps accuracy \f$\epsilon\f$
 */ 
453 454 455 456
template <int ORDER>
unsigned int getRank(const FReal singular_values[], const double eps)
{
	enum {nnodes = TensorTraits<ORDER>::nnodes};
457 458 459 460 461 462 463 464 465

	FReal nrm2(0.);
	for (unsigned int k=0; k<nnodes; ++k)
		nrm2 += singular_values[k] * singular_values[k];

	FReal nrm2k(0.);
	for (unsigned int k=nnodes; k>0; --k) {
		nrm2k += singular_values[k-1] * singular_values[k-1];
		if (nrm2k > eps*eps * nrm2)	return k;
466
	}
467
	throw std::runtime_error("rank cannot be larger than nnodes");
messner's avatar
messner committed
468
	return 0;
469 470
}

Matthias Messner's avatar
Matthias Messner committed
471 472
unsigned int getRank(const FReal singular_values[], const unsigned int size, const double eps)
{
473
	const FReal nrm2 = FBlas::scpr(size, singular_values, singular_values);
Matthias Messner's avatar
Matthias Messner committed
474 475 476 477 478 479 480 481 482
	FReal nrm2k(0.);
	for (unsigned int k=size; k>0; --k) {
		nrm2k += singular_values[k-1] * singular_values[k-1];
		if (nrm2k > eps*eps * nrm2)	return k;
	}
	throw std::runtime_error("rank cannot be larger than size");
	return 0;
}

483 484

/**
messner's avatar
messner committed
485 486 487 488
 * Compresses \f$[K_1,\dots,K_{316}]\f$ in \f$C\f$. Attention: the matrices
 * \f$U,B\f$ are not initialized, no memory is allocated as input, as output
 * they store the respective matrices. The matrix \f$C\f$ stores
 * \f$[K_1,\dots,K_{316}]\f$ as input and \f$[C_1,\dots,C_{316}]\f$ as output.
489 490 491 492 493 494 495 496
 *
 * @param[in] epsilon accuracy
 * @param[out] U matrix of size \f$\ell^3\times r\f$
 * @param[in] C matrix of size \f$\ell^3\times 316 \ell^e\f$ storing \f$[K_1,\dots,K_{316}]\f$
 * @param[out] C matrix of size \f$r\times 316 r\f$ storing \f$[C_1,\dots,C_{316}]\f$
 * @param[out] B matrix of size \f$\ell^3\times r\f$
 */
template <int ORDER>
497 498
unsigned int Compress(const FReal epsilon, const unsigned int ninteractions,
											FReal* &U,	FReal* &C, FReal* &B)
499 500 501
{
	// compile time constants
	enum {order = ORDER,
502
				nnodes = TensorTraits<ORDER>::nnodes};
503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578

	// init SVD
	const unsigned int LWORK = 2 * (3*nnodes + ninteractions*nnodes);
	FReal *const WORK = new FReal [LWORK];
	
	// K_col ///////////////////////////////////////////////////////////
	FReal *const K_col = new FReal [ninteractions * nnodes*nnodes]; 
	for (unsigned int i=0; i<ninteractions; ++i)
		for (unsigned int j=0; j<nnodes; ++j)
			FBlas::copy(nnodes,
									C     + i*nnodes*nnodes + j*nnodes,
									K_col + j*ninteractions*nnodes + i*nnodes);
	// singular value decomposition
	FReal *const Q = new FReal [nnodes*nnodes];
	FReal *const S = new FReal [nnodes];
	const unsigned int info_col
		= FBlas::gesvd(ninteractions*nnodes, nnodes, K_col, S, Q, nnodes,
									 LWORK, WORK);
	if (info_col!=0)
		throw std::runtime_error("SVD did not converge with " + info_col);
	delete [] K_col;
	const unsigned int k_col = getRank<ORDER>(S, epsilon);

	// Q' -> B 
	B = new FReal [nnodes*k_col];
	for (unsigned int i=0; i<k_col; ++i)
		FBlas::copy(nnodes, Q+i, nnodes, B+i*nnodes, 1);

	// K_row //////////////////////////////////////////////////////////////
	FReal *const K_row = C;

	const unsigned int info_row
		= FBlas::gesvdSO(nnodes, ninteractions*nnodes, K_row, S, Q, nnodes,
										 LWORK, WORK);
	if (info_row!=0)
		throw std::runtime_error("SVD did not converge with " + info_row);
	const unsigned int k_row = getRank<ORDER>(S, epsilon);
	delete [] WORK;

	// Q -> U
	U = Q;

	// V' -> V
	FReal *const V = new FReal [nnodes*ninteractions * k_row];
	for (unsigned int i=0; i<k_row; ++i)
		FBlas::copy(nnodes*ninteractions, K_row+i, nnodes,
								V+i*nnodes*ninteractions, 1);

	// rank k(epsilon) /////////////////////////////////////////////////////
	const unsigned int k = k_row < k_col ? k_row : k_col;

	// C_row ///////////////////////////////////////////////////////////
	C = new FReal [ninteractions * k*k];
	for (unsigned int i=0; i<k; ++i) {
		FBlas::scal(nnodes*ninteractions, S[i], V + i*nnodes*ninteractions);
		for (unsigned int m=0; m<ninteractions; ++m)
			for (unsigned int j=0; j<k; ++j)
				C[m*k*k + j*k + i]
					= FBlas::scpr(nnodes,
												V + i*nnodes*ninteractions + m*nnodes,
												B + j*nnodes);
	}

	delete [] V;
	delete [] S;
	delete [] K_row;

	return k;
}




#endif // FCHEBM2LHANDLER_HPP

// [--END--]