Hitmap 1.3
 All Data Structures Namespaces Files Functions Variables Typedefs Friends Macros Groups Pages
spring_csr_rows.c
Go to the documentation of this file.
1 
14 /*
15  * <license>
16  *
17  * Hitmap v1.2
18  *
19  * This software is provided to enhance knowledge and encourage progress in the scientific
20  * community. It should be used only for research and educational purposes. Any reproduction
21  * or use for commercial purpose, public redistribution, in source or binary forms, with or
22  * without modifications, is NOT ALLOWED without the previous authorization of the copyright
23  * holder. The origin of this software must not be misrepresented; you must not claim that you
24  * wrote the original software. If you use this software for any purpose (e.g. publication),
25  * a reference to the software package and the authors must be included.
26  *
27  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS "AS IS" AND ANY
28  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
29  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
30  * THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
32  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
33  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
34  * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36  *
37  * Copyright (c) 2007-2015, Trasgo Group, Universidad de Valladolid.
38  * All rights reserved.
39  *
40  * More information on http://trasgo.infor.uva.es/
41  *
42  * </license>
43 */
44 
45 #include <stdio.h>
46 #include <stdlib.h>
47 #include <hitmap.h>
48 #include <hit_com.h>
49 #include <unistd.h>
50 #include <string.h>
51 
52 
53 /* OUTPUT FUNCIONS */
54 #define printfRootInternal(...) { if( hit_Rank == 0 ) { printf(__VA_ARGS__); fflush(stdout); }}
55 #define printfRoot(...) printfRootInternal(__VA_ARGS__)
56 
58 #define D (3)
59 
61 #define K (1)
62 
63 #define L (0.001)
64 
66 #define P_FIXED 10
67 
69 #define ITER1 10
70 
71 #define ITER2 100
72 
74 int iter1 = ITER1;
76 int iter2 = ITER2;
77 
81 typedef struct {
82  double v[D];
83 } Vector;
84 
88 typedef struct {
89  double m[D][D];
90 } Matrix;
91 
92 
93 // Define the tile_Vector type.
95 
96 // Define the tile_Matrix type.
98 
99 // Define the tile_int type.
100 hit_tileNewType(int);
101 
102 // Define the tile_double type.
103 hit_tileNewType(double);
104 
105 
107 HitTile_Vector gradient;
108 
110 HitTile_Matrix hessian;
111 
113 HitTile_Vector e;
114 
116 HitTile_Vector new_e;
117 
120 
123 
126 
128 HitTile_Vector graph;
129 
131 HitTile_Vector global_graph;
132 
134 HitTile_int fixed;
135 
137 HitTile_int global_fixed;
138 
139 
142 
143 
148 void calculate_GH(int vertex);
149 
154 void solve_system_iter(int vertex);
155 
156 
157 #ifdef DEBUG
158 
162 void print_matrix(Matrix m){
163 
164  int i,j;
165  for(i=0;i<D;i++){
166  for(j=0;j<D;j++){
167  printf("%.2lf\t",m.m[i][j]);
168  }
169  printf("\n");
170  }
171 }
172 
177 void print_vector(Vector v){
178 
179  int i;
180  for(i=0;i<D;i++){
181  printf("%.2lf\t",v.v[i]);
182  }
183  printf("\n");
184 }
185 #endif
186 
187 
191 #define clearV(vv) clearVInternal((vv).v)
192 #define clearVInternal(v) \
193  ((v)[0] = (v)[1] = (v)[2] = 0.0 )
194 
198 #define cpyV(r,a) cpyVInternal((r).v,(a).v)
199 #define cpyVInternal(r,a) \
200  {(r)[0] = (a)[0]; (r)[1] = (a)[1]; (r)[2] = (a)[2]; }
201 
205 #define addV(r,a,b) addVInternal((r).v,(a).v,(b).v)
206 #define addVInternal(r,a,b) \
207  {(r)[0] = (a)[0] + (b)[0]; (r)[1] = (a)[1] + (b)[1]; (r)[2] = (a)[2] + (b)[2]; }
208 
209 
213 #define subV(r,a,b) subVInternal((r).v,(a).v,(b).v)
214 #define subVInternal(r,a,b) \
215  {(r)[0] = (a)[0] - (b)[0]; (r)[1] = (a)[1] - (b)[1]; (r)[2] = (a)[2] - (b)[2]; }
216 
220 #define norm(vv) normInternal((vv).v)
221 #define normInternal(v) \
222  (sqrt( pow((v)[0],2) + pow((v)[1],2) + pow((v)[2],2) ))
223 
227 #define multSV(r,s,vv) multSVInternal((r).v,s,(vv).v)
228 #define multSVInternal(r,s,v) \
229  {(r)[0] = (s) * (v)[0]; (r)[1] = (s) * (v)[1]; (r)[2] = (s) * (v)[2]; }
230 
234 #define clearM(mm) clearMInternal((mm).m)
235 #define clearMInternal(m) \
236  ((m)[0][0] = (m)[0][1] = (m)[0][2] = \
237  (m)[1][0] = (m)[1][1] = (m)[1][2] = \
238  (m)[2][0] = (m)[2][1] = (m)[2][2] = 0.0 )
239 
240 
250 #define det(mm) detInternal((mm).m)
251 #define detInternal(m) \
252  ( (m)[0][0] * (m)[1][1] * (m)[2][2] \
253  + (m)[0][1] * (m)[1][2] * (m)[2][0] \
254  + (m)[0][2] * (m)[1][0] * (m)[2][1] \
255  - (m)[0][0] * (m)[1][2] * (m)[2][1] \
256  - (m)[0][1] * (m)[1][0] * (m)[2][2] \
257  - (m)[0][2] * (m)[1][1] * (m)[2][0] )
258 
262 #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)
263 #define init_matrixInternal(m,a,b,c,d,e,f,g,h,i) \
264  {(m)[0][0] = a; (m)[0][1] = b; (m)[0][2] = c; \
265  (m)[1][0] = d; (m)[1][1] = e; (m)[1][2] = f; \
266  (m)[2][0] = g; (m)[2][1] = h; (m)[2][2] = i; }
267 
271 #define multSM(r,s,mm) multSMInternal((r).m,s,(mm).m)
272 #define multSMInternal(r,s,m) \
273  {(r)[0][0] = (s) * (m)[0][0]; (r)[0][1] = (s) * (m)[0][1]; (r)[0][2] = (s) * (m)[0][2]; \
274  (r)[1][0] = (s) * (m)[1][0]; (r)[1][1] = (s) * (m)[1][1]; (r)[1][2] = (s) * (m)[1][2]; \
275  (r)[2][0] = (s) * (m)[2][0]; (r)[2][1] = (s) * (m)[2][1]; (r)[2][2] = (s) * (m)[2][2]; }
276 
280 #define multSI(r,s) multSIInternal((r).m,s)
281 #define multSIInternal(r,s) \
282  {(r)[0][0] = (s); (r)[0][1] = 0; (r)[0][2] = 0; \
283  (r)[1][0] = 0; (r)[1][1] = (s); (r)[1][2] = 0; \
284  (r)[2][0] = 0; (r)[2][1] = 0; (r)[2][2] = (s); }
285 
286 
290 #define addM(r,a,b) addMInternal((r).m,(a).m,(b).m)
291 #define addMInternal(r,a,b) \
292  {(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]; \
293  (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]; \
294  (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]; }
295 
300 #define multVVt(r,a,b) multVVtInternal((r).m,(a).v,(b).v)
301 #define multVVtInternal(r,a,b) \
302  {(r)[0][0] = (a)[0] * (b)[0]; (r)[0][1] = (a)[0] * (b)[1]; (r)[0][2] = (a)[0] * (b)[2]; \
303  (r)[1][0] = (a)[1] * (b)[0]; (r)[1][1] = (a)[1] * (b)[1]; (r)[1][2] = (a)[1] * (b)[2]; \
304  (r)[2][0] = (a)[2] * (b)[0]; (r)[2][1] = (a)[2] * (b)[1]; (r)[2][2] = (a)[2] * (b)[2]; }
305 
306 
310 #define multMV(r,mm,vv) multMVInternal((r).v,(mm).m,(vv).v)
311 #define multMVInternal(r,m,v) \
312  {(r)[0] = (m)[0][0] * (v)[0] + (m)[0][1] * (v)[1] + (m)[0][2] * (v)[2]; \
313  (r)[1] = (m)[1][0] * (v)[0] + (m)[1][1] * (v)[1] + (m)[1][2] * (v)[2]; \
314  (r)[2] = (m)[2][0] * (v)[0] + (m)[2][1] * (v)[1] + (m)[2][2] * (v)[2]; }
315 
316 
320 #define inv(r,mm) invInternal((r).m,(mm).m)
321 #define invInternal(r,m) \
322  {(r)[0][0] = (m)[1][1] * (m)[2][2] - (m)[1][2] * (m)[2][1]; \
323  (r)[0][1] = - (m)[0][1] * (m)[2][2] + (m)[0][2] * (m)[2][1]; \
324  (r)[0][2] = (m)[0][1] * (m)[1][2] - (m)[0][2] * (m)[1][1]; \
325  (r)[1][0] = - (m)[1][0] * (m)[2][2] + (m)[1][2] * (m)[2][0]; \
326  (r)[1][1] = (m)[0][0] * (m)[2][2] - (m)[0][2] * (m)[2][0]; \
327  (r)[1][2] = - (m)[0][0] * (m)[1][2] + (m)[0][2] * (m)[1][0]; \
328  (r)[2][0] = (m)[1][0] * (m)[2][1] - (m)[1][1] * (m)[2][0]; \
329  (r)[2][1] = - (m)[0][0] * (m)[2][1] + (m)[0][1] * (m)[2][0]; \
330  (r)[2][2] = (m)[0][0] * (m)[1][1] - (m)[0][1] * (m)[1][0]; \
331  multSMInternal(r,1/detInternal(m),r); }
332 
336 #define multMM(r,a,b) multMMInternal((r).m,(a).m,(b).m)
337 #define multMMInternal(r,a,b) \
338  {(r)[0][0] = (a)[0][0] * (b)[0][0] + (a)[0][1] * (b)[1][0] + (a)[0][2] * (b)[2][0]; \
339  (r)[0][1] = (a)[0][0] * (b)[0][1] + (a)[0][1] * (b)[1][1] + (a)[0][2] * (b)[2][1]; \
340  (r)[0][2] = (a)[0][0] * (b)[0][2] + (a)[0][1] * (b)[1][2] + (a)[0][2] * (b)[2][2]; \
341  (r)[1][0] = (a)[1][0] * (b)[0][0] + (a)[1][1] * (b)[1][0] + (a)[1][2] * (b)[2][0]; \
342  (r)[1][1] = (a)[1][0] * (b)[0][1] + (a)[1][1] * (b)[1][1] + (a)[1][2] * (b)[2][1]; \
343  (r)[1][2] = (a)[1][0] * (b)[0][2] + (a)[1][1] * (b)[1][2] + (a)[1][2] * (b)[2][2]; \
344  (r)[2][0] = (a)[2][0] * (b)[0][0] + (a)[2][1] * (b)[1][0] + (a)[2][2] * (b)[2][0]; \
345  (r)[2][1] = (a)[2][0] * (b)[0][1] + (a)[2][1] * (b)[1][1] + (a)[2][2] * (b)[2][1]; \
346  (r)[2][2] = (a)[2][0] * (b)[0][2] + (a)[2][1] * (b)[1][2] + (a)[2][2] * (b)[2][2]; }
347 
348 
357 
358 // Number of iteration of the jacobi iterative method.
359 #define ITER_SJ 10
360 
361  Vector x;
362 
363  // The new approximation.
364  Vector x_1;
365 
366  // Clear the vectors.
367  clearV(x);
368  clearV(x_1);
369 
370  int iter;
371  int i,j;
372 
373  // Iterative method.
374  for(iter=0; iter<ITER_SJ; iter++){
375  for(i=0;i<D;i++){
376 
377  if( A.m[i][i] == 0 ) continue;
378 
379  // calculate the new value.
380  x_1.v[i] = 0.0;
381  for(j=0;j<D;j++){
382  if(i != j) x_1.v[i] += A.m[i][j] * x.v[j];
383  }
384  x_1.v[i] = (b.v[i] - x_1.v[i]) / A.m[i][i];
385 
386  }
387  cpyV(x,x_1);
388  }
389 
390  return x;
391 }
392 
393 
405 
406  Vector x;
407 
408  double detA = det(A);
409 
410  // If the matrix has a determinant, we can get the inverse.
411  if(detA != 0){
412 
413  Matrix InvA;
414  inv(InvA,A);
415  multMV(x,InvA,b);
416  } else {
417  //printf("#Warning: singular matrix\n");
418  x = solve_jacobi(A,b);
419  }
420  return x;
421 }
422 
427 
432 
434  hit_gcTileDomainShapeAlloc(&fixed, int, ext_shape, HIT_VERTICES);
435 
436 }
437 
442 
443  if(hit_Rank == 0){
446  }
447 
450 
451  hit_tileFree(e);
455 
456 }
457 
458 
463 void print_help(char * name){
464  printf("%s [-n NEWTON METHOD ITERATIONS] [-j JACOBI METHOD ITERATIONS] FILE \n",name);
465 
466  printf(" -n NEWTON METHOD ITERATIONS number of iterations (default 10)\n");
467  printf(" -j JACOBI METHOD ITERATIONS number of iterations (default 100)\n");
468  printf(" FILE input graph file\n");
469 
470 }
471 
479 char *replace_str(char *str, const char *orig, const char *rep){
480 
481  static char buffer[4096];
482  char *p;
483 
484  if(!(p = strstr(str, orig))) return NULL;
485 
486  strncpy(buffer, str, (size_t) (p-str)); // Copy characters from 'str' start to 'orig' st$
487  buffer[p-str] = '\0';
488 
489  sprintf(buffer+(p-str), "%s%s", rep, p+strlen(orig));
490 
491  return buffer;
492 }
493 
494 
495 
501 
502  srandom(0);
503  if(hit_Rank != 0) return;
504 
505  int n = hit_cShapeNvertices(global_shape);
506 
507  int j,i;
508 
509  // Set the coordinates
510  for(j=0;j<D;j++){
511  for(i=0;i<n;i++){
512  double number = (double) random();
513  hit_gcTileVertexAt(global_graph,i).v[j] = number;
514  }
515  }
516 
517  srandom(0);
518  for(i=0;i<n;i++){
519 
520  if ((random() % 100) < P_FIXED){
522  } else {
524  }
525  }
526 
527 }
528 
529 
530 
531 
532 
533 
537 void init_graph(int argc, char ** argv){
538 
539  // Define the spring structure.
540  global_shape = HIT_CSR_SHAPE_NULL;
541 
542  if(argc > 1){
543 
544  // Argument.
545  char c;
546 
547  // Parse the command-line arguments.
548  while ((c = (char) getopt (argc, argv, "hn:j:")) != -1)
549  switch (c) {
550  case 'n':
551  sscanf(optarg,"%d",&iter1);
552  break;
553  case 'j':
554  sscanf(optarg,"%d",&iter2);
555  break;
556  case 'h':
557  if(hit_Rank == 0) print_help(argv[0]);
558  hit_comFinalize();
559  exit(EXIT_SUCCESS);
560  default:
561  abort ();
562  break;
563  }
564 
565  if (hit_Rank != 0) return;
566 
567  // Graph file name.
568  char * graph_file = argv[optind];
569 
570  printfRoot("# Graph: %s\n",graph_file);
571  printfRoot("# Iterations %d x %d\n",iter1,iter2);
572 
573  global_shape = hit_fileHBRead(graph_file);
574 
575  printfRoot("# Vertices %d\n",hit_cShapeNvertices(global_shape));
576 
577 
578  // Allocate memory
581 
582  // Random coordinates
584 
585  int i;
586  int nfixed = 0;
587  for(i=0;i<hit_cShapeNvertices(global_shape);i++){
588  nfixed += hit_gcTileVertexAt(global_fixed,i);
589  }
590 
591  printf("# N fixed: %d\n",nfixed);
592 
593  } else {
594 
595  iter1 = iter2 = 3;
596 
597  if(hit_Rank != 0) return;
598 
599  hit_cShapeAddEdge2(&global_shape,0,1);
600  hit_cShapeAddEdge2(&global_shape,1,2);
601  hit_cShapeAddEdge2(&global_shape,2,3);
602  hit_cShapeAddEdge2(&global_shape,4,5);
603  hit_cShapeAddEdge2(&global_shape,5,6);
604  hit_cShapeAddEdge2(&global_shape,6,7);
605  hit_cShapeAddEdge2(&global_shape,0,4);
606  hit_cShapeAddEdge2(&global_shape,1,5);
607  hit_cShapeAddEdge2(&global_shape,2,6);
608  hit_cShapeAddEdge2(&global_shape,3,7);
609  hit_cShapeAddEdge2(&global_shape,1,4);
610  hit_cShapeAddEdge2(&global_shape,2,5);
611 
612  // Include each (i,i) edge.
613  hit_cShapeAddEdge2(&global_shape,0,0);
614  hit_cShapeAddEdge2(&global_shape,1,1);
615  hit_cShapeAddEdge2(&global_shape,2,2);
616  hit_cShapeAddEdge2(&global_shape,3,3);
617  hit_cShapeAddEdge2(&global_shape,4,4);
618  hit_cShapeAddEdge2(&global_shape,5,5);
619  hit_cShapeAddEdge2(&global_shape,6,6);
620  hit_cShapeAddEdge2(&global_shape,7,7);
621  hit_cShapeCard(global_shape,1) = 8; // @todo @javfres Previous functions should do this
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_cShapeVertexIterator(vertex,shape) {
676 
677  if (!hit_gcTileVertexAt(fixed,vertex)){
678 
679  calculate_GH(vertex);
680 
681  Vector Gi = hit_gcTileVertexAt(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_laySparseRows,topo,&global_shape);
735 
736  // Get shape and extended shape
737  shape = hit_layShape(lay);
738  ext_shape = hit_cShapeExpand(shape,global_shape,1);
739 
740  // Allocate memory of all the variables
741  init_structures();
742 
743  // Create the HitVector type.
744  HitType HitVector;
745  hit_comTypeStruct(&HitVector,Vector,1, v,3,HIT_DOUBLE);
746 
747  // Create the communication objects.
748  HitCom comA = hit_comSparseScatter(lay, &global_graph, &graph, HitVector);
750  HitCom comA3 = hit_comSparseUpdate(lay, &fixed, HIT_INT);
751  HitCom comB = hit_comSparseUpdate(lay, &graph, HitVector);
752  HitCom comC = hit_comSparseUpdate(lay, &e, HitVector);
753 
754  hit_clockStop(init_time);
755  hit_clockStart(comp_time);
756 
757  // Send all the data from the root proc to the other procs.
758  hit_comDo(&comA); hit_comDo(&comA2); hit_comDo(&comA3);
759  hit_comDo(&comB);
760 
761  // Calculate the inital gradient norm.
762  // double gnorm1 = gradient_norm();
763  // printfRoot("# Initial gradient norm: %lf\n",gnorm1);
764 
765  // Main loop for the newton method
766  for(i=0;i<iter1;i++){
767  printfRoot("# iter: %d/%d\n",i,iter1);
768 
769  // Iterate trough all the vertices and
770  // obtain the gradient and the hessian.
771  hit_cShapeVertexIterator(vertex,shape) {
772  if (!hit_gcTileVertexAt(fixed,vertex))
773  calculate_GH(vertex);
774  }
775 
776  // The initial position displacement is zero.
778 
779 
780  // Loop for the jacobi method to solve the system.
781  int j;
782  for(j=0;j<iter2;j++){
783 
784  // Perform a iteration.
785  hit_cShapeVertexIterator(vertex,shape) {
786  if (!hit_gcTileVertexAt(fixed,vertex))
787  solve_system_iter(vertex);
788  }
789 
790 
791  // Update the displacement.
793  hit_comDo(&comC);
794 
795  }
796 
797  // Update the position with the final displacement.
798  hit_cShapeVertexIterator(vertex,shape) {
799  if (!hit_gcTileVertexAt(fixed,vertex)){
800  Vector * current = &(hit_gcTileVertexAt(graph,vertex));
801  subV(*current,*current,hit_gcTileVertexAt(e,vertex));
802  }
803  }
804 
805  hit_comDo(&comB);
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(comA2);
822  hit_comFree(comA3);
823  hit_comFree(comB);
824  hit_comFree(comC);
825 
826  free_structures();
827 
828  hit_layFree( lay );
829 
830  hit_topFree( topo );
831  hit_comFinalize();
832 
833  return 0;
834 }
835 
836 
837 
838 
842 Vector g(int i, int j){
843 
844  Vector gij;
845 
846  if(i == j){
847  clearV(gij);
848  return gij;
849  }
850 
851  // Get the two vectors
852  Vector r_i = (hit_gcTileVertexAt(graph,i));
853  Vector r_j = (hit_gcTileVertexAt(graph,j));
854 
855  // Auxiliar variables.
856  Vector v;
857 
858  // v = r_i - r_j
859  subV(v,r_i,r_j);
860  // n = ||v||
861  double n = norm(v);
862 
863  if( n == 0){
864  clearV(gij);
865  return gij;
866  }
867 
868  // gij = - K * ((L - n) / n) * (v)
869  double scalar = - K * ((L - n) / n);
870  multSV(gij,scalar,v);
871 
872  return gij;
873 }
874 
875 
879 Matrix W(int i, int j){
880 
881  Matrix Wij;
882 
883  if(i == j){
884  clearM(Wij);
885  return Wij;
886  }
887 
888  // Get the two vectors
889  Vector r_i = (hit_gcTileVertexAt(graph,i));
890  Vector r_j = (hit_gcTileVertexAt(graph,j));
891 
892  // Auxiliary variables.
893  Vector v;
894  Matrix a, b;
895  double scalar_a, scalar_b;
896 
897  // v = r_i - r_j;
898  subV(v,r_i,r_j);
899  // n = ||v||
900  double n = norm(v);
901 
902  if( n == 0){
903  clearM(Wij);
904  return Wij;
905  }
906 
907  // h = -((L-n)/n) * eye(D) + (L/(n^3)) * (v * v')
908  scalar_a = -((L-n)/(n));
909  multSI(a,scalar_a);
910  scalar_b = (L/(pow(n,3)));
911  multVVt(b,v,v);
912  multSM(b,scalar_b,b);
913  addM(Wij,a,b);
914  // h = K * h
915  multSM(Wij,K,Wij);
916 
917  return Wij;
918 }
919 
920 
926 void calculate_GH(int vertex){
927 
928  Vector gs;
929  Matrix Ws;
930 
931  clearV(gs);
932  clearM(Ws);
933 
934  int edge;
935  hit_cShapeEdgeIterator(edge,ext_shape,vertex){
936 
937  // Get the neighbor.
938  int nghb_index = hit_cShapeEdgeTarget(ext_shape,edge);
939 
940  // Calculate gradient[i]
941  Vector gij = g(vertex,nghb_index);
942  addV(gs,gs,gij);
943 
944  // Calculate part of the hessian[i][i]
945  Matrix Wij = W(vertex,nghb_index);
946  addM(Ws,Ws,Wij);
947 
948  // Calculate hessian[i][j] if j is variable.
949  if(!hit_gcTileVertexAt(fixed,nghb_index)){
950 
951  Matrix Hij = W(nghb_index,vertex);
952  multSM(Hij,-1,Hij);
953 
954  hit_gcTileEdgeIteratorAt(hessian,vertex,edge) = Hij;
955  }
956  }
957 
958  hit_gcTileVertexAt(gradient,vertex) = gs;
959  hit_gcTileEdgeAt(hessian,vertex,vertex) = Ws;
960 
961 }
962 
968 void solve_system_iter(int vertex){
969 
970  Matrix Hii = hit_gcTileEdgeAt(hessian,vertex,vertex);
971  Vector Gi = hit_gcTileVertexAt(gradient,vertex);
972 
973  // Aux sum vector.
974  Vector sum;
975  clearV(sum);
976 
977  int edge;
978  hit_cShapeEdgeIterator(edge,ext_shape,vertex){
979 
980  // Get the neighbor.
981  int nghb_index = hit_cShapeEdgeTarget(ext_shape,edge);
982  if(nghb_index == vertex) continue;
983 
984  // If the node is fixed, it is no part of the equation system.
985  if(hit_gcTileVertexAt(fixed,nghb_index)) continue;
986 
987  // Get the j displacement and ij hessian part.
988  Vector ej = hit_gcTileVertexAt(e,nghb_index);
989 
990  Matrix Hij = hit_gcTileEdgeIteratorAt(hessian,vertex,edge);
991 
992  // Aux par vector.
993  Vector part;
994 
995  // Calculate and add the contribution of vertex j.
996  // SUM_{i not j} ( Hij e_j )
997  multMV(part,Hij,ej);
998  addV(sum,sum,part);
999 
1000  }
1001 
1002  // Gi - SUM_{i not j} ( Hij e_j )
1003  subV(sum,Gi,sum);
1004 
1005  // e_i = H_ii^(-1) * (Gi - SUM_{i not j} ( Hij e_j ) )
1006  hit_gcTileVertexAt(new_e,vertex) = solve(Hii,sum);
1007 
1008 }
1009 
#define multMV(r, mm, vv)
#define P_FIXED
Vector g(int i, int j)
MPI_Op HitOp
Definition: hit_com.h:118
void hit_gcTileCopyVertices(void *destP, void *srcP)
Definition: hit_gctile.c:145
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_cShapeEdgeTarget(s, edge)
HitCom hit_comSparseUpdate(HitLayout lay, const void *tileP, HitType baseType)
Definition: hit_com.c:1084
#define K
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 ITER_SJ
#define multSM(r, s, mm)
#define HIT_INT
Definition: hit_com.h:74
#define hit_Rank
Definition: hit_com.h:140
#define norm(vv)
int iter2
Definition: spring_bitmap.c:80
#define hit_cShapeVertexIterator(var, shape)
Definition: hit_cshape.h:392
int iter
Definition: mmult_bit.c:68
#define hit_tileElemAt(var, ndims,...)
Definition: hit_tile.h:519
#define HIT_LAYOUT_NULL_STATIC
Definition: hit_layout.h:300
#define D
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
HitTile_Matrix hessian
#define ITER2
#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_cShapeEdgeIterator(var, shape, vertex)
Definition: hit_cshape.h:416
#define hit_gcTileEdgeIteratorAt(var, vertex, edge_index)
HitTile_Vector graph
#define hit_comSparseScatter(lay, tilePSend, tilePRecv, baseType)
Definition: hit_com.h:852
#define HIT_EDGES
Definition: hit_tile.h:214
#define hit_cShapeCard(shape, dim)
Definition: hit_cshape.h:108
#define clearV(vv)
void hit_gcTileClearVertices(void *varP)
Definition: hit_gctile.c:132
#define subV(r, a, b)
#define multVVt(r, a, b)
#define ITER1
#define printfRoot(...)
Vector solve(Matrix A, Vector b)
void free_structures()
#define hit_gcTileDomainShapeAlloc(var, baseType, shape, allocOpts)
Definition: hit_gctile.h:99
#define L
#define cpyV(r, a)
HitTile_int fixed
#define x
#define hit_gcTileVertexAt(var, vertex)
Definition: hit_gctile.h:157
HitShape hit_cShapeExpand(HitShape shape, HitShape original, int amount)
Definition: hit_cshape.c:208
#define multSV(r, s, vv)
#define hit_cShapeNvertices(shape)
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
void hit_comInit(int *pargc, char **pargv[])
Definition: hit_com.c:111
HitShape shape
#define addM(r, a, b)
void random_coordinates()
HitShape HIT_CSR_SHAPE_NULL
Definition: hit_shape.c:70
#define hit_fileHBRead(hbfile)
Definition: hit_file.h:77
void hit_layFree(HitLayout lay)
Definition: hit_layout.c:2251
#define hit_clockStop(c)
Definition: hit_utils.h:109
#define det(mm)
#define addV(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
#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()
#define hit_cShapeAddEdge2(shape, x, y)
Definition: hit_cshape.h:358
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 inv(r, mm)
#define clearM(mm)
#define HIT_DOUBLE
Definition: hit_com.h:78
#define hit_gcTileEdgeAt(var, pos1, pos2)
Definition: hit_gctile.h:167
void hit_comFinalize()
Definition: hit_com.c:159
#define hit_clockGetSeconds(c)
Definition: hit_utils.h:190
int iter1
Definition: spring_bitmap.c:78
#define multSI(r, s)
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