#ifndef GTHREADET_C
#define GTHREADET_C GTHREADET_C
#include"row_dispenser.c"
#include<time.h>
/*
* using the GNU multiprecision library because there are no specified limits to the matrix elements (even if there were the det is going to be huge regardless)
* and because it has handy rational numbers and input functions
*/
#include"gmp.h"
/*
* posix threads for maximum portability :>
*/
#include<pthread.h>
mpq_t* threadet(mpq_t **matrix,size_t size,size_t number_of_threads,char verbose);
void* dethread(void* thr);
void calc_row(mpq_t **matrix,size_t size,size_t target,size_t source);
/*
* here the process creates and dispatches number_of_threads threads
*/
mpq_t* threadet(mpq_t **matrix,size_t size,size_t number_of_threads,char verbose)
{
size_t i;
size_t j;
size_t k;
/* shared thread parameters */
struct Pair_Dispenser thr;
/* swap */
mpq_t *hold;
/* array of the threads we are executing. We are going to join these every column of the matrix cleaned because it is a very convinient sync.*/
pthread_t* threds;
/* result is accumulated here */
mpq_t *res;
/* placeholder for thread return values. They all return NULL but join() requires a pointer to store return */
void **temp;
/* we can't clean a column with a zero so we swap a impotent column with a column with a non zero value ( if there is such) */
char sign;
/* timers REALTIME CLOCK could be fooled by an administrator changing the system time, but is the only guaranteed clock in the posix standart*/
struct timespec pass_start,pass_finish;
struct timespec start,finish;
if(verbose==1)
{
clock_gettime(CLOCK_REALTIME,&start);
}
threds=malloc(sizeof(pthread_t*) * number_of_threads);
res=malloc(sizeof(mpq_t));
temp = malloc(sizeof(void*));
Pair_Dispenser_Init(&thr,size,matrix,verbose);
mpq_init(*res);
thr.current.target_row=1;
sign=1;
for(i=0;i<size;++i)
{
/* if the row can't be used to clean the column search for one that can */
if(mpq_sgn(matrix[i][i])==0)
{
for(j=i+1;j<=size && mpq_sgn(matrix[j][i])==0;++j);
if(j>size)
{
mpq_set_d(*res,0);
free(threds);
free(temp);
if(verbose == 1)
{
clock_gettime(CLOCK_REALTIME,&finish);
finish.tv_sec-=start.tv_sec;
finish.tv_nsec-=start.tv_nsec;
printf("TIME: %f\n",(double)(finish.tv_sec + 0.000000001*finish.tv_nsec));
}
return res;
}else
{
hold = matrix[j];
matrix[j] = matrix[i];
matrix[i] = hold;
/* change the sign approprietly */
sign = (sign + (j-i)%2)%2;
}
}
if(verbose==1)
{
clock_gettime(CLOCK_REALTIME,&pass_start);
}
/* clean the i-th row this is one pass*/
k=((number_of_threads<size-i)?number_of_threads:size-i);
for(j=0;j<k;++j)
{
pthread_create(threds+j,NULL,dethread,&thr);
}
/* wait for all the threads to finish calculating this is the second 1/2 that could slow down parallel execution*/
/* critical part */
for(j=0;j<k;++j)
{
pthread_join(threds[j],temp);
}
/* critical part */
if(verbose == 1)
{
clock_gettime(CLOCK_REALTIME,&pass_finish);
pass_finish.tv_sec-=pass_start.tv_sec;
pass_finish.tv_nsec-=pass_start.tv_nsec;
printf("A pass has finished, taking: %f seconds\n\n",(double)(pass_finish.tv_sec + 0.000000001*pass_finish.tv_nsec));
}
thr.current.target_row = (++thr.current.source_row) + 1;
}
mpq_set_d(*res,((sign == 1)?1:-1));
for(i=0;i<=size;++i)
{
mpq_mul(*res,*res,matrix[i][i]);
}
free(threds);
free(temp);
if(verbose == 1)
{
clock_gettime(CLOCK_REALTIME,&finish);
finish.tv_sec-=start.tv_sec;
finish.tv_nsec-=start.tv_nsec;
printf("TIME: %f\n",(double)(finish.tv_sec + 0.000000001*finish.tv_nsec));
}
return res;
}
void* dethread(void* thr)
{
#define THR ((struct Pair_Dispenser*)thr)
clock_t start,finish;
if(THR->verbose == 1)
{
start=clock();
}
struct row_pair temp;
/* while you can get a valid pair - calculate */
while((temp=Pair_Dispenser_Get_Pair(THR)).target_row <= THR->max)
{
calc_row(THR->matrix,THR->max,temp.target_row,temp.source_row);
}
if(THR->verbose == 1)
{
finish=clock();
printf("A thread has finished, taking: %f seconds\n",((double)(finish-start))/CLOCKS_PER_SEC);
}
return NULL;
#undef THR
}
void calc_row(mpq_t **matrix,size_t size,size_t target,size_t source)
{
mpq_t mod;
mpq_t temp1;
mpq_init(mod);
mpq_init(temp1);
mpq_div(mod,matrix[target][source],matrix[source][source]);
size_t k;
for(k=source;k<=size;++k)
{
mpq_mul(temp1,matrix[source][k],mod);
mpq_sub(matrix[target][k],matrix[target][k],temp1);
}
}
#endif //#ifndef GTHREADET_C