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