Hitmap 1.3
 All Data Structures Namespaces Files Functions Variables Typedefs Friends Macros Groups Pages
spring_bitmap_iterskip.c
Go to the documentation of this file.
1 
13 /*
14  * <license>
15  *
16  * Hitmap v1.2
17  *
18  * This software is provided to enhance knowledge and encourage progress in the scientific
19  * community. It should be used only for research and educational purposes. Any reproduction
20  * or use for commercial purpose, public redistribution, in source or binary forms, with or
21  * without modifications, is NOT ALLOWED without the previous authorization of the copyright
22  * holder. The origin of this software must not be misrepresented; you must not claim that you
23  * wrote the original software. If you use this software for any purpose (e.g. publication),
24  * a reference to the software package and the authors must be included.
25  *
26  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS "AS IS" AND ANY
27  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
28  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
29  * THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
31  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
32  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
33  * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
34  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
35  *
36  * Copyright (c) 2007-2015, Trasgo Group, Universidad de Valladolid.
37  * All rights reserved.
38  *
39  * More information on http://trasgo.infor.uva.es/
40  *
41  * </license>
42 */
43 
44 #include <stdio.h>
45 #include <stdlib.h>
46 #include <hitmap.h>
47 #include <hit_com.h>
48 #include <unistd.h>
49 #include <string.h>
50 
51 
52 /* OUTPUT FUNCIONS */
53 #define printfRootInternal(...) { if( hit_Rank == 0 ) { printf(__VA_ARGS__); fflush(stdout); }}
54 #define printfRoot(...) printfRootInternal(__VA_ARGS__)
55 
57 #define D (3)
58 
60 #define K (1)
61 
62 #define L (0.001)
63 
65 #define P_FIXED 10
66 
68 #define ITER1 10
69 
70 #define ITER2 100
71 
73 int iter1 = ITER1;
75 int iter2 = ITER2;
76 
80 typedef struct {
81  double v[D];
82 } Vector;
83 
87 typedef struct {
88  double m[D][D];
89 } Matrix;
90 
91 
92 // Define the tile_Vector type.
94 
95 // Define the tile_Matrix type.
97 
98 // Define the tile_int type.
99 hit_tileNewType(int);
100 
101 // Define the tile_double type.
102 hit_tileNewType(double);
103 
104 
106 HitTile_Vector gradient;
107 
109 HitTile_Matrix hessian;
110 
112 HitTile_Vector e;
113 
115 HitTile_Vector new_e;
116 
119 
122 
125 
127 HitTile_Vector graph;
128 
130 HitTile_Vector global_graph;
131 
133 HitTile_int fixed;
134 
136 HitTile_int global_fixed;
137 
138 
141 
142 
147 void calculate_GH(int vertex);
148 
153 void solve_system_iter(int vertex);
154 
155 
156 #ifdef DEBUG
157 
161 void print_matrix(Matrix m){
162 
163  int i,j;
164  for(i=0;i<D;i++){
165  for(j=0;j<D;j++){
166  printf("%.2lf\t",m.m[i][j]);
167  }
168  printf("\n");
169  }
170 }
171 
176 void print_vector(Vector v){
177 
178  int i;
179  for(i=0;i<D;i++){
180  printf("%.2lf\t",v.v[i]);
181  }
182  printf("\n");
183 }
184 #endif
185 
186 
190 #define clearV(vv) clearVInternal((vv).v)
191 #define clearVInternal(v) \
192  ((v)[0] = (v)[1] = (v)[2] = 0.0 )
193 
197 #define cpyV(r,a) cpyVInternal((r).v,(a).v)
198 #define cpyVInternal(r,a) \
199  {(r)[0] = (a)[0]; (r)[1] = (a)[1]; (r)[2] = (a)[2]; }
200 
204 #define addV(r,a,b) addVInternal((r).v,(a).v,(b).v)
205 #define addVInternal(r,a,b) \
206  {(r)[0] = (a)[0] + (b)[0]; (r)[1] = (a)[1] + (b)[1]; (r)[2] = (a)[2] + (b)[2]; }
207 
208 
212 #define subV(r,a,b) subVInternal((r).v,(a).v,(b).v)
213 #define subVInternal(r,a,b) \
214  {(r)[0] = (a)[0] - (b)[0]; (r)[1] = (a)[1] - (b)[1]; (r)[2] = (a)[2] - (b)[2]; }
215 
219 #define norm(vv) normInternal((vv).v)
220 #define normInternal(v) \
221  (sqrt( pow((v)[0],2) + pow((v)[1],2) + pow((v)[2],2) ))
222 
226 #define multSV(r,s,vv) multSVInternal((r).v,s,(vv).v)
227 #define multSVInternal(r,s,v) \
228  {(r)[0] = (s) * (v)[0]; (r)[1] = (s) * (v)[1]; (r)[2] = (s) * (v)[2]; }
229 
233 #define clearM(mm) clearMInternal((mm).m)
234 #define clearMInternal(m) \
235  ((m)[0][0] = (m)[0][1] = (m)[0][2] = \
236  (m)[1][0] = (m)[1][1] = (m)[1][2] = \
237  (m)[2][0] = (m)[2][1] = (m)[2][2] = 0.0 )
238 
239 
249 #define det(mm) detInternal((mm).m)
250 #define detInternal(m) \
251  ( (m)[0][0] * (m)[1][1] * (m)[2][2] \
252  + (m)[0][1] * (m)[1][2] * (m)[2][0] \
253  + (m)[0][2] * (m)[1][0] * (m)[2][1] \
254  - (m)[0][0] * (m)[1][2] * (m)[2][1] \
255  - (m)[0][1] * (m)[1][0] * (m)[2][2] \
256  - (m)[0][2] * (m)[1][1] * (m)[2][0] )
257 
261 #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)
262 #define init_matrixInternal(m,a,b,c,d,e,f,g,h,i) \
263  {(m)[0][0] = a; (m)[0][1] = b; (m)[0][2] = c; \
264  (m)[1][0] = d; (m)[1][1] = e; (m)[1][2] = f; \
265  (m)[2][0] = g; (m)[2][1] = h; (m)[2][2] = i; }
266 
270 #define multSM(r,s,mm) multSMInternal((r).m,s,(mm).m)
271 #define multSMInternal(r,s,m) \
272  {(r)[0][0] = (s) * (m)[0][0]; (r)[0][1] = (s) * (m)[0][1]; (r)[0][2] = (s) * (m)[0][2]; \
273  (r)[1][0] = (s) * (m)[1][0]; (r)[1][1] = (s) * (m)[1][1]; (r)[1][2] = (s) * (m)[1][2]; \
274  (r)[2][0] = (s) * (m)[2][0]; (r)[2][1] = (s) * (m)[2][1]; (r)[2][2] = (s) * (m)[2][2]; }
275 
279 #define multSI(r,s) multSIInternal((r).m,s)
280 #define multSIInternal(r,s) \
281  {(r)[0][0] = (s); (r)[0][1] = 0; (r)[0][2] = 0; \
282  (r)[1][0] = 0; (r)[1][1] = (s); (r)[1][2] = 0; \
283  (r)[2][0] = 0; (r)[2][1] = 0; (r)[2][2] = (s); }
284 
285 
289 #define addM(r,a,b) addMInternal((r).m,(a).m,(b).m)
290 #define addMInternal(r,a,b) \
291  {(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]; \
292  (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]; \
293  (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]; }
294 
299 #define multVVt(r,a,b) multVVtInternal((r).m,(a).v,(b).v)
300 #define multVVtInternal(r,a,b) \
301  {(r)[0][0] = (a)[0] * (b)[0]; (r)[0][1] = (a)[0] * (b)[1]; (r)[0][2] = (a)[0] * (b)[2]; \
302  (r)[1][0] = (a)[1] * (b)[0]; (r)[1][1] = (a)[1] * (b)[1]; (r)[1][2] = (a)[1] * (b)[2]; \
303  (r)[2][0] = (a)[2] * (b)[0]; (r)[2][1] = (a)[2] * (b)[1]; (r)[2][2] = (a)[2] * (b)[2]; }
304 
305 
309 #define multMV(r,mm,vv) multMVInternal((r).v,(mm).m,(vv).v)
310 #define multMVInternal(r,m,v) \
311  {(r)[0] = (m)[0][0] * (v)[0] + (m)[0][1] * (v)[1] + (m)[0][2] * (v)[2]; \
312  (r)[1] = (m)[1][0] * (v)[0] + (m)[1][1] * (v)[1] + (m)[1][2] * (v)[2]; \
313  (r)[2] = (m)[2][0] * (v)[0] + (m)[2][1] * (v)[1] + (m)[2][2] * (v)[2]; }
314 
315 
319 #define inv(r,mm) invInternal((r).m,(mm).m)
320 #define invInternal(r,m) \
321  {(r)[0][0] = (m)[1][1] * (m)[2][2] - (m)[1][2] * (m)[2][1]; \
322  (r)[0][1] = - (m)[0][1] * (m)[2][2] + (m)[0][2] * (m)[2][1]; \
323  (r)[0][2] = (m)[0][1] * (m)[1][2] - (m)[0][2] * (m)[1][1]; \
324  (r)[1][0] = - (m)[1][0] * (m)[2][2] + (m)[1][2] * (m)[2][0]; \
325  (r)[1][1] = (m)[0][0] * (m)[2][2] - (m)[0][2] * (m)[2][0]; \
326  (r)[1][2] = - (m)[0][0] * (m)[1][2] + (m)[0][2] * (m)[1][0]; \
327  (r)[2][0] = (m)[1][0] * (m)[2][1] - (m)[1][1] * (m)[2][0]; \
328  (r)[2][1] = - (m)[0][0] * (m)[2][1] + (m)[0][1] * (m)[2][0]; \
329  (r)[2][2] = (m)[0][0] * (m)[1][1] - (m)[0][1] * (m)[1][0]; \
330  multSMInternal(r,1/detInternal(m),r); }
331 
335 #define multMM(r,a,b) multMMInternal((r).m,(a).m,(b).m)
336 #define multMMInternal(r,a,b) \
337  {(r)[0][0] = (a)[0][0] * (b)[0][0] + (a)[0][1] * (b)[1][0] + (a)[0][2] * (b)[2][0]; \
338  (r)[0][1] = (a)[0][0] * (b)[0][1] + (a)[0][1] * (b)[1][1] + (a)[0][2] * (b)[2][1]; \
339  (r)[0][2] = (a)[0][0] * (b)[0][2] + (a)[0][1] * (b)[1][2] + (a)[0][2] * (b)[2][2]; \
340  (r)[1][0] = (a)[1][0] * (b)[0][0] + (a)[1][1] * (b)[1][0] + (a)[1][2] * (b)[2][0]; \
341  (r)[1][1] = (a)[1][0] * (b)[0][1] + (a)[1][1] * (b)[1][1] + (a)[1][2] * (b)[2][1]; \
342  (r)[1][2] = (a)[1][0] * (b)[0][2] + (a)[1][1] * (b)[1][2] + (a)[1][2] * (b)[2][2]; \
343  (r)[2][0] = (a)[2][0] * (b)[0][0] + (a)[2][1] * (b)[1][0] + (a)[2][2] * (b)[2][0]; \
344  (r)[2][1] = (a)[2][0] * (b)[0][1] + (a)[2][1] * (b)[1][1] + (a)[2][2] * (b)[2][1]; \
345  (r)[2][2] = (a)[2][0] * (b)[0][2] + (a)[2][1] * (b)[1][2] + (a)[2][2] * (b)[2][2]; }
346 
347 
356 
357 // Number of iteration of the jacobi iterative method.
358 #define ITER_SJ 10
359 
360  Vector x;
361 
362  // The new approximation.
363  Vector x_1;
364 
365  // Clear the vectors.
366  clearV(x);
367  clearV(x_1);
368 
369  int iter;
370  int i,j;
371 
372  // Iterative method.
373  for(iter=0; iter<ITER_SJ; iter++){
374  for(i=0;i<D;i++){
375 
376  if( A.m[i][i] == 0 ) continue;
377 
378  // calculate the new value.
379  x_1.v[i] = 0.0;
380  for(j=0;j<D;j++){
381  if(i != j) x_1.v[i] += A.m[i][j] * x.v[j];
382  }
383  x_1.v[i] = (b.v[i] - x_1.v[i]) / A.m[i][i];
384 
385  }
386  cpyV(x,x_1);
387  }
388 
389  return x;
390 }
391 
392 
404 
405  Vector x;
406 
407  double detA = det(A);
408 
409  // If the matrix has a determinant, we can get the inverse.
410  if(detA != 0){
411 
412  Matrix InvA;
413  inv(InvA,A);
414  multMV(x,InvA,b);
415  } else {
416  //printf("#Warning: singular matrix\n");
417  x = solve_jacobi(A,b);
418  }
419  return x;
420 }
421 
426 
431 
433  hit_gbTileDomainShapeAlloc(&fixed, int, ext_shape, HIT_VERTICES);
434 
435 }
436 
441 
442  if(hit_Rank == 0){
445  }
446 
449 
450  hit_tileFree(e);
454 
455 }
456 
457 
462 void print_help(char * name){
463  printf("%s [-n NEWTON METHOD ITERATIONS] [-j JACOBI METHOD ITERATIONS] FILE \n",name);
464 
465  printf(" -n NEWTON METHOD ITERATIONS number of iterations (default 10)\n");
466  printf(" -j JACOBI METHOD ITERATIONS number of iterations (default 100)\n");
467  printf(" FILE input graph file\n");
468 
469 }
470 
478 char *replace_str(char *str, const char *orig, const char *rep){
479 
480  static char buffer[4096];
481  char *p;
482 
483  if(!(p = strstr(str, orig))) return NULL;
484 
485  strncpy(buffer, str, (size_t) (p-str)); // Copy characters from 'str' start to 'orig' st$
486  buffer[p-str] = '\0';
487 
488  sprintf(buffer+(p-str), "%s%s", rep, p+strlen(orig));
489 
490  return buffer;
491 }
492 
493 
494 
500 
501  srandom(0);
502  if(hit_Rank != 0) return;
503 
504  int n = hit_bShapeNvertices(global_shape);
505 
506  int j,i;
507 
508  // Set the coordinates
509  for(j=0;j<D;j++){
510  for(i=0;i<n;i++){
511  double number = (double) random();
512  hit_gbTileVertexAt(global_graph,i).v[j] = number;
513  }
514  }
515 
516  srandom(0);
517  for(i=0;i<n;i++){
518 
519  if ((random() % 100) < P_FIXED){
521  } else {
523  }
524  }
525 
526 }
527 
528 
529 
530 
531 
532 
536 void init_graph(int argc, char ** argv){
537 
538  // Define the spring structure.
539  global_shape = HIT_BITMAP_SHAPE_NULL;
540 
541  if(argc > 1){
542 
543  // Argument.
544  char c;
545 
546  // Parse the command-line arguments.
547  while ((c = (char) getopt (argc, argv, "hn:j:")) != -1)
548  switch (c) {
549  case 'n':
550  sscanf(optarg,"%d",&iter1);
551  break;
552  case 'j':
553  sscanf(optarg,"%d",&iter2);
554  break;
555  case 'h':
556  if(hit_Rank == 0) print_help(argv[0]);
557  hit_comFinalize();
558  exit(EXIT_SUCCESS);
559  default:
560  abort ();
561  break;
562  }
563 
564  if (hit_Rank != 0) return;
565 
566  // Graph file name.
567  char * graph_file = argv[optind];
568 
569  printfRoot("# Graph: %s\n",graph_file);
570  printfRoot("# Iterations %d x %d\n",iter1,iter2);
571 
572  global_shape = hit_fileHBReadBitmap(graph_file);
573 
574  printfRoot("# Vertices %d\n",hit_bShapeNvertices(global_shape));
575 
576 
577  // Allocate memory
580 
581  // Random coordinates
583 
584  int i;
585  int nfixed = 0;
586  for(i=0;i<hit_bShapeNvertices(global_shape);i++){
587  nfixed += hit_gbTileVertexAt(global_fixed,i);
588  }
589 
590  printf("# N fixed: %d\n",nfixed);
591 
592  } else {
593 
594  iter1 = iter2 = 3;
595 
596  if(hit_Rank != 0) return;
597 
598  hit_bShapeAddEdge2(&global_shape,0,1);
599  hit_bShapeAddEdge2(&global_shape,1,2);
600  hit_bShapeAddEdge2(&global_shape,2,3);
601  hit_bShapeAddEdge2(&global_shape,4,5);
602  hit_bShapeAddEdge2(&global_shape,5,6);
603  hit_bShapeAddEdge2(&global_shape,6,7);
604  hit_bShapeAddEdge2(&global_shape,0,4);
605  hit_bShapeAddEdge2(&global_shape,1,5);
606  hit_bShapeAddEdge2(&global_shape,2,6);
607  hit_bShapeAddEdge2(&global_shape,3,7);
608  hit_bShapeAddEdge2(&global_shape,1,4);
609  hit_bShapeAddEdge2(&global_shape,2,5);
610 
611  // Include each (i,i) edge.
612  hit_bShapeAddEdge2(&global_shape,0,0);
613  hit_bShapeAddEdge2(&global_shape,1,1);
614  hit_bShapeAddEdge2(&global_shape,2,2);
615  hit_bShapeAddEdge2(&global_shape,3,3);
616  hit_bShapeAddEdge2(&global_shape,4,4);
617  hit_bShapeAddEdge2(&global_shape,5,5);
618  hit_bShapeAddEdge2(&global_shape,6,6);
619  hit_bShapeAddEdge2(&global_shape,7,7);
620 
621  // Allocate memory
624 
625  // Set initial positions
626  Vector v0 = {{0,0,0}};
628  Vector v1 = {{100,0,0}};
630  Vector v2 = {{200,0,0}};
632  Vector v3 = {{300,0,0}};
634  Vector v4 = {{0,100,200}};
636  Vector v5 = {{100,100,0}};
638  Vector v6 = {{200,100,0}};
640  Vector v7 = {{300,100,0}};
642 
643  // Set the fix/free status
652 
653  }
654 }
655 
656 
657 
662 double gradient_norm(){
663 
664  HitTile_double norm, sum_norm;
665 
666  hit_tileDomainAlloc(&norm, double, 1,1);
667  hit_tileDomainAlloc(&sum_norm, double, 1,1);
668 
669  hit_tileElemAt(norm,1,0) = 0.0;
670 
671  int vertex;
672 
673  hit_bShapeVertexIterator(vertex,shape) {
674 
675  if (!hit_gbTileVertexAt(fixed,vertex)){
676 
677  calculate_GH(vertex);
678 
679  Vector Gi = hit_gbTileVertexAt(gradient,vertex);
680 
681  int d;
682  for(d=0;d<D;d++){
683  hit_tileElemAt(norm,1,0) += pow(Gi.v[d],2);
684  }
685  }
686  }
687 
688  HitOp op_sum;
690 
691 
692  HitRanks root = {{0,0,0,-1}};
693  HitCom com_sum = hit_comReduce(lay, root, &norm, &sum_norm, HIT_DOUBLE, op_sum);
694  hit_comDo(&com_sum);
695  hit_comFree(com_sum);
696 
697  double res = 0;
698 
699  if( hit_Rank == 0 ){
700  res = sqrt(hit_tileElemAt(sum_norm,1,0));
701  }
702 
703  hit_tileFree(norm);
704  hit_tileFree(sum_norm);
705 
706 
707  return res;
708 
709 }
710 
711 
712 
713 
717 int main(int argc, char ** argv) {
718 
719  int i;
720  int vertex;
721  HitClock init_time, comp_time;
722 
723  hit_comInit(&argc,&argv);
724 
725  hit_clockStart(init_time);
726  init_graph(argc,argv);
727 
728  // Create the topology object.
729  HitTopology topo = hit_topology(plug_topPlain);
730 
731  // Distribute the graph among the processors.
732  lay = hit_layout(plug_layBitmap,topo,&global_shape);
733 
734 
735  // Get shape and extended shape
736  shape = hit_layShape(lay);
737  ext_shape = hit_bShapeExpand(shape,global_shape,1);
738 
739  // Allocate memory of all the variables
740  init_structures();
741 
742  // Create the HitVector type.
743  HitType HitVector;
744  hit_comTypeStruct(&HitVector,Vector,1, v,3,HIT_DOUBLE);
745 
746  // Create the communication objects.
747  HitCom comA = hit_comSparseScatter(lay, &global_graph, &graph, HitVector);
749  HitCom comA3 = hit_comSparseUpdate(lay, &fixed, HIT_INT);
750  HitCom comB = hit_comSparseUpdate(lay, &graph, HitVector);
751  HitCom comC = hit_comSparseUpdate(lay, &e, HitVector);
752 
753  hit_clockStop(init_time);
754  hit_clockStart(comp_time);
755 
756  // Send all the data from the root proc to the other procs.
757  hit_comDo(&comA); hit_comDo(&comA2); hit_comDo(&comA3);
758  hit_comDo(&comB);
759 
760  // Calculate the inital gradient norm.
761  // double gnorm1 = gradient_norm();
762  // printfRoot("# Initial gradient norm: %lf\n",gnorm1);
763 
764  // Main loop for the newton method
765  for(i=0;i<iter1;i++){
766  printfRoot("# iter: %d/%d\n",i,iter1);
767 
768  // Iterate trough all the vertices and
769  // obtain the gradient and the hessian.
770  hit_bShapeVertexIterator(vertex,shape) {
771  if (!hit_gbTileVertexAt(fixed,vertex))
772  calculate_GH(vertex);
773  }
774 
775  // The initial position displacement is zero.
777 
778 
779  // Loop for the jacobi method to solve the system.
780  int j;
781  for(j=0;j<iter2;j++){
782 
783  // Perform a iteration.
784  hit_bShapeVertexIterator(vertex,shape) {
785  if (!hit_gbTileVertexAt(fixed,vertex))
786  solve_system_iter(vertex);
787  }
788 
789 
790  // Update the displacement.
792  hit_comDo(&comC);
793 
794  }
795 
796  // Update the position with the final displacement.
797  hit_bShapeVertexIterator(vertex,shape) {
798  if (!hit_gbTileVertexAt(fixed,vertex)){
799  Vector * current = &(hit_gbTileVertexAt(graph,vertex));
800  subV(*current,*current,hit_gbTileVertexAt(e,vertex));
801  }
802  }
803 
804  hit_comDo(&comB);
805 
806  }
807 
808 
809  hit_clockStop(comp_time);
810  hit_clockWorldReduce(init_time);
811  hit_clockWorldReduce(comp_time);
812 
813  double gnorm = gradient_norm();
814  printfRoot("# Final gradient norm: %lf\n",gnorm);
815  printfRoot("# Init time: %lf\n",hit_clockGetSeconds(init_time));
816  printfRoot("# Comp time: %lf\n",hit_clockGetSeconds(comp_time));
817  printfRoot("# Total time: %lf\n",hit_clockGetSeconds(init_time)+hit_clockGetSeconds(comp_time));
818 
819 
820  hit_comFree(comA);
821  hit_comFree(comB);
822  hit_comFree(comC);
823 
824  free_structures();
825  hit_layFree( lay );
826  hit_topFree( topo );
827  hit_comFinalize();
828 
829  return 0;
830 }
831 
832 
833 
834 
838 Vector g(int i, int j){
839 
840  Vector gij;
841 
842  if(i == j){
843  clearV(gij);
844  return gij;
845  }
846 
847  // Get the two vectors
848  Vector r_i = (hit_gbTileVertexAt(graph,i));
849  Vector r_j = (hit_gbTileVertexAt(graph,j));
850 
851  // Auxiliar variables.
852  Vector v;
853 
854  // v = r_i - r_j
855  subV(v,r_i,r_j);
856  // n = ||v||
857  double n = norm(v);
858 
859  if( n == 0){
860  clearV(gij);
861  return gij;
862  }
863 
864  // gij = - K * ((L - n) / n) * (v)
865  double scalar = - K * ((L - n) / n);
866  multSV(gij,scalar,v);
867 
868  return gij;
869 }
870 
871 
875 Matrix W(int i, int j){
876 
877  Matrix Wij;
878 
879  if(i == j){
880  clearM(Wij);
881  return Wij;
882  }
883 
884  // Get the two vectors
885  Vector r_i = (hit_gbTileVertexAt(graph,i));
886  Vector r_j = (hit_gbTileVertexAt(graph,j));
887 
888  // Auxiliary variables.
889  Vector v;
890  Matrix a, b;
891  double scalar_a, scalar_b;
892 
893  // v = r_i - r_j;
894  subV(v,r_i,r_j);
895  // n = ||v||
896  double n = norm(v);
897 
898  if( n == 0){
899  clearM(Wij);
900  return Wij;
901  }
902 
903  // h = -((L-n)/n) * eye(D) + (L/(n^3)) * (v * v')
904  scalar_a = -((L-n)/(n));
905  multSI(a,scalar_a);
906  scalar_b = (L/(pow(n,3)));
907  multVVt(b,v,v);
908  multSM(b,scalar_b,b);
909  addM(Wij,a,b);
910  // h = K * h
911  multSM(Wij,K,Wij);
912 
913  return Wij;
914 }
915 
916 
922 void calculate_GH(int vertex){
923 
924  Vector gs;
925  Matrix Ws;
926 
927  clearV(gs);
928  clearM(Ws);
929 
930  int edge;
931  hit_bShapeEdgeIteratorSkip(edge,ext_shape,vertex){
932 
933  // Get the neighbor.
934  int nghb_index = hit_bShapeEdgeTargetSkip(ext_shape,edge);
935 
936  // Calculate gradient[i]
937  Vector gij = g(vertex,nghb_index);
938  addV(gs,gs,gij);
939 
940  // Calculate part of the hessian[i][i]
941  Matrix Wij = W(vertex,nghb_index);
942  addM(Ws,Ws,Wij);
943 
944  // Calculate hessian[i][j] if j is variable.
945  if(!hit_gbTileVertexAt(fixed,nghb_index)){
946 
947  Matrix Hij = W(nghb_index,vertex);
948  multSM(Hij,-1,Hij);
949 
950  hit_gbTileEdgeIteratorSkipAt(hessian,vertex,edge) = Hij;
951  }
952  }
953 
954  hit_gbTileVertexAt(gradient,vertex) = gs;
955  hit_gbTileEdgeAt(hessian,vertex,vertex) = Ws;
956 
957 }
958 
964 void solve_system_iter(int vertex){
965 
966  Matrix Hii = hit_gbTileEdgeAt(hessian,vertex,vertex);
967  Vector Gi = hit_gbTileVertexAt(gradient,vertex);
968 
969  // Aux sum vector.
970  Vector sum;
971  clearV(sum);
972 
973  int edge;
974  hit_bShapeEdgeIteratorSkip(edge,ext_shape,vertex){
975 
976  // Get the neighbor.
977  int nghb_index = hit_bShapeEdgeTargetSkip(ext_shape,edge);
978  if(nghb_index == vertex) continue;
979 
980  // If the node is fixed, it is no part of the equation system.
981  if(hit_gbTileVertexAt(fixed,nghb_index)) continue;
982 
983  // Get the j displacement and ij hessian part.
984  Vector ej = hit_gbTileVertexAt(e,nghb_index);
985 
986  Matrix Hij = hit_gbTileEdgeIteratorSkipAt(hessian,vertex,edge);
987 
988  // Aux par vector.
989  Vector part;
990 
991  // Calculate and add the contribution of vertex j.
992  // SUM_{i not j} ( Hij e_j )
993  multMV(part,Hij,ej);
994  addV(sum,sum,part);
995 
996  }
997 
998  // Gi - SUM_{i not j} ( Hij e_j )
999  subV(sum,Gi,sum);
1000 
1001  // e_i = H_ii^(-1) * (Gi - SUM_{i not j} ( Hij e_j ) )
1002  hit_gbTileVertexAt(new_e,vertex) = solve(Hii,sum);
1003 
1004 }
1005 
1006 
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_bShapeEdgeTargetSkip(s, edge)
#define hit_layShape(lay)
Definition: hit_layout.h:650
HitTile_Vector global_graph
#define hit_tileNewType(baseType)
Definition: hit_tile.h:163
#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 D
#define hit_fileHBReadBitmap(hbfile)
Definition: hit_file.h:138
#define HIT_INT
Definition: hit_com.h:74
#define ITER2
#define hit_Rank
Definition: hit_com.h:140
#define clearM(mm)
#define cpyV(r, a)
int iter2
Definition: spring_bitmap.c:80
int iter
Definition: mmult_bit.c:68
#define ITER_SJ
#define hit_tileElemAt(var, ndims,...)
Definition: hit_tile.h:519
#define multVVt(r, a, b)
HitShape hit_bShapeExpand(HitShape shape, HitShape original, int amount)
Definition: hit_bshape.c:193
#define HIT_LAYOUT_NULL_STATIC
Definition: hit_layout.h:300
#define printfRoot(...)
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
#define multSI(r, s)
#define addM(r, a, b)
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
#define hit_bShapeEdgeIteratorSkip(var, shape, vertex)
Definition: hit_bshape.h:504
void hit_gbTileClearVertices(void *varP)
Definition: hit_gbtile.c:146
HitTile_Matrix hessian
#define hit_tileDomainAlloc(newVarP, baseType, numDims,...)
Definition: hit_tile.h:336
#define det(mm)
void solve_system_iter(int vertex)
#define hit_topology(name,...)
Definition: hit_topology.h:308
#define hit_gbTileEdgeIteratorSkipAt(var, vertex, edge_index)
Definition: hit_gbtile.h:206
#define K
HitTile_Vector graph
#define hit_comSparseScatter(lay, tilePSend, tilePRecv, baseType)
Definition: hit_com.h:852
#define HIT_EDGES
Definition: hit_tile.h:214
#define norm(vv)
#define hit_gbTileVertexAt(var, vertex)
Definition: hit_gbtile.h:154
Vector solve(Matrix A, Vector b)
void free_structures()
#define hit_bShapeVertexIterator(var, shape)
Definition: hit_bshape.h:455
HitTile_int fixed
#define x
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 multMV(r, mm, vv)
#define subV(r, a, b)
HitShape ext_shape
#define inv(r, mm)
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 hit_bShapeNvertices(shape)
#define multSM(r, s, mm)
#define L
void random_coordinates()
#define hit_gbTileEdgeAt(var, pos1, pos2)
Definition: hit_gbtile.h:164
#define P_FIXED
#define addV(r, a, b)
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 ITER1
#define m(a, b, c)
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 clearV(vv)
#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
void calculate_GH(int vertex)
#define v(a, b, c)
#define hit_layout(name, topo,...)
Definition: hit_layout.h:415
double gradient_norm()
HitLayout lay
Definition: heat.c:107
#define hit_tileFree(var)
Definition: hit_tile.h:369
#define HIT_DOUBLE
Definition: hit_com.h:78
void hit_comFinalize()
Definition: hit_com.c:159
#define hit_clockGetSeconds(c)
Definition: hit_utils.h:190
#define multSV(r, s, vv)
#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 hit_clockStart(c)
Definition: hit_utils.h:87
double m[D][D]
Definition: spring_bitmap.c:93