diff --git a/Api/src/IbGMResDrEngine.hpp b/Api/src/IbGMResDrEngine.hpp index 626bc45bdae02d6a557b8a9da4d089bc4a669323..6987bbfa085e84f4b1f3e1a3efbc8c6d45308d05 100644 --- a/Api/src/IbGMResDrEngine.hpp +++ b/Api/src/IbGMResDrEngine.hpp @@ -95,7 +95,7 @@ public: log = new Logger<S>(); //Call to BlockGMRes with parameters int res = BlockGMRes<UserMatrix<T>,T,S>(*matrix,B,X_init,max_iter,*sol, - max_iter,*log,tolerance); + restart,*log,tolerance); return res; } diff --git a/Api/tests/testApiGMresUserCase.c b/Api/tests/testApiGMresUserCase.c new file mode 100644 index 0000000000000000000000000000000000000000..39a5ccbb5e3eb6912ffa0560520d338cc45d8f2b --- /dev/null +++ b/Api/tests/testApiGMresUserCase.c @@ -0,0 +1,136 @@ +#include <stdlib.h> +#include <stdio.h> +#include <string.h> +#include "Ib-GMRes-Dr.h" +#include <complex.h> +#include <cblas.h> + +/** + * This test case solve the system provided by G. Sylvand. + * + */ + +//Structure User Env +typedef struct env{ + int nbRHS; +}Env; + + +//Structure for storing the matrix +typedef struct matrix{ + int dim1; //nb line + int dim2; //nb col + double complex * datas; //is an array of dim*dim doubles complex +}Matrix; + + +//Define callbacks +//Matrix block of vector products +void MatxBlockVect(void * mat, int nbVect, void * input, void **output, void * env){ + Matrix * my_mat = (Matrix*)mat; + Env * my_env = (Env*)env; + + int nbLine = my_mat->dim1; + int nbRHS = my_env->nbRHS; + int nbCol = my_mat->dim1; + double complex alpha = 1; + double complex beta = 0; + cblas_zgemm(CblasColMajor,CblasNoTrans,CblasNoTrans, + nbLine,nbRHS,nbCol,&alpha,my_mat->datas,nbCol, + input,nbCol,&beta,*output,nbCol); +} + +void load_datas(Matrix * mat,FILE * fd){ + int ari; + fread(&ari,sizeof(int),1,fd); + //Check : + if(ari != 3){ //Check complex double + printf("EUH ....\n"); + exit(0); + } + //Then + int dim[2]; + fread(dim,sizeof(int),2,fd); + printf("Mat is %dx%d\n",dim[0],dim[1]); + int size[2]; //will contains SizeofEl and 0 + fread(size,sizeof(int),2,fd); + printf("\tSizeOfEl : %d, and Zero :%d\n",size[0],size[1]); + + //Allocate what's needed + mat->datas = malloc(size[0]*dim[0]*dim[1]); + //set dim + mat->dim1 = dim[0]; + mat->dim2 = dim[1]; + + //Fill our structure + fread(mat->datas,size[0],dim[0]*dim[1],fd); +} + + +int main(int ac , char ** av){ + //Open Matrix file + if(ac < 3){ + printf("Need to provide 2 files : Matrice and RHS\n"); + return 0; + } + FILE * fd = fopen(av[1], "r"); + //Create our structure + Matrix * mat = malloc(sizeof(Matrix)); + load_datas(mat,fd); + fclose(fd); + + //Same with RHS : + FILE * fd2 = fopen(av[2], "r"); + Matrix * RHS = malloc(sizeof(Matrix)); + load_datas(RHS,fd2); + fclose(fd2); + + //Create env + Env * userEnv = malloc(sizeof(Env)); + userEnv->nbRHS = RHS->dim2; + + //Create Starting Block + double * X0 = malloc(sizeof(double)* userEnv->nbRHS * mat->dim1); + + //Init library + IbGMResDr_handle handle = Init_ibgmresdr(CMPLX_DBL,mat,mat->dim1,userEnv); + //Set my product Matrix block of vectors + ibgmresdr_set_usergemm(MatxBlockVect,handle); + //Tolerance + double tol = 0.01; + //set param : nb_iter, restart, ptr to tolerance + ibgmresdr_set_parameters((mat->dim1)/userEnv->nbRHS, + ((mat->dim1)/userEnv->nbRHS)/2, + &tol, handle); + + //Solve env + int nb_iter_needed = ibgmresdr_solve(userEnv->nbRHS,RHS->datas,X0,handle); + + //Won't be used (remark : we won't free the memory since it will + //be done by calling dealloc on the library) + double complex * res = ibgmresdr_get_results(handle); + + //Get the convergence history + double * historic = ibgmresdr_get_logs(handle); + + //Write down the convergence + FILE * fd3 = fopen("/home/cyrille/IbGMResDr/matrices_BEM/ConvergenceHisto.txt", + "w"); + fprintf(fd3,"Ite\tMin\tMax\tTime\n"); + for(int i=0 ; i<nb_iter_needed ; ++i){ + fprintf(fd3,"%d\t%e\t%e\t%e\n",i,historic[3*i],historic[3*i+1],historic[3*i+2]); + } + fclose(fd3); + + //Free the handle + ibgmresdr_dealloc(handle); + + + free(X0); + free(userEnv); + free(RHS->datas); + free(RHS); + free(mat->datas); + free(mat); + return 0; +}