Commit fafa16dc authored by COULAUD Olivier's avatar COULAUD Olivier

Fix some typos in the API

Add SVD in Blas
parent 077ce97b
// ===================================================================================
// Logiciel initial: ScalFmm Version 0.5
// Co-auteurs : Olivier Coulaud, Bérenger Bramas.
// Propriétaires : INRIA.
// Copyright © 2011-2012, diffusé sous les termes et conditions d’une licence propriétaire.
// Initial software: ScalFmm Version 0.5
// Co-authors: Olivier Coulaud, Bérenger Bramas.
// Owners: INRIA.
// Copyright © 2011-2012, spread under the terms and conditions of a proprietary license.
// Copyright ScalFmm 2014 I
// 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".
// ===================================================================================
//
#ifndef CKERNELAPI_H
#define CKERNELAPI_H
......@@ -31,7 +36,7 @@ typedef void* Scalfmm_Handle;
//< it gives the level of the cell, its morton index, it position in term of box at that level
//< and the spatial position of its center
typedef void* (*Callback_init_cell)(int level, long long morton_index, int* tree_position, double* spatial_position);
//< Function to destroy what have bee initalized by the user (should be give in Scalfmm_dealloc_handle)
//< Function to destroy what have bee initialized by the user (should be give in Scalfmm_dealloc_handle)
typedef void (*Callback_free_cell)(void*);
//< This function init an handle (and an tree based on the given properties)
......@@ -41,7 +46,7 @@ void Scalfmm_dealloc_handle(Scalfmm_Handle handle, Callback_free_cell cellDestro
//< This function should be used to insert an array of particle in the tree
//< The indexes are the one used on the particles operator
//< The posision of the particles should be composed of one triple per particle:
//< The position of the particles should be composed of one triple per particle:
//< xyzxyzxyz...
void Scalfmm_insert_array_of_particles(Scalfmm_Handle handle, int nbParticles, int* particleIndexes, double* particleXYZ);
//< To insert one particle only
......
// ===================================================================================
// Logiciel initial: ScalFmm Version 0.5
// Co-auteurs : Olivier Coulaud, Bérenger Bramas.
// Propriétaires : INRIA.
// Copyright © 2011-2012, diffusé sous les termes et conditions d’une licence propriétaire.
// Initial software: ScalFmm Version 0.5
// Co-authors: Olivier Coulaud, Bérenger Bramas.
// Owners: INRIA.
// Copyright © 2011-2012, spread under the terms and conditions of a proprietary license.
// Copyright ScalFmm 2014 I
// 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 "../../Src/Containers/FOctree.hpp"
#include "../../Src/Containers/FVector.hpp"
......@@ -156,10 +161,10 @@ typedef FFmmAlgorithmThread<OctreeClass, CoreCell, CoreContainerClass, CoreKerne
// Our scalfmm handle
struct ScalFmmCoreHandle {
struct ScalFmmCoreConfig {
//paramètres en lecture/écriture :
int treeHeight; // hombre de niveaux de l'arbre (int)
FReal boxWidth; // taille de la boîte racine (FReal)
FPoint boxCenter; // position du centre de la boîte racine (FReal[3])
// Read/Write parameter
int treeHeight; // Number of level in the octree
FReal boxWidth; // Simulation box size (root level)
FPoint boxCenter; // Center position of the box simulation(FReal[3])
};
ScalFmmCoreConfig config;
......
......@@ -20,7 +20,7 @@
#include "Utils/FMemUtils.hpp"
#include "Utils/FBlas.hpp"
#include "Utils/FComplexe.hpp"
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* This class is a spherical harmonic kernels using blas
......@@ -98,6 +98,8 @@ protected:
matrix[idxCol * FF_MATRIX_ROW_DIM + idxRow] = workMatrix[idxCol + idxRow * FF_MATRIX_COLUMN_DIM];
}
}
#ifdef BLAS_M2L_P
//
// Single Layer
//
......@@ -124,6 +126,7 @@ protected:
++idxRow ;
}
} //
#endif
preM2LTransitions[idxLevel][indexM2LTransition(idxX,idxY,idxZ)] = matrix;
}
}
......@@ -243,20 +246,41 @@ public:
*/
void multipoleToLocal(FComplexe*const FRestrict local_exp, const FComplexe* const FRestrict multipole_exp_src,
const FComplexe* const FRestrict M2L_Outer_transfer){
//
// Copy original vector and compute exp2nexp
FMemUtils::copyall<FComplexe>(temporaryMultiSource, multipole_exp_src, CellClass::GetPoleSize());
// Get a computable vector
preExpNExp(temporaryMultiSource);
#ifdef BLAS_SPHERICAL_COMPRESS
FComplexe * res = new FComplexe [FF_MATRIX_ROW_DIM] ;
this->CompressM2LMatrix(M2L_Outer_transfer,res,local_exp, temporaryMultiSource) ;
#endif
const FReal one[2] = {1.0 , 0.0};
// Matrix-vector product
FBlas::c_gemva(
FF_MATRIX_ROW_DIM,
FF_MATRIX_COLUMN_DIM,
one,
FComplexe::ToFReal(M2L_Outer_transfer),
FComplexe::ToFReal(temporaryMultiSource),
FComplexe::ToFReal(local_exp));
FComplexe::ToFReal(local_exp) );
#ifdef BLAS_SPHERICAL_COMPRESS
// Compare the two computations
std::cout << "\n ====== Compare ====== \n"<<std::endl;
FComplexe tmp ;
for(int idxCol =0 ;idxCol< FF_MATRIX_ROW_DIM ; ++idxCol ){ // Col
tmp.setRealImag( local_exp[idxCol].getReal() - res[idxCol].getReal(), local_exp[idxCol].getImag() - res[idxCol].getImag()) ;
std::cout << idxCol << " " << local_exp[idxCol] <<" " << res[idxCol] << "Error: " << tmp.norm()<< std::endl;
}
std::cout << std::endl ;
std::cout << "============================ \n"<<std::endl;
delete [] res;
#endif
#ifdef DEBUG_SPHERICAL_M2L
......@@ -266,6 +290,7 @@ public:
std::cout << temporaryMultiSource[idxCol] <<" ";
}
std::cout << std::endl ;
#ifdef PTINT_MATRIX
std::cout << "\n ====== MultipolToLocal MatrixTransfer ====== \n"<<std::endl;
std::cout << " FF_MATRIX_ROW_DIM: " << FF_MATRIX_ROW_DIM<<std::endl;
......@@ -278,7 +303,7 @@ public:
std::cout << std::endl ;
}
std::cout << std::endl ;
#endif
std::cout << "\n ====== Local expansion ====== \n"<<std::endl;
for(int idxCol =0 ;idxCol< FF_MATRIX_ROW_DIM ; ++idxCol ){ // Col
......@@ -290,6 +315,110 @@ public:
#endif
}
/** M2L
* Compress by SVD all the M2L transfer matrix
*/
void CompressM2LMatrix(const FComplexe* const FRestrict M2L_Matrix,FComplexe*const res,
FComplexe*const FRestrict local_exp, const FComplexe* const FRestrict multipole_exp_src)
{
int LWORK = 5*FF_MATRIX_ROW_DIM;
FComplexe *const work = new FComplexe [LWORK];
FReal *const rwork = new FReal [ LWORK];
// singular value decomposition
FComplexe *const matrix = new FComplexe [FF_MATRIX_ROW_DIM*FF_MATRIX_COLUMN_DIM];
FComplexe *const Q = new FComplexe [FF_MATRIX_COLUMN_DIM*FF_MATRIX_COLUMN_DIM];
FReal *const S = new FReal [ FF_MATRIX_ROW_DIM];
//
// Work on a copy of the M2L matrix
FMemUtils::copyall<FComplexe>(matrix, M2L_Matrix, FF_MATRIX_COLUMN_DIM*FF_MATRIX_COLUMN_DIM);
//
#ifdef PRINT_MATRIX
std::cout << "\n ====== matrix copy of MatrixTransfer ====== \n"<<std::endl;
std::cout << " FF_MATRIX_ROW_DIM: " << FF_MATRIX_ROW_DIM<<std::endl;
std::cout << " FF_MATRIX_COLUMN_DIM: " << FF_MATRIX_COLUMN_DIM<<std::endl;
for(int idxRow =0 ;idxRow< FF_MATRIX_ROW_DIM ; ++idxRow ){ // Row
for(int idxCol =0 ;idxCol< FF_MATRIX_COLUMN_DIM ; ++idxCol ){ // Col
std::cout << idxRow<< " " << idxCol << " " << matrix[idxCol * FF_MATRIX_ROW_DIM + idxRow] << std::endl ;
}
}
std::cout << std::endl ;
#endif
// Perform SVD
FBlas::c_gesvd(
FF_MATRIX_ROW_DIM,
FF_MATRIX_COLUMN_DIM,
FComplexe::ToFReal(matrix),
S,
FComplexe::ToFReal(Q),
FF_MATRIX_COLUMN_DIM,
LWORK,
FComplexe::ToFReal(work),
rwork );
//
std::cout << "\n ====== V**H Matrix ====== \n"<<std::endl;
std::cout << " FF_MATRIX_ROW_DIM: " << FF_MATRIX_ROW_DIM<<std::endl;
std::cout << " FF_MATRIX_COLUMN_DIM: " << FF_MATRIX_COLUMN_DIM<<std::endl;
for(int idxRow =0 ;idxRow< FF_MATRIX_ROW_DIM ; ++idxRow ){ // Row
std::cout << "Row="<<idxRow<<" : " ;
for(int idxCol =0 ;idxCol< FF_MATRIX_COLUMN_DIM ; ++idxCol ){ // Col
std::cout << matrix[idxCol * FF_MATRIX_ROW_DIM + idxRow] <<" ";
}
std::cout << std::endl ;
}
std::cout << std::endl ;
int rank = FF_MATRIX_ROW_DIM ;
/*
double eps= 1.0e-6 ;
for( int idxRow =0 ; idxRow< FF_MATRIX_ROW_DIM ; ++idxRow ){
if(S[idxRow]/ S[0] < eps) {
S[idxRow] = 0.0 ;
rank = rank - 1 ;
}
std::cout << " idxRow=" << idxRow << " "<< S[idxRow]/ S[0] <<std::endl;
}
*/
std::cout << " New rank " << rank << std::endl ;
//
// store matrices U,V and sigma.
//
//
// Compute the product
// call zgemv ('N',rank, N, ALPHA, VT, N, X, 1, BETA, YC, 1)
// Matrix-vector product
const FReal one[2] = {1.0 , 0.0};
FComplexe * tempory = new FComplexe[FF_MATRIX_ROW_DIM];
FBlas::c_gemva(
FF_MATRIX_ROW_DIM,
FF_MATRIX_COLUMN_DIM,
one,
FComplexe::ToFReal(matrix),
FComplexe::ToFReal(temporaryMultiSource),
FComplexe::ToFReal(tempory) );
// Apply S
for (int i = 0; i < rank ; ++i){
std::cout << i << " "<< S[i] << " "<< tempory[i] << std::endl;
res[i] *= 0.0;
tempory[i] *= S[i];
}
// Apply the U matrix
FBlas::c_gemva(
FF_MATRIX_ROW_DIM,
rank,
one,
FComplexe::ToFReal(Q),
FComplexe::ToFReal(tempory),
FComplexe::ToFReal(res) );
//
delete [] matrix ;
delete [] work ;
delete [] rwork ;
std::cout << " \n End Compress M2L" <<std::endl;
}
};
#endif // FSPHERICALBLASKERNEL_HPP
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment