FChebM2LHandler.hpp 18.7 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 64 65 66 67 68 69 70 71
 *
 * @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)
	{
72
		const char precision_type = (typeid(FReal)==typeid(double) ? 'd' : 'f');
73
		std::stringstream stream;
74
		stream << "m2l_k"<< MatrixKernelClass::Identifier << "_" << precision_type
75 76 77 78 79 80 81
					 << "_o" << order << "_e" << epsilon << ".bin";
		return stream.str();
	}

	
public:
	FChebM2LHandler(const FReal _epsilon)
messner's avatar
messner committed
82
		: MatrixKernel(), U(NULL), C(NULL), B(NULL), epsilon(_epsilon), rank(0)
83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101
	{}

	~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);
102 103 104 105

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


106
		// write info
107
		std::cout << "Compressed and set M2L operators (" << long(sizeM2L) << " B) in "
108
							<< time.tacAndElapsed() << "sec."	<< std::endl;
109 110 111 112 113 114 115 116 117
	}

	/**
	 * 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();
118 119 120 121 122 123 124 125 126 127
	}

	/**
	 * 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
128
	static unsigned int ComputeAndCompress(const FReal epsilon, FReal* &U, FReal* &C, FReal* &B);
129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155

	/**
	 * 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
156
    FBlas::gemva(nnodes, rank, 1., U, const_cast<FReal*>(X), x);
157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175
  }

  /**
	 * 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
176
    FBlas::gemva(rank, rank, scale, C + idx*rank*rank, const_cast<FReal*>(Y), X);
177 178 179 180 181
  }
  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
182
    FBlas::gemva(rank, rank, scale, C + idx*rank*rank, const_cast<FReal*>(Y), X);
183 184 185 186 187
  }
  void applyC(FReal CellWidth,
							const FReal *const Y, FReal *const X) const
  {
		const FReal scale(MatrixKernel.getScaleFactor(CellWidth));
messner's avatar
messner committed
188
    FBlas::gemva(rank, rank * 343, scale, C, const_cast<FReal*>(Y), X);
189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220
  }

  /**
	 * 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
221
unsigned int
222 223 224 225 226
FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompress(const FReal epsilon,
																															FReal* &U,
																															FReal* &C,
																															FReal* &B)
{
227 228 229
	// allocate memory and store compressed M2L operators
	if (U||C||B) throw std::runtime_error("Compressed M2L operators are already set");

230
	// interpolation points of source (Y) and target (X) cell
231
	FPoint X[nnodes], Y[nnodes];
232
	// set roots of target cell (X)
233
	FChebTensor<order>::setRoots(FPoint(0.,0.,0.), FReal(2.), X);
234 235 236 237 238 239 240 241 242 243 244 245 246
	// 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)
247
					const FPoint cy(FReal(2.*i), FReal(2.*j), FReal(2.*k));
248 249 250 251 252 253 254 255 256 257 258 259 260 261
					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");
262 263 264 265 266 267 268 269 270 271 272 273


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

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

278

279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296
	// 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++;
297
				} else FBlas::setzero(rank*rank, C + idx*rank*rank);
298 299 300 301 302
			}
	if (counter != ninteractions)
		throw std::runtime_error("Number of interactions must correspond to 316");
	delete [] _C;
		
303 304 305 306 307 308 309 310 311

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


312 313 314 315 316 317 318 319 320 321 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
	// 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
356
	std::cout << "Compressed M2L operators ("<< rank << ") stored in binary file "	<< filename
357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372
						<< " 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();
373 374 375 376 377 378
	if (size<=0) {
		std::cout << "Info: The requested binary file " << filename
							<< " does not yet exist. Compute it now ... " << std::endl;
		this->ComputeAndCompressAndStoreInBinaryFileAndReadFromFileAndSet();
		return;
	} 
379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398
	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
399 400
	std::cout << "Compressed M2L operators (" << rank << ") read from binary file "
						<< filename << " in " << time.tacAndElapsed() << "sec."	<< std::endl;
401 402
}

messner's avatar
messner committed
403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423
/*
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;
}
*/
424 425 426 427 428 429


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

430
/**
messner's avatar
messner committed
431
 * Computes the low-rank \f$k\f$ based on \f$\|K-U\Sigma_kV^\top\|_F \le
432 433 434
 * \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
435
 * the low-rank follows as \f$\|K-U\Sigma_kV^\top\|_F^2 = \sum_{i=k+1}^N
436 437 438 439 440 441 442
 * \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$
 */ 
443 444 445 446
template <int ORDER>
unsigned int getRank(const FReal singular_values[], const double eps)
{
	enum {nnodes = TensorTraits<ORDER>::nnodes};
447 448 449 450 451 452 453 454 455

	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;
456
	}
457
	throw std::runtime_error("rank cannot be larger than nnodes");
messner's avatar
messner committed
458
	return 0;
459 460
}

Matthias Messner's avatar
Matthias Messner committed
461 462
unsigned int getRank(const FReal singular_values[], const unsigned int size, const double eps)
{
463
	const FReal nrm2 = FBlas::scpr(size, singular_values, singular_values);
Matthias Messner's avatar
Matthias Messner committed
464 465 466 467 468 469 470 471 472
	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;
}

473 474

/**
messner's avatar
messner committed
475 476 477 478
 * 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.
479 480 481 482 483 484 485 486
 *
 * @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>
487 488
unsigned int Compress(const FReal epsilon, const unsigned int ninteractions,
											FReal* &U,	FReal* &C, FReal* &B)
489 490 491
{
	// compile time constants
	enum {order = ORDER,
492
				nnodes = TensorTraits<ORDER>::nnodes};
493 494 495 496 497 498 499 500 501 502 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

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