Hitmap 1.3
 All Data Structures Namespaces Files Functions Variables Typedefs Friends Macros Groups Pages
spring_full_dense.c
Go to the documentation of this file.
1 
12 /*
13  * <license>
14  *
15  * Hitmap v1.2
16  *
17  * This software is provided to enhance knowledge and encourage progress in the scientific
18  * community. It should be used only for research and educational purposes. Any reproduction
19  * or use for commercial purpose, public redistribution, in source or binary forms, with or
20  * without modifications, is NOT ALLOWED without the previous authorization of the copyright
21  * holder. The origin of this software must not be misrepresented; you must not claim that you
22  * wrote the original software. If you use this software for any purpose (e.g. publication),
23  * a reference to the software package and the authors must be included.
24  *
25  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS "AS IS" AND ANY
26  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
27  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
28  * THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
29  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
30  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
31  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
32  * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
33  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
34  *
35  * Copyright (c) 2007-2015, Trasgo Group, Universidad de Valladolid.
36  * All rights reserved.
37  *
38  * More information on http://trasgo.infor.uva.es/
39  *
40  * </license>
41 */
42 
43 #include <stdio.h>
44 #include <stdlib.h>
45 #include <hitmap.h>
46 #include <hit_com.h>
47 #include <unistd.h>
48 #include <string.h>
49 
50 
51 
52 
53 /* OUTPUT FUNCIONS */
54 #define printfRootInternal(...) { if( hit_Rank == 0 ) { printf(__VA_ARGS__); fflush(stdout); }}
55 #define printfRoot(...) printfRootInternal(__VA_ARGS__)
56 
57 
59 #define D (3)
60 
62 #define K (1)
63 
64 #define L (0.001)
65 
67 #define P_FIXED 10
68 
70 #define ITER1 10
71 
72 #define ITER2 100
73 
75 int iter1 = ITER1;
77 int iter2 = ITER2;
78 
82 typedef struct {
83  double v[D];
84 } Vector;
85 
89 typedef struct {
90  double m[D][D];
91 } Matrix;
92 
96 typedef struct {
97  Vector r;
98  int fixed;
99 } Node;
100 
101 
102 
103 // Define the tile_Node.
105 
106 // Define the tile_VectorA.
108 
109 // Define the tile_double
110 hit_tileNewType(double);
111 
112 // Define the tile_Matrix
114 
115 hit_tileNewType(int);
116 
118 
119 
123 
124 
126 HitTile_Vector gradient;
127 
129 HitTile_Matrix hessian;
130 
132 HitTile_Vector e;
133 
135 HitTile_Vector new_e;
136 
139 
141 HitTile_Node graph;
142 
143 
144 
147 
148 
149 
154 void calculate_GH(int vertex);
155 
160 void solve_system_iter(int vertex);
161 
162 
163 #ifdef DEBUG
164 
168 void print_matrix(Matrix m){
169 
170  int i,j;
171  for(i=0;i<D;i++){
172  for(j=0;j<D;j++){
173  printf("%.2lf\t",m.m[i][j]);
174  }
175  printf("\n");
176  }
177 }
178 
183 void print_vector(Vector v){
184 
185  int i;
186  for(i=0;i<D;i++){
187  printf("%.2lf\t",v.v[i]);
188  }
189  printf("\n");
190 }
191 #endif
192 
193 
197 #define clearV(vv) clearVInternal((vv).v)
198 #define clearVInternal(v) \
199  ((v)[0] = (v)[1] = (v)[2] = 0.0 )
200 
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]; }
207 
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]; }
214 
215 
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]; }
222 
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) ))
229 
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]; }
236 
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 )
245 
246 
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] )
264 
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; }
273 
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]; }
282 
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); }
291 
292 
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]; }
301 
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]; }
311 
312 
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]; }
321 
322 
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); }
338 
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]; }
353 
354 
363 
364 // Number of iteration of the jacobi iterative method.
365 #define ITER_SJ 10
366 
367  Vector x;
368 
369  // The new approximation.
370  Vector x_1;
371 
372  // Clear the vectors.
373  clearV(x);
374  clearV(x_1);
375 
376  int iter;
377  int i,j;
378 
379  // Iterative method.
380  for(iter=0; iter<ITER_SJ; iter++){
381  for(i=0;i<D;i++){
382 
383  if( A.m[i][i] == 0 ) continue;
384 
385  // calculate the new value.
386  x_1.v[i] = 0.0;
387  for(j=0;j<D;j++){
388  if(i != j) x_1.v[i] += A.m[i][j] * x.v[j];
389  }
390  x_1.v[i] = (b.v[i] - x_1.v[i]) / A.m[i][i];
391 
392  }
393  cpyV(x,x_1);
394  }
395 
396  return x;
397 }
398 
399 
411 
412  Vector x;
413 
414  double detA = det(A);
415 
416  // If the matrix has a determinant, we can get the inverse.
417  if(detA != 0){
418 
419  Matrix InvA;
420  inv(InvA,A);
421  multMV(x,InvA,b);
422  } else {
423  //printf("#Warning: singular matrix\n");
424  x = solve_jacobi(A,b);
425  }
426  return x;
427 }
428 
429 
434 
435 
436  if(hit_Rank != 0){
437  hit_tileDomainShapeAlloc(&graph,Node,global_shape);
438  }
439 
440  hit_tileDomainShapeAlloc(&e,Vector,global_shape);
441  hit_tileDomainShapeAlloc(&new_e,Vector,global_shape);
442 
444  hit_tileDomainShapeAlloc(&hessian,Matrix,incident_shape);
445 }
446 
447 
448 
453 
455  hit_tileFree(e);
459 
460 }
461 
462 
463 
464 
469 void print_help(char * name){
470  printf("%s [-n NEWTON METHOD ITERATIONS] [-j JACOBI METHOD ITERATIONS] FILE \n",name);
471 
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");
475 
476 }
477 
478 
479 
487 char *replace_str(char *str, const char *orig, const char *rep){
488 
489  static char buffer[4096];
490  char *p;
491 
492  if(!(p = strstr(str, orig))) return NULL;
493 
494  strncpy(buffer, str, (size_t) (p-str)); // Copy characters from 'str' start to 'orig' st$
495  buffer[p-str] = '\0';
496 
497  sprintf(buffer+(p-str), "%s%s", rep, p+strlen(orig));
498 
499  return buffer;
500 }
501 
502 
503 
504 void read_coordinates(char * name, HitShape shape_rc, HitTile_Node graph_rc){
505 
506  if(hit_Rank != 0) return;
507 
508 
509  FILE * file = fopen(name,"r");
510 
511 
512  if(file == NULL){
513  printf("File not exists\n");
514  return;
515  }
516 
517 
518 
519  char c = (char) getc(file);
520 
521  // Skip first comments
522  if(c == '%'){
523  while( c != '\n' ) c = (char) getc(file);
524  }
525  ungetc(c,file);
526 
527  int res,n,d,i,j;
528 
529  int n_shape = hit_sigCard(hit_shapeSig(shape_rc,0));
530  res = fscanf(file,"%d %d\n",&n,&d);
531  HIT_NOT_USED(res);
532 
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);
535 
536  srandom(0);
537 
538 
539  // Read the coordinates
540  for(j=0;j<D;j++){
541  for(i=0;i<n;i++){
542  if(j<d){
543  double number;
544  res = fscanf(file,"%lf\n",&number);
545  hit_tileElemAt(graph_rc,1,i).r.v[j] = number;
546  } else {
547  hit_tileElemAt(graph_rc,1,i).r.v[j] = 0.0;
548  }
549  }
550  }
551 
552 
553  for(i=0;i<n;i++){
554 
555  if ((random() % 100) < P_FIXED){
556  hit_tileElemAt(graph_rc,1,i).fixed = 1;
557  } else {
558  hit_tileElemAt(graph_rc,1,i).fixed = 0;
559  }
560 
561  }
562 
563  fclose(file);
564 
565 }
566 
567 
568 
572 void init_graph(int argc, char ** argv){
573 
574 
575  if(argc > 1){
576 
577  // Argument.
578  char c;
579 
580  // Parse the command-line arguments.
581  while ((c = (char) getopt (argc, argv, "hn:j:")) != -1)
582  switch (c) {
583  case 'n':
584  sscanf(optarg,"%d",&iter1);
585  break;
586  case 'j':
587  sscanf(optarg,"%d",&iter2);
588  break;
589  case 'h':
590  if(hit_Rank == 0) print_help(argv[0]);
591  hit_comFinalize();
592  exit(EXIT_SUCCESS);
593  default:
594  abort ();
595  break;
596  }
597 
598  if (hit_Rank != 0) return;
599 
600  // Graph file name.
601  char * graph_file = argv[optind];
602 
603  printfRoot("# Graph: %s\n",graph_file);
604  printfRoot("# Iterations %dx%d\n",iter1,iter2);
605 
606  Nvertices = hit_fileHBVertices(graph_file);
607 
608  printfRoot("# Vertices %d\n",Nvertices);
609 
610 
611  HitSig sig = hit_sig(0,Nvertices-1,1);
612  incident_shape = hit_shape(2,sig,sig);
613 
614 
615  global_shape = hit_shape(1,hit_sig(0,Nvertices-1,1));
616  hit_tileDomainShapeAlloc(&graph,Node,global_shape);
617 
618  // Coordinates file
619  char * coord_file = replace_str(graph_file,".rb","_coord.mtx");
620  printfRoot("# Coordinates: %s\n",coord_file);
621 
622  read_coordinates(coord_file,global_shape,graph);
623 
624 
625  int i;
626  int nfixed = 0;
627  for(i=0;i<Nvertices;i++){
628 
629  nfixed += hit_tileElemAt(graph,1,i).fixed;
630  }
631 
632  printf("# N fixed: %d\n",nfixed);
633 
634  } else {
635 
636  iter1 = iter2 = 3;
637  Nvertices = 8;
638  HitSig sig = hit_sig(0,Nvertices-1,1);
639  incident_shape = hit_shape(2,sig,sig);
640 
641  if(hit_Rank != 0) return;
642 
643 
644  // Allocate memory
645  global_shape = hit_shape(1,hit_sig(0,Nvertices-1,1));
646  hit_tileDomainShapeAlloc(&graph,Node,global_shape);
647 
648  // Set initial positions
649  Node n0 = {{{0,0,0}},1};
650  hit_tileElemAt(graph,1,0) = n0;
651  Node n1 = {{{100,0,0}},0};
652  hit_tileElemAt(graph,1,1) = n1;
653  Node n2 = {{{200,0,0}},0};
654  hit_tileElemAt(graph,1,2) = n2;
655  Node n3 = {{{300,0,0}},1};
656  hit_tileElemAt(graph,1,3) = n3;
657  Node n4 = {{{0,100,200}},1};
658  hit_tileElemAt(graph,1,4) = n4;
659  Node n5 = {{{100,100,0}},0};
660  hit_tileElemAt(graph,1,5) = n5;
661  Node n6 = {{{200,100,0}},0};
662  hit_tileElemAt(graph,1,6) = n6;
663  Node n7 = {{{300,100,0}},1};
664  hit_tileElemAt(graph,1,7) = n7;
665 
666  }
667 
668 
669 }
670 
671 
672 
677 double gradient_norm(){
678 
679  HitTile_double norm, sum_norm;
680 
681  hit_tileDomainAlloc(&norm, double, 1,1);
682  hit_tileDomainAlloc(&sum_norm, double, 1,1);
683 
684  hit_tileElemAt(norm,1,0) = 0.0;
685 
686  int vertex;
687 
688  hit_shapeIterator(vertex,shape,0){
689 
690  Node current = hit_tileElemAt(graph,1,vertex);
691 
692  if (!current.fixed){
693 
694  calculate_GH(vertex);
695 
696  Vector Gi = hit_tileElemAt(gradient,1,vertex);
697 
698  int d;
699  for(d=0;d<D;d++){
700  hit_tileElemAt(norm,1,0) += pow(Gi.v[d],2);
701  }
702 
703  //printf("vertex: %d: %.3f - %.3f - %.3f \n",vertex,pow(Gi.v[0],2),pow(Gi.v[1],2),pow(Gi.v[2],2));
704 
705  }
706  }
707 
708 
709  HitOp op_sum;
711 
712 
713  HitRanks root = {{0,0,0,-1}};
714  HitCom com_sum = hit_comReduce(lay, root, &norm, &sum_norm, HIT_DOUBLE, op_sum);
715  hit_comDo(&com_sum);
716  hit_comFree(com_sum);
717 
718  double res = 0;
719 
720  if( hit_Rank == 0 ){
721  res = sqrt(hit_tileElemAt(sum_norm,1,0));
722  }
723 
724  hit_tileFree(norm);
725  hit_tileFree(sum_norm);
726 
727 
728  return res;
729 
730 }
731 
732 
733 
734 
738 int main(int argc, char ** argv) {
739 
740  int i;
741  int vertex;
742  HitClock init_time, comp_time;
743 
744  hit_comInit(&argc,&argv);
745 
746  hit_clockStart(init_time);
747 
748 
749  if(hit_NProcs == 1){
750  printf("There is only one processor\n");
751  exit(EXIT_FAILURE);
752  }
753 
754  init_graph(argc,argv);
755 
756 
757 
758  // Create the topology object.
759  HitTopology topo = hit_topology(plug_topPlain);
760 
761  // Create the layout to send the incident matrix
762  HitLayout lay_all = hit_layout(plug_layBlocks,topo,hit_shape(1,hit_sig(0,hit_NProcs-1,1)));
763 
764  // Root processor.
765  HitRanks root = {{0,-1,-1,-1}};
766 
767  // Send the number of vertices to allocate the tile.
768  HitTile_int vertices;
769  hit_tileDomainAlloc(&vertices, int, 1,1);
770  hit_tileElemAt(vertices,1,0) = Nvertices;
771  HitCom sendNVertices = hit_comBroadcast(lay_all,root,&vertices,HIT_INT);
772  hit_comDo(&sendNVertices);
773  Nvertices = hit_tileElemAt(vertices,1,0);
774 
775  // Allocate the tile for the incident matrix
776  if(hit_Rank != 0){
777  HitSig sig = hit_sig(0,Nvertices-1,1);
778  incident_shape = hit_shape(2,sig,sig);
779  }
780 
781  // Create the global shape
782  global_shape = hit_shape(1,hit_sig(0,Nvertices-1,1));
783 
784  // Distribute the graph vertices among the processors.
785  lay = hit_layout(plug_layBlocks,topo,global_shape);
786  // Get shape
787  shape = hit_layShape(lay);
788 
789  // Allocate memory of all the variables
790  init_structures();
791 
792 
793 
794  // Create the HitVector type.
795  HitType HitVector;
796  hit_comTypeStruct(&HitVector,Vector,1, v,3,HIT_DOUBLE);
797 
798  // Create the HitNode type
799  HitType HitNode;
800  hit_comTypeStruct(&HitNode,Node,2, r,1,HitVector,fixed,1,HIT_INT);
801 
802 
803  // Create the communication objects.
804  HitCom comA = hit_comBroadcast(lay_all,root,&graph,HitNode);
805  HitPattern patB = hit_pattern( HIT_PAT_UNORDERED );
806  HitPattern patC = hit_pattern( HIT_PAT_UNORDERED );
807 
808  // @javfres 2015-10-09 Support for more processes than nodes
809  if(hit_layImActive( lay ))
810  for(i=0;i<hit_layNumActives(lay);i++){
811  //for(i=0;i<hit_NProcs;i++){
812 
813  HitRanks ranks = {{i,-1,-1,-1}};
814  HitShape otherShape = hit_layShapeOther(lay, ranks);
815  hit_patternAdd(&patB,
816  hit_comBroadcastSelect(lay,ranks,&graph,otherShape,HIT_COM_ARRAYCOORDS,HitNode));
817  hit_patternAdd(&patC,
818  hit_comBroadcastSelect(lay,ranks,&e,otherShape,HIT_COM_ARRAYCOORDS,HitVector));
819 
820  }
821 
822 
823  hit_clockStop(init_time);
824  hit_clockStart(comp_time);
825 
826  // Send all the data from the root proc to the other procs.
827  hit_comDo(&comA);
828 
829  // Calculate the inital gradient norm.
830  double ignorm = gradient_norm();
831  printfRoot("# Initial gradient norm: %lf\n",ignorm);
832 
833 
834  // Main loop for the newton method
835  for(i=0;i<iter1;i++){
836  printfRoot("# iter: %d/%d\n",i,iter1);
837 
838  // Iterate trough all the vertices and
839  // obtain the gradient and the hessian.
840 
841  hit_shapeIterator(vertex,shape,0){
842 
843  Node current = hit_tileElemAt(graph,1,vertex);
844  if (!current.fixed)
845  calculate_GH(vertex);
846  }
847 
848  // The initial position displacement is zero.
849  bzero(e.data,sizeof(Vector) * (size_t) Nvertices);
850 
851  // Loop for the jacobi method to solve the system.
852  int j;
853  for(j=0;j<iter2;j++){
854 
855  // Perform a iteration.
856  hit_shapeIterator(vertex,shape,0){
857  Node current = hit_tileElemAt(graph,1,vertex);
858  if (!current.fixed)
859  solve_system_iter(vertex);
860  }
861 
862 
863  // Update the displacement.
864  memcpy(e.data,new_e.data,sizeof(Vector)* (size_t) Nvertices);
865  hit_patternDo(patC);
866 
867  }
868 
869  // Update the position with the final displacement.
870  hit_shapeIterator(vertex,shape,0){
871 
872  Node * current = &(hit_tileElemAt(graph,1,vertex));
873  if (!current->fixed){
874  int d;
875  for(d=0;d<D;d++)
876  current->r.v[d] -= hit_tileElemAt(e,1,vertex).v[d];
877  }
878  }
879 
880  hit_patternDo(patB);
881  }
882 
883 
884  double gnorm = gradient_norm();
885  printfRoot("# Final gradient norm: %lf\n",gnorm);
886 
887  hit_clockStop(comp_time);
888  hit_clockWorldReduce(init_time);
889  hit_clockWorldReduce(comp_time);
890 
891  printfRoot("# Init time: %lf\n",hit_clockGetSeconds(init_time));
892  printfRoot("# Comp time: %lf\n",hit_clockGetSeconds(comp_time));
893  printfRoot("# Total time: %lf\n",hit_clockGetSeconds(init_time)+hit_clockGetSeconds(comp_time));
894 
895  hit_comFree(comA);
896  hit_comFree(sendNVertices);
897  hit_patternFree(&patB);
898  hit_patternFree(&patC);
899  free_structures();
900 
901  hit_layFree( lay );
902  hit_layFree( lay_all );
903  hit_topFree( topo );
904  hit_comFinalize();
905 
906  return 0;
907 }
908 
909 
913 Vector g(int i, int j){
914 
915  Vector gij;
916 
917  if(i == j){
918  clearV(gij);
919  return gij;
920  }
921 
922  // Get the two vectors
923  Vector r_i = (hit_tileElemAt(graph,1,i)).r;
924  Vector r_j = (hit_tileElemAt(graph,1,j)).r;
925 
926  // Auxiliar variables.
927  Vector v;
928 
929  // v = r_i - r_j
930  subV(v,r_i,r_j);
931  // n = ||v||
932  double n = norm(v);
933 
934  if( n == 0){
935  clearV(gij);
936  return gij;
937  }
938 
939  // gij = - K * ((L - n) / n) * (v)
940  double scalar = - K * ((L - n) / n);
941  multSV(gij,scalar,v);
942 
943  return gij;
944 }
945 
946 
947 
948 
952 Matrix W(int i, int j){
953 
954  Matrix Wij;
955 
956  if(i == j){
957  clearM(Wij);
958  return Wij;
959  }
960 
961  // Get the two vectors
962  Vector r_i = (hit_tileElemAt(graph,1,i)).r;
963  Vector r_j = (hit_tileElemAt(graph,1,j)).r;
964 
965  // Auxiliary variables.
966  Vector v;
967  Matrix a, b;
968  double scalar_a, scalar_b;
969 
970  // v = r_i - r_j;
971  subV(v,r_i,r_j);
972  // n = ||v||
973  double n = norm(v);
974 
975  if( n == 0){
976  clearM(Wij);
977  return Wij;
978  }
979 
980  // h = -((L-n)/n) * eye(D) + (L/(n^3)) * (v * v')
981  scalar_a = -((L-n)/(n));
982  multSI(a,scalar_a);
983  scalar_b = (L/(pow(n,3)));
984  multVVt(b,v,v);
985  multSM(b,scalar_b,b);
986  addM(Wij,a,b);
987  // h = K * h
988  multSM(Wij,K,Wij);
989 
990  return Wij;
991 }
992 
993 
998 void calculate_GH(int vertex){
999 
1000  Vector gs;
1001  Matrix Ws;
1002 
1003  clearV(gs);
1004  clearM(Ws);
1005 
1006  int edge;
1007  for(edge=0;edge<Nvertices;edge++){
1008 
1009  // Get the neighbor.
1010  int nghb_index = edge;
1011 
1012  // Calculate gradient[i]
1013  Vector gij = g(vertex,nghb_index);
1014  addV(gs,gs,gij);
1015 
1016  // Calculate part of the hessian[i][i]
1017  Matrix Wij = W(vertex,nghb_index);
1018  addM(Ws,Ws,Wij);
1019 
1020 
1021  // Calculate hessian[i][j] if j is variable.
1022  Node nj = hit_tileElemAt(graph,1,nghb_index);
1023 
1024  if(!nj.fixed){
1025  Matrix Hij = W(nghb_index,vertex);
1026  multSM(Hij,-1,Hij);
1027 
1028  hit_tileElemAt(hessian,2,vertex,nghb_index) = Hij;
1029 
1030 
1031  }
1032 
1033  }
1034 
1035  hit_tileElemAt(gradient,1,vertex) = gs;
1036  hit_tileElemAt(hessian,2,vertex,vertex) = Ws;
1037 }
1038 
1039 
1040 
1046 void solve_system_iter(int vertex){
1047 
1048  Matrix Hii = hit_tileElemAt(hessian,2,vertex,vertex);
1049  Vector Gi = hit_tileElemAt(gradient,1,vertex);
1050 
1051  Vector new_ei;
1052 
1053  // Aux sum vector.
1054  Vector sum;
1055  clearV(sum);
1056 
1057  int edge;
1058 
1059  for(edge=0;edge<Nvertices;edge++){
1060 
1061  // Get the neighbor.
1062  int nghb_index = edge;
1063  if(nghb_index == vertex) continue;
1064 
1065  // If the node is fixed, it is no part of the ecuation system.
1066  Node nj = hit_tileElemAt(graph,1,nghb_index);
1067  if(nj.fixed) continue;
1068 
1069  // Get the j displacement and ij hessian part.
1070  Vector ej = hit_tileElemAt(e,1,nghb_index);
1071  Matrix Hij = hit_tileElemAt(hessian,2,vertex,nghb_index);
1072 
1073  // Aux par vector.
1074  Vector part;
1075 
1076  // Calculate and add the contribution of vertex j.
1077  // SUM_{i not j} ( Hij e_j )
1078  multMV(part,Hij,ej);
1079  addV(sum,sum,part);
1080 
1081  }
1082 
1083 
1084  // Gi - SUM_{i not j} ( Hij e_j )
1085  subV(sum,Gi,sum);
1086 
1087  // e_i = H_ii^(-1) * (Gi - SUM_{i not j} ( Hij e_j ) )
1088  new_ei = solve(Hii,sum);
1089 
1090  hit_tileElemAt(new_e,1,vertex) = new_ei;
1091 
1092 
1093 }
Vector g(int i, int j)
#define hit_shape(nd,...)
Definition: hit_sshape.h:175
MPI_Op HitOp
Definition: hit_com.h:118
void print_help(char *name)
Definition: heat.c:234
#define hit_layShape(lay)
Definition: hit_layout.h:650
#define hit_tileNewType(baseType)
Definition: hit_tile.h:163
void hit_comOpSumDouble(void *, void *, int *, HitType *)
Definition: hit_com.c:2665
#define A
Definition: mg.c:64
HitOp op_sum
Definition: mg.c:176
double v[D]
Definition: spring_bitmap.c:86
#define ITER2
#define HIT_INT
Definition: hit_com.h:74
int Nvertices
Definition: spring_dense.c:121
HitCom hit_comBroadcastSelect(HitLayout lay, HitRanks root, const void *tile, HitShape selection, int modeSelect, HitType baseType)
Definition: hit_com.c:731
#define hit_Rank
Definition: hit_com.h:140
int iter2
Definition: spring_bitmap.c:80
int iter
Definition: mmult_bit.c:68
#define hit_tileElemAt(var, ndims,...)
Definition: hit_tile.h:519
#define addV(r, a, b)
#define HIT_LAYOUT_NULL_STATIC
Definition: hit_layout.h:300
void hit_comFree(HitCom issue)
Definition: hit_com.c:1995
#define K
char * replace_str(char *str, const char *orig, const char *rep)
Matrix W(int i, int j)
void hit_comDo(HitCom *issue)
Definition: hit_com.c:2408
#define multSI(r, s)
HitTile_Matrix hessian
#define subV(r, a, b)
#define hit_tileDomainAlloc(newVarP, baseType, numDims,...)
Definition: hit_tile.h:336
void solve_system_iter(int vertex)
#define hit_topology(name,...)
Definition: hit_topology.h:308
#define hit_layShapeOther(lay, ranks)
Definition: hit_layout.h:732
#define hit_sigCard(sig)
Definition: hit_sig.h:162
Definition: hit_sig.h:79
#define hit_fileHBVertices(hbfile)
Definition: hit_file.h:157
HitTile_Vector graph
#define ITER_SJ
#define r(a, b, c)
#define HIT_COM_ARRAYCOORDS
Definition: hit_com.h:345
#define clearM(mm)
#define HIT_PAT_UNORDERED
Definition: hit_pattern.h:81
void hit_patternAdd(HitPattern *pattern, HitCom comm)
Definition: hit_pattern.c:61
#define multSV(r, s, vv)
#define D
Vector solve(Matrix A, Vector b)
void free_structures()
#define hit_NProcs
Definition: hit_com.h:142
#define P_FIXED
#define hit_tileDomainShapeAlloc(newVarP, baseType, shape)
Definition: hit_tile.h:354
HitTile_int fixed
#define ITER1
#define x
int fixed
Definition: spring_dense.c:101
#define printfRoot(...)
#define det(mm)
void hit_topFree(HitTopology topo)
Definition: hit_topology.c:129
#define hit_comOp(function, operation)
Definition: hit_com.h:1036
MPI_Datatype HitType
Definition: hit_com.h:70
#define multSM(r, s, mm)
#define HIT_NOT_USED(x)
Definition: hit_error.h:50
#define addM(r, a, b)
Vector solve_jacobi(Matrix A, Vector b)
#define hit_clockWorldReduce(c)
Definition: hit_utils.h:156
void hit_patternFree(HitPattern *pattern)
Definition: hit_pattern.c:94
#define inv(r, mm)
void hit_comInit(int *pargc, char **pargv[])
Definition: hit_com.c:111
HitShape shape
#define hit_layImActive(lay)
Definition: hit_layout.h:797
HitShape incident_shape
Definition: spring_dense.c:128
#define L
#define norm(vv)
#define cpyV(r, a)
void hit_layFree(HitLayout lay)
Definition: hit_layout.c:2251
#define hit_clockStop(c)
Definition: hit_utils.h:109
HitShape global_shape
int main(int argc, char *argv[])
Definition: cannonAsync.c:62
#define m(a, b, c)
HitTile_Vector new_e
#define multMV(r, mm, vv)
Vector r
Definition: spring_dense.c:100
#define hit_comReduce(lay, root, tilePSend, tilePRecv, baseType, operation)
Definition: hit_com.h:725
#define hit_comTypeStruct(new_type, Nstruct, n,...)
Definition: hit_com.h:208
void init_graph(HitTile_double graph, HitShape shape_global)
Definition: heat.c:244
HitTile_Vector gradient
void init_structures()
HitTile_Vector e
#define hit_shapeIterator(var, shape, dim)
Definition: hit_sshape.h:603
int hit_layNumActives(HitLayout lay)
Definition: hit_layout.c:2213
void calculate_GH(int vertex)
#define v(a, b, c)
#define hit_comBroadcast(lay, root, tile, baseType)
Definition: hit_com.h:673
void read_coordinates(char *name, HitShape shape_rc, HitTile_Node graph_rc)
#define hit_layout(name, topo,...)
Definition: hit_layout.h:415
double gradient_norm()
HitLayout lay
Definition: heat.c:107
#define hit_shapeSig(shape, dim)
Definition: hit_sshape.h:400
#define hit_min(a, b)
Definition: hit_funcop.h:65
#define hit_tileFree(var)
Definition: hit_tile.h:369
#define clearV(vv)
#define multVVt(r, a, b)
#define HIT_DOUBLE
Definition: hit_com.h:78
#define hit_patternDo(pattern)
Definition: hit_pattern.h:157
void hit_comFinalize()
Definition: hit_com.c:159
#define hit_clockGetSeconds(c)
Definition: hit_utils.h:190
int iter1
Definition: spring_bitmap.c:78
double ** buffer
Definition: refMPIluBack.c:103
#define hit_clockStart(c)
Definition: hit_utils.h:87
double m[D][D]
Definition: spring_bitmap.c:93