Hitmap 1.3
 All Data Structures Namespaces Files Functions Variables Typedefs Friends Macros Groups Pages
spring_bitmap.c
Go to the documentation of this file.
1 
18 /*
19  * <license>
20  *
21  * Hitmap v1.2
22  *
23  * This software is provided to enhance knowledge and encourage progress in the scientific
24  * community. It should be used only for research and educational purposes. Any reproduction
25  * or use for commercial purpose, public redistribution, in source or binary forms, with or
26  * without modifications, is NOT ALLOWED without the previous authorization of the copyright
27  * holder. The origin of this software must not be misrepresented; you must not claim that you
28  * wrote the original software. If you use this software for any purpose (e.g. publication),
29  * a reference to the software package and the authors must be included.
30  *
31  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS "AS IS" AND ANY
32  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
33  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
34  * THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
35  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
36  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
37  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
38  * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
39  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
40  *
41  * Copyright (c) 2007-2015, Trasgo Group, Universidad de Valladolid.
42  * All rights reserved.
43  *
44  * More information on http://trasgo.infor.uva.es/
45  *
46  * </license>
47 */
48 
49 #include <stdio.h>
50 #include <stdlib.h>
51 #include <hitmap.h>
52 #include <hit_com.h>
53 #include <unistd.h>
54 #include <string.h>
55 
56 
57 /* OUTPUT FUNCIONS */
58 #define printfRootInternal(...) { if( hit_Rank == 0 ) { printf(__VA_ARGS__); fflush(stdout); }}
59 #define printfRoot(...) printfRootInternal(__VA_ARGS__)
60 
62 #define D (3)
63 
65 #define K (1)
66 
67 #define L (0.001)
68 
70 #define P_FIXED 10
71 
73 #define ITER1 10
74 
75 #define ITER2 100
76 
78 int iter1 = ITER1;
80 int iter2 = ITER2;
81 
85 typedef struct {
86  double v[D];
87 } Vector;
88 
92 typedef struct {
93  double m[D][D];
94 } Matrix;
95 
96 
97 // Define the tile_Vector type.
99 
100 // Define the tile_Matrix type.
102 
103 // Define the tile_int type.
104 hit_tileNewType(int);
105 
106 // Define the tile_double type.
107 hit_tileNewType(double);
108 
109 
111 HitTile_Vector gradient;
112 
114 HitTile_Matrix hessian;
115 
117 HitTile_Vector e;
118 
120 HitTile_Vector new_e;
121 
124 
127 
130 
132 HitTile_Vector graph;
133 
135 HitTile_Vector global_graph;
136 
138 HitTile_int fixed;
139 
141 HitTile_int global_fixed;
142 
143 
146 
147 
152 void calculate_GH(int vertex);
153 
158 void solve_system_iter(int vertex);
159 
160 
161 #ifdef DEBUG
162 
166 void print_matrix(Matrix m){
167 
168  int i,j;
169  for(i=0;i<D;i++){
170  for(j=0;j<D;j++){
171  printf("%.2lf\t",m.m[i][j]);
172  }
173  printf("\n");
174  }
175 }
176 
181 void print_vector(Vector v){
182 
183  int i;
184  for(i=0;i<D;i++){
185  printf("%.2lf\t",v.v[i]);
186  }
187  printf("\n");
188 }
189 #endif
190 
191 
195 #define clearV(vv) clearVInternal((vv).v)
196 #define clearVInternal(v) \
197  ((v)[0] = (v)[1] = (v)[2] = 0.0 )
198 
202 #define cpyV(r,a) cpyVInternal((r).v,(a).v)
203 #define cpyVInternal(r,a) \
204  {(r)[0] = (a)[0]; (r)[1] = (a)[1]; (r)[2] = (a)[2]; }
205 
209 #define addV(r,a,b) addVInternal((r).v,(a).v,(b).v)
210 #define addVInternal(r,a,b) \
211  {(r)[0] = (a)[0] + (b)[0]; (r)[1] = (a)[1] + (b)[1]; (r)[2] = (a)[2] + (b)[2]; }
212 
213 
217 #define subV(r,a,b) subVInternal((r).v,(a).v,(b).v)
218 #define subVInternal(r,a,b) \
219  {(r)[0] = (a)[0] - (b)[0]; (r)[1] = (a)[1] - (b)[1]; (r)[2] = (a)[2] - (b)[2]; }
220 
224 #define norm(vv) normInternal((vv).v)
225 #define normInternal(v) \
226  (sqrt( pow((v)[0],2) + pow((v)[1],2) + pow((v)[2],2) ))
227 
231 #define multSV(r,s,vv) multSVInternal((r).v,s,(vv).v)
232 #define multSVInternal(r,s,v) \
233  {(r)[0] = (s) * (v)[0]; (r)[1] = (s) * (v)[1]; (r)[2] = (s) * (v)[2]; }
234 
238 #define clearM(mm) clearMInternal((mm).m)
239 #define clearMInternal(m) \
240  ((m)[0][0] = (m)[0][1] = (m)[0][2] = \
241  (m)[1][0] = (m)[1][1] = (m)[1][2] = \
242  (m)[2][0] = (m)[2][1] = (m)[2][2] = 0.0 )
243 
244 
254 #define det(mm) detInternal((mm).m)
255 #define detInternal(m) \
256  ( (m)[0][0] * (m)[1][1] * (m)[2][2] \
257  + (m)[0][1] * (m)[1][2] * (m)[2][0] \
258  + (m)[0][2] * (m)[1][0] * (m)[2][1] \
259  - (m)[0][0] * (m)[1][2] * (m)[2][1] \
260  - (m)[0][1] * (m)[1][0] * (m)[2][2] \
261  - (m)[0][2] * (m)[1][1] * (m)[2][0] )
262 
266 #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)
267 #define init_matrixInternal(m,a,b,c,d,e,f,g,h,i) \
268  {(m)[0][0] = a; (m)[0][1] = b; (m)[0][2] = c; \
269  (m)[1][0] = d; (m)[1][1] = e; (m)[1][2] = f; \
270  (m)[2][0] = g; (m)[2][1] = h; (m)[2][2] = i; }
271 
275 #define multSM(r,s,mm) multSMInternal((r).m,s,(mm).m)
276 #define multSMInternal(r,s,m) \
277  {(r)[0][0] = (s) * (m)[0][0]; (r)[0][1] = (s) * (m)[0][1]; (r)[0][2] = (s) * (m)[0][2]; \
278  (r)[1][0] = (s) * (m)[1][0]; (r)[1][1] = (s) * (m)[1][1]; (r)[1][2] = (s) * (m)[1][2]; \
279  (r)[2][0] = (s) * (m)[2][0]; (r)[2][1] = (s) * (m)[2][1]; (r)[2][2] = (s) * (m)[2][2]; }
280 
284 #define multSI(r,s) multSIInternal((r).m,s)
285 #define multSIInternal(r,s) \
286  {(r)[0][0] = (s); (r)[0][1] = 0; (r)[0][2] = 0; \
287  (r)[1][0] = 0; (r)[1][1] = (s); (r)[1][2] = 0; \
288  (r)[2][0] = 0; (r)[2][1] = 0; (r)[2][2] = (s); }
289 
290 
294 #define addM(r,a,b) addMInternal((r).m,(a).m,(b).m)
295 #define addMInternal(r,a,b) \
296  {(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]; \
297  (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]; \
298  (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]; }
299 
304 #define multVVt(r,a,b) multVVtInternal((r).m,(a).v,(b).v)
305 #define multVVtInternal(r,a,b) \
306  {(r)[0][0] = (a)[0] * (b)[0]; (r)[0][1] = (a)[0] * (b)[1]; (r)[0][2] = (a)[0] * (b)[2]; \
307  (r)[1][0] = (a)[1] * (b)[0]; (r)[1][1] = (a)[1] * (b)[1]; (r)[1][2] = (a)[1] * (b)[2]; \
308  (r)[2][0] = (a)[2] * (b)[0]; (r)[2][1] = (a)[2] * (b)[1]; (r)[2][2] = (a)[2] * (b)[2]; }
309 
310 
314 #define multMV(r,mm,vv) multMVInternal((r).v,(mm).m,(vv).v)
315 #define multMVInternal(r,m,v) \
316  {(r)[0] = (m)[0][0] * (v)[0] + (m)[0][1] * (v)[1] + (m)[0][2] * (v)[2]; \
317  (r)[1] = (m)[1][0] * (v)[0] + (m)[1][1] * (v)[1] + (m)[1][2] * (v)[2]; \
318  (r)[2] = (m)[2][0] * (v)[0] + (m)[2][1] * (v)[1] + (m)[2][2] * (v)[2]; }
319 
320 
324 #define inv(r,mm) invInternal((r).m,(mm).m)
325 #define invInternal(r,m) \
326  {(r)[0][0] = (m)[1][1] * (m)[2][2] - (m)[1][2] * (m)[2][1]; \
327  (r)[0][1] = - (m)[0][1] * (m)[2][2] + (m)[0][2] * (m)[2][1]; \
328  (r)[0][2] = (m)[0][1] * (m)[1][2] - (m)[0][2] * (m)[1][1]; \
329  (r)[1][0] = - (m)[1][0] * (m)[2][2] + (m)[1][2] * (m)[2][0]; \
330  (r)[1][1] = (m)[0][0] * (m)[2][2] - (m)[0][2] * (m)[2][0]; \
331  (r)[1][2] = - (m)[0][0] * (m)[1][2] + (m)[0][2] * (m)[1][0]; \
332  (r)[2][0] = (m)[1][0] * (m)[2][1] - (m)[1][1] * (m)[2][0]; \
333  (r)[2][1] = - (m)[0][0] * (m)[2][1] + (m)[0][1] * (m)[2][0]; \
334  (r)[2][2] = (m)[0][0] * (m)[1][1] - (m)[0][1] * (m)[1][0]; \
335  multSMInternal(r,1/detInternal(m),r); }
336 
340 #define multMM(r,a,b) multMMInternal((r).m,(a).m,(b).m)
341 #define multMMInternal(r,a,b) \
342  {(r)[0][0] = (a)[0][0] * (b)[0][0] + (a)[0][1] * (b)[1][0] + (a)[0][2] * (b)[2][0]; \
343  (r)[0][1] = (a)[0][0] * (b)[0][1] + (a)[0][1] * (b)[1][1] + (a)[0][2] * (b)[2][1]; \
344  (r)[0][2] = (a)[0][0] * (b)[0][2] + (a)[0][1] * (b)[1][2] + (a)[0][2] * (b)[2][2]; \
345  (r)[1][0] = (a)[1][0] * (b)[0][0] + (a)[1][1] * (b)[1][0] + (a)[1][2] * (b)[2][0]; \
346  (r)[1][1] = (a)[1][0] * (b)[0][1] + (a)[1][1] * (b)[1][1] + (a)[1][2] * (b)[2][1]; \
347  (r)[1][2] = (a)[1][0] * (b)[0][2] + (a)[1][1] * (b)[1][2] + (a)[1][2] * (b)[2][2]; \
348  (r)[2][0] = (a)[2][0] * (b)[0][0] + (a)[2][1] * (b)[1][0] + (a)[2][2] * (b)[2][0]; \
349  (r)[2][1] = (a)[2][0] * (b)[0][1] + (a)[2][1] * (b)[1][1] + (a)[2][2] * (b)[2][1]; \
350  (r)[2][2] = (a)[2][0] * (b)[0][2] + (a)[2][1] * (b)[1][2] + (a)[2][2] * (b)[2][2]; }
351 
352 
361 
362 // Number of iteration of the jacobi iterative method.
363 #define ITER_SJ 10
364 
365  Vector x;
366 
367  // The new approximation.
368  Vector x_1;
369 
370  // Clear the vectors.
371  clearV(x);
372  clearV(x_1);
373 
374  int iter;
375  int i,j;
376 
377  // Iterative method.
378  for(iter=0; iter<ITER_SJ; iter++){
379  for(i=0;i<D;i++){
380 
381  if( A.m[i][i] == 0 ) continue;
382 
383  // calculate the new value.
384  x_1.v[i] = 0.0;
385  for(j=0;j<D;j++){
386  if(i != j) x_1.v[i] += A.m[i][j] * x.v[j];
387  }
388  x_1.v[i] = (b.v[i] - x_1.v[i]) / A.m[i][i];
389 
390  }
391  cpyV(x,x_1);
392  }
393 
394  return x;
395 }
396 
397 
409 
410  Vector x;
411 
412  double detA = det(A);
413 
414  // If the matrix has a determinant, we can get the inverse.
415  if(detA != 0){
416 
417  Matrix InvA;
418  inv(InvA,A);
419  multMV(x,InvA,b);
420  } else {
421  //printf("#Warning: singular matrix\n");
422  x = solve_jacobi(A,b);
423  }
424  return x;
425 }
426 
431 
436 
438  hit_gbTileDomainShapeAlloc(&fixed, int, ext_shape, HIT_VERTICES);
439 
440 }
441 
446 
447  if(hit_Rank == 0){
450  }
451 
454 
455  hit_tileFree(e);
459 
460 }
461 
462 
467 void print_help(char * name){
468  printf("%s [-n NEWTON METHOD ITERATIONS] [-j JACOBI METHOD ITERATIONS] FILE \n",name);
469 
470  printf(" -n NEWTON METHOD ITERATIONS number of iterations (default 10)\n");
471  printf(" -j JACOBI METHOD ITERATIONS number of iterations (default 100)\n");
472  printf(" FILE input graph file\n");
473 
474 }
475 
483 char *replace_str(char *str, const char *orig, const char *rep){
484 
485  static char buffer[4096];
486  char *p;
487 
488  if(!(p = strstr(str, orig))) return NULL;
489 
490  strncpy(buffer, str, (size_t) (p-str)); // Copy characters from 'str' start to 'orig' st$
491  buffer[p-str] = '\0';
492 
493  sprintf(buffer+(p-str), "%s%s", rep, p+strlen(orig));
494 
495  return buffer;
496 }
497 
498 
499 
505 
506  srandom(0);
507  if(hit_Rank != 0) return;
508 
509  int n = hit_bShapeNvertices(global_shape);
510 
511  int j,i;
512 
513  // Set the coordinates
514  for(j=0;j<D;j++){
515  for(i=0;i<n;i++){
516  double number = (double) random();
517  hit_gbTileVertexAt(global_graph,i).v[j] = number;
518  }
519  }
520 
521  srandom(0);
522  for(i=0;i<n;i++){
523 
524  if ((random() % 100) < P_FIXED){
526  } else {
528  }
529  }
530 
531 }
532 
533 
534 
535 
536 
537 
541 void init_graph(int argc, char ** argv){
542 
543  // Define the spring structure.
544  global_shape = HIT_BITMAP_SHAPE_NULL;
545 
546  if(argc > 1){
547 
548  // Argument.
549  char c;
550 
551  // Parse the command-line arguments.
552  while ((c = (char) getopt (argc, argv, "hn:j:")) != -1)
553  switch (c) {
554  case 'n':
555  sscanf(optarg,"%d",&iter1);
556  break;
557  case 'j':
558  sscanf(optarg,"%d",&iter2);
559  break;
560  case 'h':
561  if(hit_Rank == 0) print_help(argv[0]);
562  hit_comFinalize();
563  exit(EXIT_SUCCESS);
564  default:
565  abort ();
566  break;
567  }
568 
569  if (hit_Rank != 0) return;
570 
571  // Graph file name.
572  char * graph_file = argv[optind];
573 
574  printfRoot("# Graph: %s\n",graph_file);
575  printfRoot("# Iterations %d x %d\n",iter1,iter2);
576 
577  global_shape = hit_fileHBReadBitmap(graph_file);
578 
579  printfRoot("# Vertices %d\n",hit_bShapeNvertices(global_shape));
580 
581 
582  // Allocate memory
585 
586  // Random coordinates
588 
589  int i;
590  int nfixed = 0;
591  for(i=0;i<hit_bShapeNvertices(global_shape);i++){
592  nfixed += hit_gbTileVertexAt(global_fixed,i);
593  }
594 
595  printf("# N fixed: %d\n",nfixed);
596 
597  } else {
598 
599  iter1 = iter2 = 3;
600 
601  if(hit_Rank != 0) return;
602 
603  hit_bShapeAddEdge2(&global_shape,0,1);
604  hit_bShapeAddEdge2(&global_shape,1,2);
605  hit_bShapeAddEdge2(&global_shape,2,3);
606  hit_bShapeAddEdge2(&global_shape,4,5);
607  hit_bShapeAddEdge2(&global_shape,5,6);
608  hit_bShapeAddEdge2(&global_shape,6,7);
609  hit_bShapeAddEdge2(&global_shape,0,4);
610  hit_bShapeAddEdge2(&global_shape,1,5);
611  hit_bShapeAddEdge2(&global_shape,2,6);
612  hit_bShapeAddEdge2(&global_shape,3,7);
613  hit_bShapeAddEdge2(&global_shape,1,4);
614  hit_bShapeAddEdge2(&global_shape,2,5);
615 
616  // Include each (i,i) edge.
617  hit_bShapeAddEdge2(&global_shape,0,0);
618  hit_bShapeAddEdge2(&global_shape,1,1);
619  hit_bShapeAddEdge2(&global_shape,2,2);
620  hit_bShapeAddEdge2(&global_shape,3,3);
621  hit_bShapeAddEdge2(&global_shape,4,4);
622  hit_bShapeAddEdge2(&global_shape,5,5);
623  hit_bShapeAddEdge2(&global_shape,6,6);
624  hit_bShapeAddEdge2(&global_shape,7,7);
625 
626  // Allocate memory
629 
630  // Set initial positions
631  Vector v0 = {{0,0,0}};
633  Vector v1 = {{100,0,0}};
635  Vector v2 = {{200,0,0}};
637  Vector v3 = {{300,0,0}};
639  Vector v4 = {{0,100,200}};
641  Vector v5 = {{100,100,0}};
643  Vector v6 = {{200,100,0}};
645  Vector v7 = {{300,100,0}};
647 
648  // Set the fix/free status
657 
658  }
659 }
660 
661 
662 
667 double gradient_norm(){
668 
669  HitTile_double norm, sum_norm;
670 
671  hit_tileDomainAlloc(&norm, double, 1,1);
672  hit_tileDomainAlloc(&sum_norm, double, 1,1);
673 
674  hit_tileElemAt(norm,1,0) = 0.0;
675 
676  int vertex;
677 
678  hit_bShapeVertexIterator(vertex,shape) {
679 
680  if (!hit_gbTileVertexAt(fixed,vertex)){
681 
682  calculate_GH(vertex);
683 
684  Vector Gi = hit_gbTileVertexAt(gradient,vertex);
685 
686  int d;
687  for(d=0;d<D;d++){
688  hit_tileElemAt(norm,1,0) += pow(Gi.v[d],2);
689  }
690  }
691  }
692 
693  HitOp op_sum;
695 
696 
697  HitRanks root = {{0,0,0,-1}};
698  HitCom com_sum = hit_comReduce(lay, root, &norm, &sum_norm, HIT_DOUBLE, op_sum);
699  hit_comDo(&com_sum);
700  hit_comFree(com_sum);
701 
702  double res = 0;
703 
704  if( hit_Rank == 0 ){
705  res = sqrt(hit_tileElemAt(sum_norm,1,0));
706  }
707 
708  hit_tileFree(norm);
709  hit_tileFree(sum_norm);
710 
711 
712  return res;
713 
714 }
715 
716 
717 
718 
722 int main(int argc, char ** argv) {
723 
724  int i;
725  int vertex;
726  HitClock init_time, comp_time;
727 
728  hit_comInit(&argc,&argv);
729 
730  hit_clockStart(init_time);
731  init_graph(argc,argv);
732 
733  // Create the topology object.
734  HitTopology topo = hit_topology(plug_topPlain);
735 
736  // Distribute the graph among the processors.
737  lay = hit_layout(plug_layBitmap,topo,&global_shape);
738 
739 
740  // Get shape and extended shape
741  shape = hit_layShape(lay);
742  ext_shape = hit_bShapeExpand(shape,global_shape,1);
743 
744  // Allocate memory of all the variables
745  init_structures();
746 
747  // Create the HitVector type.
748  HitType HitVector;
749  hit_comTypeStruct(&HitVector,Vector,1, v,3,HIT_DOUBLE);
750 
751  // Create the communication objects.
752  HitCom comA = hit_comSparseScatter(lay, &global_graph, &graph, HitVector);
754  HitCom comA3 = hit_comSparseUpdate(lay, &fixed, HIT_INT);
755  HitCom comB = hit_comSparseUpdate(lay, &graph, HitVector);
756  HitCom comC = hit_comSparseUpdate(lay, &e, HitVector);
757 
758  hit_clockStop(init_time);
759  hit_clockStart(comp_time);
760 
761  // Send all the data from the root proc to the other procs.
762  hit_comDo(&comA); hit_comDo(&comA2); hit_comDo(&comA3);
763  hit_comDo(&comB);
764 
765  // Calculate the inital gradient norm.
766  // double gnorm1 = gradient_norm();
767  // printfRoot("# Initial gradient norm: %lf\n",gnorm1);
768 
769  // Main loop for the newton method
770  for(i=0;i<iter1;i++){
771  printfRoot("# iter: %d/%d\n",i,iter1);
772 
773  // Iterate trough all the vertices and
774  // obtain the gradient and the hessian.
775  hit_bShapeVertexIterator(vertex,shape) {
776  if (!hit_gbTileVertexAt(fixed,vertex))
777  calculate_GH(vertex);
778  }
779 
780  // The initial position displacement is zero.
782 
783 
784  // Loop for the jacobi method to solve the system.
785  int j;
786  for(j=0;j<iter2;j++){
787 
788  // Perform a iteration.
789  hit_bShapeVertexIterator(vertex,shape) {
790  if (!hit_gbTileVertexAt(fixed,vertex))
791  solve_system_iter(vertex);
792  }
793 
794 
795  // Update the displacement.
797  hit_comDo(&comC);
798 
799  }
800 
801  // Update the position with the final displacement.
802  hit_bShapeVertexIterator(vertex,shape) {
803  if (!hit_gbTileVertexAt(fixed,vertex)){
804  Vector * current = &(hit_gbTileVertexAt(graph,vertex));
805  subV(*current,*current,hit_gbTileVertexAt(e,vertex));
806  }
807  }
808 
809  hit_comDo(&comB);
810 
811  }
812 
813 
814  hit_clockStop(comp_time);
815  hit_clockWorldReduce(init_time);
816  hit_clockWorldReduce(comp_time);
817 
818  double gnorm = gradient_norm();
819  printfRoot("# Final gradient norm: %lf\n",gnorm);
820  printfRoot("# Init time: %lf\n",hit_clockGetSeconds(init_time));
821  printfRoot("# Comp time: %lf\n",hit_clockGetSeconds(comp_time));
822  printfRoot("# Total time: %lf\n",hit_clockGetSeconds(init_time)+hit_clockGetSeconds(comp_time));
823 
824 
825  hit_comFree(comA);
826  hit_comFree(comB);
827  hit_comFree(comC);
828 
829  free_structures();
830  hit_layFree( lay );
831  hit_topFree( topo );
832  hit_comFinalize();
833 
834  return 0;
835 }
836 
837 
838 
839 
843 Vector g(int i, int j){
844 
845  Vector gij;
846 
847  if(i == j){
848  clearV(gij);
849  return gij;
850  }
851 
852  // Get the two vectors
853  Vector r_i = (hit_gbTileVertexAt(graph,i));
854  Vector r_j = (hit_gbTileVertexAt(graph,j));
855 
856  // Auxiliar variables.
857  Vector v;
858 
859  // v = r_i - r_j
860  subV(v,r_i,r_j);
861  // n = ||v||
862  double n = norm(v);
863 
864  if( n == 0){
865  clearV(gij);
866  return gij;
867  }
868 
869  // gij = - K * ((L - n) / n) * (v)
870  double scalar = - K * ((L - n) / n);
871  multSV(gij,scalar,v);
872 
873  return gij;
874 }
875 
876 
880 Matrix W(int i, int j){
881 
882  Matrix Wij;
883 
884  if(i == j){
885  clearM(Wij);
886  return Wij;
887  }
888 
889  // Get the two vectors
890  Vector r_i = (hit_gbTileVertexAt(graph,i));
891  Vector r_j = (hit_gbTileVertexAt(graph,j));
892 
893  // Auxiliary variables.
894  Vector v;
895  Matrix a, b;
896  double scalar_a, scalar_b;
897 
898  // v = r_i - r_j;
899  subV(v,r_i,r_j);
900  // n = ||v||
901  double n = norm(v);
902 
903  if( n == 0){
904  clearM(Wij);
905  return Wij;
906  }
907 
908  // h = -((L-n)/n) * eye(D) + (L/(n^3)) * (v * v')
909  scalar_a = -((L-n)/(n));
910  multSI(a,scalar_a);
911  scalar_b = (L/(pow(n,3)));
912  multVVt(b,v,v);
913  multSM(b,scalar_b,b);
914  addM(Wij,a,b);
915  // h = K * h
916  multSM(Wij,K,Wij);
917 
918  return Wij;
919 }
920 
921 
927 void calculate_GH(int vertex){
928 
929  Vector gs;
930  Matrix Ws;
931 
932  clearV(gs);
933  clearM(Ws);
934 
935  int edge;
936  hit_bShapeEdgeIterator(edge,ext_shape,vertex){
937 
938  // Get the neighbor.
939  int nghb_index = hit_bShapeEdgeTarget(ext_shape,edge);
940 
941  // Calculate gradient[i]
942  Vector gij = g(vertex,nghb_index);
943  addV(gs,gs,gij);
944 
945  // Calculate part of the hessian[i][i]
946  Matrix Wij = W(vertex,nghb_index);
947  addM(Ws,Ws,Wij);
948 
949  // Calculate hessian[i][j] if j is variable.
950  if(!hit_gbTileVertexAt(fixed,nghb_index)){
951 
952  Matrix Hij = W(nghb_index,vertex);
953  multSM(Hij,-1,Hij);
954 
955  hit_gbTileEdgeIteratorAt(hessian,vertex,edge) = Hij;
956  }
957  }
958 
959  hit_gbTileVertexAt(gradient,vertex) = gs;
960  hit_gbTileEdgeAt(hessian,vertex,vertex) = Ws;
961 
962 }
963 
969 void solve_system_iter(int vertex){
970 
971  Matrix Hii = hit_gbTileEdgeAt(hessian,vertex,vertex);
972  Vector Gi = hit_gbTileVertexAt(gradient,vertex);
973 
974  // Aux sum vector.
975  Vector sum;
976  clearV(sum);
977 
978  int edge;
979  hit_bShapeEdgeIterator(edge,ext_shape,vertex){
980 
981  // Get the neighbor.
982  int nghb_index = hit_bShapeEdgeTarget(ext_shape,edge);
983  if(nghb_index == vertex) continue;
984 
985  // If the node is fixed, it is no part of the equation system.
986  if(hit_gbTileVertexAt(fixed,nghb_index)) continue;
987 
988  // Get the j displacement and ij hessian part.
989  Vector ej = hit_gbTileVertexAt(e,nghb_index);
990 
991  Matrix Hij = hit_gbTileEdgeIteratorAt(hessian,vertex,edge);
992 
993  // Aux par vector.
994  Vector part;
995 
996  // Calculate and add the contribution of vertex j.
997  // SUM_{i not j} ( Hij e_j )
998  multMV(part,Hij,ej);
999  addV(sum,sum,part);
1000 
1001  }
1002 
1003  // Gi - SUM_{i not j} ( Hij e_j )
1004  subV(sum,Gi,sum);
1005 
1006  // e_i = H_ii^(-1) * (Gi - SUM_{i not j} ( Hij e_j ) )
1007  hit_gbTileVertexAt(new_e,vertex) = solve(Hii,sum);
1008 
1009 }
1010 
1011 
Vector g(int i, int j)
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
HitTile_Vector global_graph
#define hit_tileNewType(baseType)
Definition: hit_tile.h:163
#define norm(vv)
#define hit_bShapeAddEdge2(shape, x, y)
Definition: hit_bshape.h:432
HitCom hit_comSparseUpdate(HitLayout lay, const void *tileP, HitType baseType)
Definition: hit_com.c:1084
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 hit_fileHBReadBitmap(hbfile)
Definition: hit_file.h:138
#define HIT_INT
Definition: hit_com.h:74
#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 hit_bShapeEdgeTarget(s, edge)
#define hit_bShapeEdgeIterator(var, shape, vertex)
Definition: hit_bshape.h:494
HitShape hit_bShapeExpand(HitShape shape, HitShape original, int amount)
Definition: hit_bshape.c:193
#define HIT_LAYOUT_NULL_STATIC
Definition: hit_layout.h:300
void hit_comFree(HitCom issue)
Definition: hit_com.c:1995
char * replace_str(char *str, const char *orig, const char *rep)
HitTile_int global_fixed
Matrix W(int i, int j)
void hit_comDo(HitCom *issue)
Definition: hit_com.c:2408
void hit_gbTileCopyVertices(void *destP, void *srcP)
Definition: hit_gbtile.c:159
void hit_gbTileClearVertices(void *varP)
Definition: hit_gbtile.c:146
HitTile_Matrix hessian
#define ITER_SJ
#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 clearV(vv)
HitTile_Vector graph
#define hit_comSparseScatter(lay, tilePSend, tilePRecv, baseType)
Definition: hit_com.h:852
#define HIT_EDGES
Definition: hit_tile.h:214
#define subV(r, a, b)
#define det(mm)
#define hit_gbTileVertexAt(var, vertex)
Definition: hit_gbtile.h:154
#define L
Definition: spring_bitmap.c:67
Vector solve(Matrix A, Vector b)
void free_structures()
#define hit_bShapeVertexIterator(var, shape)
Definition: hit_bshape.h:455
HitTile_int fixed
#define x
#define multSI(r, s)
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 clearM(mm)
HitShape ext_shape
#define ITER2
Definition: spring_bitmap.c:75
Vector solve_jacobi(Matrix A, Vector b)
#define hit_clockWorldReduce(c)
Definition: hit_utils.h:156
void hit_comInit(int *pargc, char **pargv[])
Definition: hit_com.c:111
HitShape shape
#define K
Definition: spring_bitmap.c:65
#define hit_bShapeNvertices(shape)
void random_coordinates()
#define hit_gbTileEdgeAt(var, pos1, pos2)
Definition: hit_gbtile.h:164
#define cpyV(r, a)
#define multSV(r, s, vv)
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)
#define multVVt(r, a, b)
HitTile_Vector new_e
HitShape HIT_BITMAP_SHAPE_NULL
Definition: hit_shape.c:71
#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 addM(r, a, b)
void calculate_GH(int vertex)
#define multMV(r, mm, vv)
#define v(a, b, c)
#define printfRoot(...)
Definition: spring_bitmap.c:59
#define hit_layout(name, topo,...)
Definition: hit_layout.h:415
#define inv(r, mm)
double gradient_norm()
HitLayout lay
Definition: heat.c:107
#define D
Definition: spring_bitmap.c:62
#define ITER1
Definition: spring_bitmap.c:73
#define hit_tileFree(var)
Definition: hit_tile.h:369
#define HIT_DOUBLE
Definition: hit_com.h:78
#define P_FIXED
Definition: spring_bitmap.c:70
void hit_comFinalize()
Definition: hit_com.c:159
#define hit_clockGetSeconds(c)
Definition: hit_utils.h:190
#define hit_gbTileDomainShapeAlloc(var, baseType, shape, allocOpts)
Definition: hit_gbtile.h:99
int iter1
Definition: spring_bitmap.c:78
double ** buffer
Definition: refMPIluBack.c:103
#define HIT_VERTICES
Definition: hit_tile.h:208
#define multSM(r, s, mm)
#define addV(r, a, b)
#define hit_gbTileEdgeIteratorAt(var, vertex, edge_index)
#define hit_clockStart(c)
Definition: hit_utils.h:87
double m[D][D]
Definition: spring_bitmap.c:93