Commit ec31a124 authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille

User test Case added (datas from G. Sylvand)

parent 8c81a7dc
......@@ -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;
}
......
#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;
}
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