54 #define printfRootInternal(...) { if( hit_Rank == 0 ) { printf(__VA_ARGS__); fflush(stdout); }}
55 #define printfRoot(...) printfRootInternal(__VA_ARGS__)
173 printf(
"%.2lf\t",m.
m[i][j]);
187 printf(
"%.2lf\t",v.
v[i]);
197 #define clearV(vv) clearVInternal((vv).v)
198 #define clearVInternal(v) \
199 ((v)[0] = (v)[1] = (v)[2] = 0.0 )
204 #define cpyV(r,a) cpyVInternal((r).v,(a).v)
205 #define cpyVInternal(r,a) \
206 {(r)[0] = (a)[0]; (r)[1] = (a)[1]; (r)[2] = (a)[2]; }
211 #define addV(r,a,b) addVInternal((r).v,(a).v,(b).v)
212 #define addVInternal(r,a,b) \
213 {(r)[0] = (a)[0] + (b)[0]; (r)[1] = (a)[1] + (b)[1]; (r)[2] = (a)[2] + (b)[2]; }
219 #define subV(r,a,b) subVInternal((r).v,(a).v,(b).v)
220 #define subVInternal(r,a,b) \
221 {(r)[0] = (a)[0] - (b)[0]; (r)[1] = (a)[1] - (b)[1]; (r)[2] = (a)[2] - (b)[2]; }
226 #define norm(vv) normInternal((vv).v)
227 #define normInternal(v) \
228 (sqrt( pow((v)[0],2) + pow((v)[1],2) + pow((v)[2],2) ))
233 #define multSV(r,s,vv) multSVInternal((r).v,s,(vv).v)
234 #define multSVInternal(r,s,v) \
235 {(r)[0] = (s) * (v)[0]; (r)[1] = (s) * (v)[1]; (r)[2] = (s) * (v)[2]; }
240 #define clearM(mm) clearMInternal((mm).m)
241 #define clearMInternal(m) \
242 ((m)[0][0] = (m)[0][1] = (m)[0][2] = \
243 (m)[1][0] = (m)[1][1] = (m)[1][2] = \
244 (m)[2][0] = (m)[2][1] = (m)[2][2] = 0.0 )
256 #define det(mm) detInternal((mm).m)
257 #define detInternal(m) \
258 ( (m)[0][0] * (m)[1][1] * (m)[2][2] \
259 + (m)[0][1] * (m)[1][2] * (m)[2][0] \
260 + (m)[0][2] * (m)[1][0] * (m)[2][1] \
261 - (m)[0][0] * (m)[1][2] * (m)[2][1] \
262 - (m)[0][1] * (m)[1][0] * (m)[2][2] \
263 - (m)[0][2] * (m)[1][1] * (m)[2][0] )
268 #define init_matrix(mm,a,b,c,d,e,f,g,h,i) init_matrixInternal((mm).m,a,b,c,d,e,f,g,h,i)
269 #define init_matrixInternal(m,a,b,c,d,e,f,g,h,i) \
270 {(m)[0][0] = a; (m)[0][1] = b; (m)[0][2] = c; \
271 (m)[1][0] = d; (m)[1][1] = e; (m)[1][2] = f; \
272 (m)[2][0] = g; (m)[2][1] = h; (m)[2][2] = i; }
277 #define multSM(r,s,mm) multSMInternal((r).m,s,(mm).m)
278 #define multSMInternal(r,s,m) \
279 {(r)[0][0] = (s) * (m)[0][0]; (r)[0][1] = (s) * (m)[0][1]; (r)[0][2] = (s) * (m)[0][2]; \
280 (r)[1][0] = (s) * (m)[1][0]; (r)[1][1] = (s) * (m)[1][1]; (r)[1][2] = (s) * (m)[1][2]; \
281 (r)[2][0] = (s) * (m)[2][0]; (r)[2][1] = (s) * (m)[2][1]; (r)[2][2] = (s) * (m)[2][2]; }
286 #define multSI(r,s) multSIInternal((r).m,s)
287 #define multSIInternal(r,s) \
288 {(r)[0][0] = (s); (r)[0][1] = 0; (r)[0][2] = 0; \
289 (r)[1][0] = 0; (r)[1][1] = (s); (r)[1][2] = 0; \
290 (r)[2][0] = 0; (r)[2][1] = 0; (r)[2][2] = (s); }
296 #define addM(r,a,b) addMInternal((r).m,(a).m,(b).m)
297 #define addMInternal(r,a,b) \
298 {(r)[0][0] = (a)[0][0] + (b)[0][0]; (r)[0][1] = (a)[0][1] + (b)[0][1]; (r)[0][2] = (a)[0][2] + (b)[0][2]; \
299 (r)[1][0] = (a)[1][0] + (b)[1][0]; (r)[1][1] = (a)[1][1] + (b)[1][1]; (r)[1][2] = (a)[1][2] + (b)[1][2]; \
300 (r)[2][0] = (a)[2][0] + (b)[2][0]; (r)[2][1] = (a)[2][1] + (b)[2][1]; (r)[2][2] = (a)[2][2] + (b)[2][2]; }
306 #define multVVt(r,a,b) multVVtInternal((r).m,(a).v,(b).v)
307 #define multVVtInternal(r,a,b) \
308 {(r)[0][0] = (a)[0] * (b)[0]; (r)[0][1] = (a)[0] * (b)[1]; (r)[0][2] = (a)[0] * (b)[2]; \
309 (r)[1][0] = (a)[1] * (b)[0]; (r)[1][1] = (a)[1] * (b)[1]; (r)[1][2] = (a)[1] * (b)[2]; \
310 (r)[2][0] = (a)[2] * (b)[0]; (r)[2][1] = (a)[2] * (b)[1]; (r)[2][2] = (a)[2] * (b)[2]; }
316 #define multMV(r,mm,vv) multMVInternal((r).v,(mm).m,(vv).v)
317 #define multMVInternal(r,m,v) \
318 {(r)[0] = (m)[0][0] * (v)[0] + (m)[0][1] * (v)[1] + (m)[0][2] * (v)[2]; \
319 (r)[1] = (m)[1][0] * (v)[0] + (m)[1][1] * (v)[1] + (m)[1][2] * (v)[2]; \
320 (r)[2] = (m)[2][0] * (v)[0] + (m)[2][1] * (v)[1] + (m)[2][2] * (v)[2]; }
326 #define inv(r,mm) invInternal((r).m,(mm).m)
327 #define invInternal(r,m) \
328 {(r)[0][0] = (m)[1][1] * (m)[2][2] - (m)[1][2] * (m)[2][1]; \
329 (r)[0][1] = - (m)[0][1] * (m)[2][2] + (m)[0][2] * (m)[2][1]; \
330 (r)[0][2] = (m)[0][1] * (m)[1][2] - (m)[0][2] * (m)[1][1]; \
331 (r)[1][0] = - (m)[1][0] * (m)[2][2] + (m)[1][2] * (m)[2][0]; \
332 (r)[1][1] = (m)[0][0] * (m)[2][2] - (m)[0][2] * (m)[2][0]; \
333 (r)[1][2] = - (m)[0][0] * (m)[1][2] + (m)[0][2] * (m)[1][0]; \
334 (r)[2][0] = (m)[1][0] * (m)[2][1] - (m)[1][1] * (m)[2][0]; \
335 (r)[2][1] = - (m)[0][0] * (m)[2][1] + (m)[0][1] * (m)[2][0]; \
336 (r)[2][2] = (m)[0][0] * (m)[1][1] - (m)[0][1] * (m)[1][0]; \
337 multSMInternal(r,1/detInternal(m),r); }
342 #define multMM(r,a,b) multMMInternal((r).m,(a).m,(b).m)
343 #define multMMInternal(r,a,b) \
344 {(r)[0][0] = (a)[0][0] * (b)[0][0] + (a)[0][1] * (b)[1][0] + (a)[0][2] * (b)[2][0]; \
345 (r)[0][1] = (a)[0][0] * (b)[0][1] + (a)[0][1] * (b)[1][1] + (a)[0][2] * (b)[2][1]; \
346 (r)[0][2] = (a)[0][0] * (b)[0][2] + (a)[0][1] * (b)[1][2] + (a)[0][2] * (b)[2][2]; \
347 (r)[1][0] = (a)[1][0] * (b)[0][0] + (a)[1][1] * (b)[1][0] + (a)[1][2] * (b)[2][0]; \
348 (r)[1][1] = (a)[1][0] * (b)[0][1] + (a)[1][1] * (b)[1][1] + (a)[1][2] * (b)[2][1]; \
349 (r)[1][2] = (a)[1][0] * (b)[0][2] + (a)[1][1] * (b)[1][2] + (a)[1][2] * (b)[2][2]; \
350 (r)[2][0] = (a)[2][0] * (b)[0][0] + (a)[2][1] * (b)[1][0] + (a)[2][2] * (b)[2][0]; \
351 (r)[2][1] = (a)[2][0] * (b)[0][1] + (a)[2][1] * (b)[1][1] + (a)[2][2] * (b)[2][1]; \
352 (r)[2][2] = (a)[2][0] * (b)[0][2] + (a)[2][1] * (b)[1][2] + (a)[2][2] * (b)[2][2]; }
380 for(iter=0; iter<
ITER_SJ; iter++){
383 if( A.
m[i][i] == 0 )
continue;
388 if(i != j) x_1.
v[i] += A.
m[i][j] * x.
v[j];
390 x_1.
v[i] = (b.
v[i] - x_1.
v[i]) / A.
m[i][i];
414 double detA =
det(A);
470 printf(
"%s [-n NEWTON METHOD ITERATIONS] [-j JACOBI METHOD ITERATIONS] FILE \n",name);
472 printf(
" -n NEWTON METHOD ITERATIONS number of iterations (default 10)\n");
473 printf(
" -j JACOBI METHOD ITERATIONS number of iterations (default 100)\n");
474 printf(
" FILE input graph file\n");
492 if(!(p = strstr(str, orig)))
return NULL;
494 strncpy(buffer, str, (
size_t) (p-str));
495 buffer[p-str] =
'\0';
497 sprintf(buffer+(p-str),
"%s%s", rep, p+strlen(orig));
509 FILE * file = fopen(name,
"r");
513 printf(
"File not exists\n");
519 char c = (char) getc(file);
523 while( c !=
'\n' ) c = (char) getc(file);
530 res = fscanf(file,
"%d %d\n",&n,&d);
533 if(n != n_shape) printf(
"Coordinates: Number of vertices do not match %d\n",n); n =
hit_min(n,n_shape);
534 if(d != D) printf(
"Coordinates: Number of dimensions do not match %d\n",d); d =
hit_min(d,D);
544 res = fscanf(file,
"%lf\n",&number);
555 if ((random() % 100) <
P_FIXED){
581 while ((c = (
char) getopt (argc, argv,
"hn:j:")) != -1)
584 sscanf(optarg,
"%d",&
iter1);
587 sscanf(optarg,
"%d",&
iter2);
601 char * graph_file = argv[optind];
619 char * coord_file =
replace_str(graph_file,
".rb",
"_coord.mtx");
632 printf(
"# N fixed: %d\n",nfixed);
649 Node n0 = {{{0,0,0}},1};
651 Node n1 = {{{100,0,0}},0};
653 Node n2 = {{{200,0,0}},0};
655 Node n3 = {{{300,0,0}},1};
657 Node n4 = {{{0,100,200}},1};
659 Node n5 = {{{100,100,0}},0};
661 Node n6 = {{{200,100,0}},0};
663 Node n7 = {{{300,100,0}},1};
679 HitTile_double
norm, sum_norm;
738 int main(
int argc,
char ** argv) {
750 printf(
"There is only one processor\n");
768 HitTile_int vertices;
785 lay =
hit_layout(plug_layBlocks,topo,global_shape);
831 printfRoot(
"# Initial gradient norm: %lf\n",ignorm);
835 for(i=0;i<
iter1;i++){
853 for(j=0;j<
iter2;j++){
864 memcpy(
e.data,
new_e.data,
sizeof(
Vector)* (
size_t) Nvertices);
873 if (!current->
fixed){
885 printfRoot(
"# Final gradient norm: %lf\n",gnorm);
940 double scalar = -
K * ((
L - n) / n);
968 double scalar_a, scalar_b;
981 scalar_a = -((
L-n)/(n));
983 scalar_b = (
L/(pow(n,3)));
1010 int nghb_index = edge;
1013 Vector gij =
g(vertex,nghb_index);
1017 Matrix Wij =
W(vertex,nghb_index);
1025 Matrix Hij =
W(nghb_index,vertex);
1062 int nghb_index = edge;
1063 if(nghb_index == vertex)
continue;
1067 if(nj.
fixed)
continue;
1088 new_ei =
solve(Hii,sum);
#define hit_shape(nd,...)
void print_help(char *name)
#define hit_layShape(lay)
#define hit_tileNewType(baseType)
void hit_comOpSumDouble(void *, void *, int *, HitType *)
HitCom hit_comBroadcastSelect(HitLayout lay, HitRanks root, const void *tile, HitShape selection, int modeSelect, HitType baseType)
#define hit_tileElemAt(var, ndims,...)
#define HIT_LAYOUT_NULL_STATIC
void hit_comFree(HitCom issue)
char * replace_str(char *str, const char *orig, const char *rep)
void hit_comDo(HitCom *issue)
#define hit_tileDomainAlloc(newVarP, baseType, numDims,...)
void solve_system_iter(int vertex)
#define hit_topology(name,...)
#define hit_layShapeOther(lay, ranks)
#define hit_fileHBVertices(hbfile)
#define HIT_COM_ARRAYCOORDS
#define HIT_PAT_UNORDERED
void hit_patternAdd(HitPattern *pattern, HitCom comm)
Vector solve(Matrix A, Vector b)
#define hit_tileDomainShapeAlloc(newVarP, baseType, shape)
void hit_topFree(HitTopology topo)
#define hit_comOp(function, operation)
Vector solve_jacobi(Matrix A, Vector b)
#define hit_clockWorldReduce(c)
void hit_patternFree(HitPattern *pattern)
void hit_comInit(int *pargc, char **pargv[])
#define hit_layImActive(lay)
void hit_layFree(HitLayout lay)
int main(int argc, char *argv[])
#define multMV(r, mm, vv)
#define hit_comReduce(lay, root, tilePSend, tilePRecv, baseType, operation)
#define hit_comTypeStruct(new_type, Nstruct, n,...)
void init_graph(HitTile_double graph, HitShape shape_global)
#define hit_shapeIterator(var, shape, dim)
int hit_layNumActives(HitLayout lay)
void calculate_GH(int vertex)
#define hit_comBroadcast(lay, root, tile, baseType)
void read_coordinates(char *name, HitShape shape_rc, HitTile_Node graph_rc)
#define hit_layout(name, topo,...)
#define hit_shapeSig(shape, dim)
#define hit_tileFree(var)
#define hit_patternDo(pattern)
#define hit_clockGetSeconds(c)
#define hit_clockStart(c)