Commit 5aafb270 authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

add methods to call first level API of Chameleon( with data conversions)

parent b9ba4f12
......@@ -8,8 +8,13 @@ extern "C" {
* Implement Init function by creating an instance of the Engine
* class. The ptr to the instance is the Handle.
*/
extern "C" IbGMResDr_handle Init_ibgmresdr(ibgmresdr_Arithmetic ARI,
void * Matrix, int dim, void * userEnv){
extern "C"
IbGMResDr_handle Init_ibgmresdr(ibgmresdr_Arithmetic ARI,
void* Matrix,
int dim,
void* userEnv)
{
switch(ARI){
case FLOAT:{
IbGMResDrEngine<float>* Engine =
......@@ -37,7 +42,8 @@ extern "C" IbGMResDr_handle Init_ibgmresdr(ibgmresdr_Arithmetic ARI,
}
default:
std::cout<<"Arithmetic not recognised, please see the documentation"
<< "in order to choose something we propose\nExiting\n";
<< "in order to choose something we propose\n"
<< "Exiting\n";
exit(0);
break;
}
......
This diff is collapsed.
......@@ -20,7 +20,8 @@ typedef struct matrix{
//Define callbacks
//Matrix block of vector products
void MatxBlockVect(void * mat, int nbVect, void * input, void **output, void * env){
void MatxBlockVect(void * mat, int nbVect, void * input, void **output, void * env)
{
Matrix * inMat = (Matrix *) mat;
Env * inEnv = (Env *) env;
double * toWrite = (double*) (*output);
......@@ -37,7 +38,8 @@ void MatxBlockVect(void * mat, int nbVect, void * input, void **output, void * e
}
//Dot product callback
void dot_product(int size, void * vectA, void * vectB , void * res, void * env){
void dot_product(int size, void * vectA, void * vectB , void * res, void * env)
{
Env * inEnv = (Env *) env; //Not used
//cast arrays ::
double * v1 = (double*) vectA;
......@@ -51,7 +53,8 @@ void dot_product(int size, void * vectA, void * vectB , void * res, void * env){
}
//Display Matrix
void display(Matrix * mat){
void display(Matrix * mat)
{
int dim = mat->dim;
int i,j;
for(i=0 ; i<dim ; ++i){
......@@ -62,20 +65,21 @@ void display(Matrix * mat){
}
}
int main(int ac, char ** av){
int main(int argc, char *argv[])
{
MORSE_Init(2, 0);
MORSE_Enable(MORSE_ERRORS);
int dim = 25;
if(ac > 1){
dim = atoi(av[1]);
}
if (argc > 1)
dim = atoi(argv[1]);
int seed = time(NULL);
srand(seed);
/* int seed = time(NULL); */
/* srand(seed); */
//Create Matrix
Matrix * my_matrix = malloc(sizeof(Matrix));
Matrix *my_matrix = malloc(sizeof(Matrix));
my_matrix->datas = malloc(sizeof(double)*dim*dim);
my_matrix->dim = dim;
......
......@@ -23,6 +23,11 @@ pkg_check_modules(CHAMELEON chameleon REQUIRED)
pkg_check_modules(STARPU starpu-1.1 REQUIRED)
pkg_check_modules(STARPU_MPI starpumpi-1.1 REQUIRED)
if(CHAMELEON_FOUND)
add_definitions( -DMORSE_COMPLEX_CPP )
endif()
if(CBLAS_FOUND)
message(STATUS "CBLAS_FOUND = ${CBLAS_FOUND}")
message(STATUS "CBLAS_LIBRARIES = ${CBLAS_LIBRARIES}")
......@@ -58,9 +63,9 @@ option(IBBGMRESDR_BUILD_WARNINGS "Build with warning options" ON)
option(IBBGMRESDR_BUILD_DOC "Set to ON to build the Doxygen documentation " OFF )
if(IBBGMRESDR_DEBUG_MODE)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g3 -O0")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g3 -O0 -fmax-errors=2")
else()
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -fmax-errors=2")
endif()
if(IBBGMRESDR_BUILD_WARNINGS)
......
......@@ -88,39 +88,39 @@ CBLAS_INDEX cblas_izamax(const int N, const void *X, const int incX);
* ===========================================================================
*/
/*
/*
* Routines with standard 4 prefixes (s, d, c, z)
*/
void cblas_sswap(const int N, float *X, const int incX,
void cblas_sswap(const int N, float *X, const int incX,
float *Y, const int incY);
void cblas_scopy(const int N, const float *X, const int incX,
void cblas_scopy(const int N, const float *X, const int incX,
float *Y, const int incY);
void cblas_saxpy(const int N, const float alpha, const float *X,
const int incX, float *Y, const int incY);
void cblas_dswap(const int N, double *X, const int incX,
void cblas_dswap(const int N, double *X, const int incX,
double *Y, const int incY);
void cblas_dcopy(const int N, const double *X, const int incX,
void cblas_dcopy(const int N, const double *X, const int incX,
double *Y, const int incY);
void cblas_daxpy(const int N, const double alpha, const double *X,
const int incX, double *Y, const int incY);
void cblas_cswap(const int N, void *X, const int incX,
void cblas_cswap(const int N, void *X, const int incX,
void *Y, const int incY);
void cblas_ccopy(const int N, const void *X, const int incX,
void cblas_ccopy(const int N, const void *X, const int incX,
void *Y, const int incY);
void cblas_caxpy(const int N, const void *alpha, const void *X,
const int incX, void *Y, const int incY);
void cblas_zswap(const int N, void *X, const int incX,
void cblas_zswap(const int N, void *X, const int incX,
void *Y, const int incY);
void cblas_zcopy(const int N, const void *X, const int incX,
void cblas_zcopy(const int N, const void *X, const int incX,
void *Y, const int incY);
void cblas_zaxpy(const int N, const void *alpha, const void *X,
const int incX, void *Y, const int incY);
/*
/*
* Routines with S and D prefix only
*/
void cblas_srotg(float *a, float *b, float *c, float *s);
......@@ -138,7 +138,7 @@ void cblas_drotm(const int N, double *X, const int incX,
double *Y, const int incY, const double *P);
/*
/*
* Routines with S D C Z CS and ZD prefixes
*/
void cblas_sscal(const int N, const float alpha, float *X, const int incX);
......@@ -154,7 +154,7 @@ void cblas_zdscal(const int N, const double alpha, void *X, const int incX);
* ===========================================================================
*/
/*
/*
* Routines with standard 4 prefixes (S, D, C, Z)
*/
void cblas_sgemv(const CBLAS_LAYOUT layout,
......@@ -169,11 +169,11 @@ void cblas_sgbmv(CBLAS_LAYOUT layout,
const int incX, const float beta, float *Y, const int incY);
void cblas_strmv(CBLAS_LAYOUT layout, CBLAS_UPLO Uplo,
CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag,
const int N, const float *A, const int lda,
const int N, const float *A, const int lda,
float *X, const int incX);
void cblas_stbmv(CBLAS_LAYOUT layout, CBLAS_UPLO Uplo,
CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag,
const int N, const int K, const float *A, const int lda,
const int N, const int K, const float *A, const int lda,
float *X, const int incX);
void cblas_stpmv(CBLAS_LAYOUT layout, CBLAS_UPLO Uplo,
CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag,
......@@ -202,11 +202,11 @@ void cblas_dgbmv(CBLAS_LAYOUT layout,
const int incX, const double beta, double *Y, const int incY);
void cblas_dtrmv(CBLAS_LAYOUT layout, CBLAS_UPLO Uplo,
CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag,
const int N, const double *A, const int lda,
const int N, const double *A, const int lda,
double *X, const int incX);
void cblas_dtbmv(CBLAS_LAYOUT layout, CBLAS_UPLO Uplo,
CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag,
const int N, const int K, const double *A, const int lda,
const int N, const int K, const double *A, const int lda,
double *X, const int incX);
void cblas_dtpmv(CBLAS_LAYOUT layout, CBLAS_UPLO Uplo,
CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag,
......@@ -235,11 +235,11 @@ void cblas_cgbmv(CBLAS_LAYOUT layout,
const int incX, const void *beta, void *Y, const int incY);
void cblas_ctrmv(CBLAS_LAYOUT layout, CBLAS_UPLO Uplo,
CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag,
const int N, const void *A, const int lda,
const int N, const void *A, const int lda,
void *X, const int incX);
void cblas_ctbmv(CBLAS_LAYOUT layout, CBLAS_UPLO Uplo,
CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag,
const int N, const int K, const void *A, const int lda,
const int N, const int K, const void *A, const int lda,
void *X, const int incX);
void cblas_ctpmv(CBLAS_LAYOUT layout, CBLAS_UPLO Uplo,
CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag,
......@@ -268,11 +268,11 @@ void cblas_zgbmv(CBLAS_LAYOUT layout,
const int incX, const void *beta, void *Y, const int incY);
void cblas_ztrmv(CBLAS_LAYOUT layout, CBLAS_UPLO Uplo,
CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag,
const int N, const void *A, const int lda,
const int N, const void *A, const int lda,
void *X, const int incX);
void cblas_ztbmv(CBLAS_LAYOUT layout, CBLAS_UPLO Uplo,
CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag,
const int N, const int K, const void *A, const int lda,
const int N, const int K, const void *A, const int lda,
void *X, const int incX);
void cblas_ztpmv(CBLAS_LAYOUT layout, CBLAS_UPLO Uplo,
CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag,
......@@ -290,7 +290,7 @@ void cblas_ztpsv(CBLAS_LAYOUT layout, CBLAS_UPLO Uplo,
const int N, const void *Ap, void *X, const int incX);
/*
/*
* Routines with S and D prefixes only
*/
void cblas_ssymv(CBLAS_LAYOUT layout, CBLAS_UPLO Uplo,
......@@ -352,7 +352,7 @@ void cblas_dspr2(CBLAS_LAYOUT layout, CBLAS_UPLO Uplo,
const int incX, const double *Y, const int incY, double *A);
/*
/*
* Routines with C and Z prefixes only
*/
void cblas_chemv(CBLAS_LAYOUT layout, CBLAS_UPLO Uplo,
......@@ -423,7 +423,7 @@ void cblas_zhpr2(CBLAS_LAYOUT layout, CBLAS_UPLO Uplo, const int N,
* ===========================================================================
*/
/*
/*
* Routines with standard 4 prefixes (S, D, C, Z)
*/
void cblas_sgemm(CBLAS_LAYOUT layout, CBLAS_TRANSPOSE TransA,
......@@ -547,7 +547,7 @@ void cblas_ztrsm(CBLAS_LAYOUT layout, CBLAS_SIDE Side,
void *B, const int ldb);
/*
/*
* Routines with prefixes C and Z only
*/
void cblas_chemm(CBLAS_LAYOUT layout, CBLAS_SIDE Side,
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -30,16 +30,22 @@
*
* Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
*/
template<class Matrix,class Scalar,class Primary=Scalar>
template<class Matrix,
class Scalar,
class Primary=Scalar,
class KI=LapackKernI /* Kernel Interface Concept */
>
struct Arnoldi_QRInc{
static ArnReturn Exec_Procedure(Matrix& A,Block<Scalar,Primary>& X0,
Block<Scalar,Primary>& B,
static ArnReturn Exec_Procedure(Matrix& A,
Block<Scalar,Primary,KI>& X0,
Block<Scalar,Primary,KI>& B,
int MaxKSize,
std::vector<Primary>& epsilon,
Logger<Primary>& log,
DRParam<Scalar,Primary>& DR_struct,
ArnOrtho ArnChoice,
OrthoChoice ortho = OrthoChoice::MGS){
OrthoChoice ortho = OrthoChoice::MGS)
{
std::cout<<"################# Arnoldi with QR Inc ################\n";
std::cout<<"################# Mat Vect product scheduled "<<MaxKSize<<"\n";
......@@ -54,18 +60,18 @@ struct Arnoldi_QRInc{
//Storage needed for computating
HessQR<Scalar,Primary> Hess{MaxKSize,SizeBlock};
Block<Scalar,Primary> Yj(SizeBlock,MaxKSize+SizeBlock);
Block<Scalar,Primary> LS;
Block<Scalar,Primary,KI> Yj(SizeBlock,MaxKSize+SizeBlock);
Block<Scalar,Primary,KI> LS;
//Temporary
DR_struct.base.reset();
Base<Scalar,Primary> & base = DR_struct.base;
Base<Scalar,Primary,KI> & base = DR_struct.base;
Block<Scalar,Primary> R0(SizeBlock,dim);
Block<Scalar,Primary,KI> R0(SizeBlock,dim);
//Fill matrix R0 with each vector r0_i = B[i] - A * X0[i]
A.MatBlockVect(X0,R0); //R0 = A*X0
for(int i=0 ; i<SizeBlock ; ++i){ //R0 = B - R0
for(int j=0 ; j<dim ; ++j){
for(int i=0 ; i<SizeBlock ; ++i) { //R0 = B - R0
for(int j=0 ; j<dim ; ++j) {
R0.getPtr(i)[j] = B.getPtr(i)[j] - R0.getPtr(i)[j];
}
}
......@@ -73,7 +79,7 @@ struct Arnoldi_QRInc{
std::pair<Primary,Primary> minmaxR0 = R0.getMinMaxNorm(A);
double timestamp = Logger<Primary>::GetTime();
//Init Log with norm of R0 (only if log is empty)
if(0 == log.getNbIteLogged()){
if (0 == log.getNbIteLogged()) {
log.add_ite(-1,B.getSizeBlock(),
minmaxR0.first,minmaxR0.second,timestamp);
}
......@@ -81,8 +87,8 @@ struct Arnoldi_QRInc{
//V1 will be first block of base, and R1, will be used for least
//square pb, in criteria verification
Block<Scalar,Primary> V1(SizeBlock,dim);
Block<Scalar,Primary> R1(SizeBlock,SizeBlock);
Block<Scalar,Primary,KI> V1(SizeBlock,dim);
Block<Scalar,Primary,KI> R1(SizeBlock,SizeBlock);
//Q-R Decomposition of (R0)
R0.ComputeQRFacto_own(V1,R1,A);
//First block of base added
......@@ -95,24 +101,24 @@ struct Arnoldi_QRInc{
bool convergence = false;
int j = 0;
int ite_meter = 0;
while(j<MaxKSize){
while(j<MaxKSize) {
ite_meter++;
timestamp = Logger<Primary>::GetTime();
std::cout<<"\t\tStart of Iteration\n";
//Compute block Wj
Block<Scalar,Primary> Wj(SizeBlock,dim);
Block<Scalar,Primary,KI> Wj(SizeBlock,dim);
//Where to write in hess
Scalar * toWrite = Hess.getNewCol(Wj.getSizeBlock());
if(nullptr == toWrite){ //No more room for this iteration
if (nullptr == toWrite) { //No more room for this iteration
break;
}
if(A.useRightPreCond()){//Wj = A*preCond*V[j] if preCond
Block<Scalar,Primary> temp(SizeBlock,dim);
if (A.useRightPreCond()) {//Wj = A*preCond*V[j] if preCond
Block<Scalar,Primary,KI> temp(SizeBlock,dim);
A.preCondBaseProduct(base.getPtr(base.getCurrentBlock()-1),temp);
A.MatBlockVect(temp,Wj);
}else{//Wj = A*V[j] without preCond
} else {//Wj = A*V[j] without preCond
A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),
Wj.getSizeBlock(),Wj.getPtr());
}
......@@ -124,15 +130,15 @@ struct Arnoldi_QRInc{
//Check if Y need to be computed:
std::tuple<Primary,Primary,bool> MinMax = Hess.CheckRes(epsilon);
if(std::get<2>(MinMax)){
if (std::get<2>(MinMax)) {
Hess.ComputeCurrentSol(Yj,LS);
if(CheckCurrentSolution(A,B,X0,base,Yj,Hess.getNbColUsed(),epsilon)){
if (CheckCurrentSolution(A,B,X0,base,Yj,Hess.getNbColUsed(),epsilon)) {
//Adding last ite
log.add_ite(j,SizeBlock,std::get<0>(MinMax),
std::get<1>(MinMax),timestamp);
convergence = true;
break;
}else{
} else {
std::cout<<"Residual is small enough, but not the real's one \n";
}
}
......@@ -141,19 +147,19 @@ struct Arnoldi_QRInc{
log.add_ite(j,SizeBlock,std::get<0>(MinMax),std::get<1>(MinMax),timestamp);
j += SizeBlock;
}
if(!convergence){ //No convergence, we need to write current sol into X0
if (!convergence) { //No convergence, we need to write current sol into X0
//First, call the triangular solver to get the solution
Hess.ComputeCurrentSol(Yj,LS);
Block<Scalar,Primary> Sol{SizeBlock,dim};
if(A.useRightPreCond()){
Block<Scalar,Primary> temp(SizeBlock,dim);
Block<Scalar,Primary,KI> Sol{SizeBlock,dim};
if (A.useRightPreCond()) {
Block<Scalar,Primary,KI> temp(SizeBlock,dim);
base.ComputeProduct(Yj,temp);
A.preCondBlockVect(temp,Sol);
}else{//No pre Cond used
} else {//No pre Cond used
base.ComputeProduct(Yj,Sol);
}
for(int i=0 ; i<SizeBlock ; ++i){
for(int k=0 ; k<dim ; ++k){
for(int i=0 ; i<SizeBlock ; ++i) {
for(int k=0 ; k<dim ; ++k) {
X0.getPtr(i)[k] += Sol.getPtr(i)[k];
}
}
......
......@@ -28,13 +28,23 @@
*
* Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
*/
template<typename Arn,class Matrix,class Scalar,class Primary>
int BGMRes(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
int totalMVP, int MaxKSize,
Block<Scalar,Primary>& Sol,Logger<Primary>& log,
std::vector<Primary> & epsilon,
OrthoChoice ortho = OrthoChoice::MGS,
ArnOrtho arnChoice = ArnOrtho::BLOCK){
template<typename Arn,
class Matrix,
class Scalar,
class Primary,
class KI=LapackKernI /* Kernel Interface Concept */
>
int BGMRes(Matrix& A,
Block<Scalar,Primary,KI>& B,
Block<Scalar,Primary,KI>& X0,
int totalMVP,
int MaxKSize,
Block<Scalar,Primary,KI>& Sol,
Logger<Primary>& log,
std::vector<Primary>& epsilon,
OrthoChoice ortho = OrthoChoice::MGS,
ArnOrtho arnChoice = ArnOrtho::BLOCK)
{
//Cste Size
int SizeBlock = B.getSizeBlock();
......@@ -48,15 +58,15 @@ int BGMRes(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
std::cout<<"Maximum nb of Mat Vect product:\t"<<totalMVP<<"\n";
std::cout<<"Max Size Krylov space before restart :\t"<<MaxKSize<<"\n";
std::cout<<"Orthogonalization Scheme : "<<ortho<<"\n";
if(epsilon.size() == 1)
if (epsilon.size() == 1)
std::cout<<"Targeted Tolerance : "<<epsilon.front()<<"\n";
else
std::cout<<"Multi precision is being used\n";
std::cout<<"#######################################################\n";
if(A.useRightPreCond()){
if (A.useRightPreCond()) {
std::cout<<"\nRight Pre Conditionner used\n\n";
}else{
} else {
std::cout<<"\nNO - Right Pre Conditionner used\n\n";
}
......@@ -77,18 +87,18 @@ int BGMRes(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
int nbIteDone = 0;
int KrylovSpaceSpanned = 0;
ArnReturn ret{0,0,false};
while(KrylovSpaceSpanned < totalMVP && !ret.hasConverged){
while(KrylovSpaceSpanned < totalMVP && !ret.hasConverged) {
//Compute nb of mat vect product to give to Arnoldi procedure
int SizeToSpan = std::min(MaxKSize,totalMVP-KrylovSpaceSpanned);
if(SizeToSpan < SizeBlock){
if (SizeToSpan < SizeBlock) {
SizeToSpan = SizeBlock;
}
std::cout<<"\n################## Restart number "<<nb_restart
<<" #################\n";
HessStorage<Scalar> hessSave;
Base<Scalar,Primary> base{B.getSizeBlock(),SizeToSpan,dim};
Block<Scalar,Primary> LS;
Base<Scalar,Primary,KI> base{B.getSizeBlock(),SizeToSpan,dim};
Block<Scalar,Primary,KI> LS;
//No DR here, so empty struct with flag to false
DRParam<Scalar,Primary> DR_struct(false,hessSave,base,0,Scalar(0),LS);
......
......@@ -32,14 +32,25 @@
* Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
* Restarted Block GMRes with Deflation of eigenvalues (Morgan)
*/
template<typename Arn,class Matrix,class Scalar,class Primary>
int BGMResDR(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
int totalMVP, int MaxKSize,
Block<Scalar,Primary>& Sol,Logger<Primary>& log,
std::vector<Primary> & epsilon,
int nb_ev, Scalar target,
OrthoChoice ortho = OrthoChoice::MGS,
ArnOrtho arnChoice = ArnOrtho::BLOCK){
template<typename Arn,
class Matrix,
class Scalar,
class Primary,
class KI=LapackKernI /* Kernel Interface Concept */
>
int BGMResDR(Matrix& A,
Block<Scalar,Primary,KI>& B,
Block<Scalar,Primary,KI>& X0,
int totalMVP,
int MaxKSize,
Block<Scalar,Primary,KI>& Sol,
Logger<Primary>& log,
std::vector<Primary>& epsilon,
int nb_ev,
Scalar target,
OrthoChoice ortho = OrthoChoice::MGS,
ArnOrtho arnChoice = ArnOrtho::BLOCK)
{
//Cste Size
int SizeBlock = B.getSizeBlock();
......@@ -53,7 +64,7 @@ int BGMResDR(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
std::cout<<"Maximum nb of Matrix Vect product:\t"<<totalMVP<<"\n";
std::cout<<"Max Size Krylov space before restart :\t"<<MaxKSize<<"\n";
std::cout<<"Orthogonalization Scheme : "<<ortho<<"\n";
if(epsilon.size() == 1)
if (epsilon.size() == 1)
std::cout<<"Targeted Tolerance : "<<epsilon.front()<<"\n";
else
std::cout<<"Multi precision is being used\n";
......@@ -62,9 +73,9 @@ int BGMResDR(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
std::cout<<"#######################################################\n";
if(A.useRightPreCond()){
if (A.useRightPreCond()) {
std::cout<<"\n Right Pre Conditionner used\n\n";
}else{
} else {
std::cout<<"\n NO - Right Pre Conditionner used\n\n";
}
......@@ -88,19 +99,18 @@ int BGMResDR(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
//Here we define all the structure that need to be passed from 1 call to
//arnoldi to the following one during restart
Base<Scalar,Primary> base{B.getSizeBlock(),MaxKSize+SizeBlock+1,dim}; //Base
Base<Scalar,Primary,KI> base{B.getSizeBlock(),MaxKSize+SizeBlock+1,dim}; //Base
HessStorage<Scalar> hessSave; //Struct holding the hessenberg
Block<Scalar,Primary> LS; //LS <- (H*Y-R1)
Block<Scalar,Primary,KI> LS; //LS <- (H*Y-R1)
//Struct to hold them all
DRParam<Scalar,Primary> DR_struct(true,hessSave,base,nb_ev,target,LS);
while(KrylovSpaceSpanned < totalMVP && !ret.hasConverged){ //Main loop
while(KrylovSpaceSpanned < totalMVP && !ret.hasConverged) { //Main loop
//Compute nb of mat vect product to give to Arnoldi procedure
int SizeToSpan = std::min(MaxKSize,totalMVP-KrylovSpaceSpanned);
//Test if SizeTospan is big enough to hold 1 block + Deflation directions
if(SizeToSpan < (SizeBlock+nb_ev)){
if (SizeToSpan < (SizeBlock+nb_ev)) {
std::cout<<"No restart because we reach the maximum size allowed\n";
break;
}
......@@ -119,7 +129,6 @@ int BGMResDR(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
KrylovSpaceSpanned += ret.sizeKSpace;
nb_restart++;
std::cout<<"######################################################\n";
std::cout<<"\t\t Restart Done : "<<nb_restart
<<"\n";
......
This diff is collapsed.
This diff is collapsed.
......@@ -8,58 +8,50 @@
* concatenated.
* Thus a cursor to know where P ends and W starts is needed.
*/
template<typename Scalar,typename Primary>
class BlockWP : public Block<Scalar,Primary>{
template<typename Scalar,
typename Primary,
typename KI=LapackKernI /* KernelInterface Concept */
>
class BlockWP : public Block<Scalar,Primary, KI>
{
int cursor;
using Parent=Block<Scalar,Primary>;
typedef Block<Scalar,Primary,KI> Parent;
typedef Block<Scalar,Primary,KI> BlockSPK;
public:
BlockWP(int sizeBlock, int sizeVector) : Parent(sizeBlock,sizeVector),cursor(0){
BlockWP(int sizeBlock, int sizeVector):
Parent(sizeBlock,sizeVector),
cursor(0)
{
}
~BlockWP(){
}
Scalar * getW(){
return Parent::getPtr(cursor);
}
Scalar * getP(){
return Parent::getPtr();
}
int getSizeP(){
return cursor;
}
int getSizeW(){
return Parent::getSizeBlock()-cursor;
}
void incSizeP(int inc){
cursor += inc;
}
Scalar* getW() { return Parent::getPtr(cursor); }
Scalar* getP() { return Parent::getPtr(); }
int getSizeP() const { return cursor; }
int getSizeW() const { return Parent::getSizeBlock()-cursor; }
void incSizeP(int inc) { cursor += inc; }
/**
* @brief Compute C = P^{H} * W and then W = W-C*P
*
* @param
*/
void OrthoWP(Scalar * toWrite, int ldToWrite){
if(cursor>0){
void OrthoWP(Scalar *toWrite, int ldToWrite) {
if (cursor>0) {
std::cout<<"Ortho of "<<getSizeW()<<" vectors in W against "<<getSizeP()
<<" vectors in P\n";
//Generate coef inside C
call_Tgemm<>(getSizeP(),getSizeW(),Parent::getLeadingDim(),
getP(),Parent::getLeadingDim(),
getW(),Parent::getLeadingDim(),
toWrite,ldToWrite,Scalar(1),Scalar(0));
//Modify W part relatively to C and P
//Wj=(-1) * Pj * toWrite + Wj;
call_gemm<>(Parent::getLeadingDim(),getSizeW(),getSizeP(),
KI::Tgemm(getSizeP(),getSizeW(),Parent::getLeadingDim(),
getP(),Parent::getLeadingDim(),
toWrite,ldToWrite,
getW(),Parent::getLeadingDim(),
Scalar(-1),Scalar(1));
toWrite,ldToWrite,Scalar(1),Scalar(0));
//Modify W part relatively to C and P
//Wj=(-1) * Pj * toWrite + Wj;
KI::gemm(Parent::getLeadingDim(),getSizeW(),getSizeP(),
getP(),Parent::getLeadingDim(),
toWrite,ldToWrite,
getW(),Parent::getLeadingDim(),
Scalar(-1),Scalar(1));
std::cout<<"WP :: Ortho done\n";
}
}
......@@ -68,31 +60,30 @@ public:
* @brief Compute QR facto of W part. W will be overwritten by
* Q. R part generated will go inside input toWrite ptr
*/
void QRofW(Scalar * toWriteR, int ldR){
void QRofW(Scalar *toWriteR, int ldR) {
//Create temporary tau
Scalar * tau = new Scalar[getSizeW()];
memset(tau,0,sizeof(Scalar)*getSizeW());
typename KI:: template TauHolder<Scalar> tau;
//3 steps::
//1 - Call qr facto on W bloc
int geqrf = callLapack_geqrf<>(Parent::getLeadingDim(),getSizeW(),getW(),
Parent::getLeadingDim(),tau);
if(geqrf != 0){
int geqrf = KI::geqrf(Parent::getLeadingDim(),getSizeW(),getW(),
Parent::getLeadingDim(), tau);
if (geqrf != 0) {
std::cout<<"Problem on QR facto\n";
}
//2 - get the R part and write it in input
for(int i=0 ; i<getSizeW() ; ++i){
for(int i=0 ; i<getSizeW() ; ++i) {
memcpy(&toWriteR[i*ldR],&getW()[i*Parent::getLeadingDim()],
(i+1)*sizeof(Scalar));
}
//3 - Call orgqr on W block
int orgqr = callLapack_orgqr(Parent::getLeadingDim(),getSizeW(),
getW(),Parent::getLeadingDim(),tau);
int orgqr = KI::orgqr(Parent::getLeadingDim(),getSizeW(),
getW(),Parent::getLeadingDim(),tau);
//check