1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
|
#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
|