Hitmap 1.3
 All Data Structures Namespaces Files Functions Variables Typedefs Friends Macros Groups Pages
spring_dense.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 <time.h>
50 #include <string.h>
51 #include <sys/time.h>
52 
53 
54 
55 
56 /* OUTPUT FUNCIONS */
57 #define printfRootInternal(...) { if( hit_Rank == 0 ) { printf(__VA_ARGS__); fflush(stdout); }}
58 #define printfRoot(...) printfRootInternal(__VA_ARGS__)
59 
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 
99 typedef struct {
101  int fixed;
102 } Node;
103 
104 
105 
106 // Define the tile_Node.
108 
109 // Define the tile_VectorA.
111 
112 // Define the tile_double
113 hit_tileNewType(double);
114 
115 // Define the tile_Matrix
117 
118 hit_tileNewType(int);
119 
120 
122 HitTile_int incident;
123 
124 
125 
129 
130 
132 HitTile_Vector gradient;
133 
135 HitTile_Matrix hessian;
136 
138 HitTile_Vector e;
139 
141 HitTile_Vector new_e;
142 
145 
147 HitTile_Node graph;
148 
149 
150 
153 
154 
155 
160 void calculate_GH(int vertex);
161 
166 void solve_system_iter(int vertex);
167 
168 
169 #ifdef DEBUG
170 
174 void print_matrix(Matrix m){
175 
176  int i,j;
177  for(i=0;i<D;i++){
178  for(j=0;j<D;j++){
179  printf("%.2lf\t",m.m[i][j]);
180  }
181  printf("\n");
182  }
183 }
184 
189 void print_vector(Vector v){
190 
191  int i;
192  for(i=0;i<D;i++){
193  printf("%.2lf\t",v.v[i]);
194  }
195  printf("\n");
196 }
197 #endif
198 
199 
203 #define clearV(vv) clearVInternal((vv).v)
204 #define clearVInternal(v) \
205  ((v)[0] = (v)[1] = (v)[2] = 0.0 )
206 
210 #define cpyV(r,a) cpyVInternal((r).v,(a).v)
211 #define cpyVInternal(r,a) \
212  {(r)[0] = (a)[0]; (r)[1] = (a)[1]; (r)[2] = (a)[2]; }
213 
217 #define addV(r,a,b) addVInternal((r).v,(a).v,(b).v)
218 #define addVInternal(r,a,b) \
219  {(r)[0] = (a)[0] + (b)[0]; (r)[1] = (a)[1] + (b)[1]; (r)[2] = (a)[2] + (b)[2]; }
220 
221 
225 #define subV(r,a,b) subVInternal((r).v,(a).v,(b).v)
226 #define subVInternal(r,a,b) \
227  {(r)[0] = (a)[0] - (b)[0]; (r)[1] = (a)[1] - (b)[1]; (r)[2] = (a)[2] - (b)[2]; }
228 
232 #define norm(vv) normInternal((vv).v)
233 #define normInternal(v) \
234  (sqrt( pow((v)[0],2) + pow((v)[1],2) + pow((v)[2],2) ))
235 
239 #define multSV(r,s,vv) multSVInternal((r).v,s,(vv).v)
240 #define multSVInternal(r,s,v) \
241  {(r)[0] = (s) * (v)[0]; (r)[1] = (s) * (v)[1]; (r)[2] = (s) * (v)[2]; }
242 
246 #define clearM(mm) clearMInternal((mm).m)
247 #define clearMInternal(m) \
248  ((m)[0][0] = (m)[0][1] = (m)[0][2] = \
249  (m)[1][0] = (m)[1][1] = (m)[1][2] = \
250  (m)[2][0] = (m)[2][1] = (m)[2][2] = 0.0 )
251 
252 
262 #define det(mm) detInternal((mm).m)
263 #define detInternal(m) \
264  ( (m)[0][0] * (m)[1][1] * (m)[2][2] \
265  + (m)[0][1] * (m)[1][2] * (m)[2][0] \
266  + (m)[0][2] * (m)[1][0] * (m)[2][1] \
267  - (m)[0][0] * (m)[1][2] * (m)[2][1] \
268  - (m)[0][1] * (m)[1][0] * (m)[2][2] \
269  - (m)[0][2] * (m)[1][1] * (m)[2][0] )
270 
274 #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)
275 #define init_matrixInternal(m,a,b,c,d,e,f,g,h,i) \
276  {(m)[0][0] = a; (m)[0][1] = b; (m)[0][2] = c; \
277  (m)[1][0] = d; (m)[1][1] = e; (m)[1][2] = f; \
278  (m)[2][0] = g; (m)[2][1] = h; (m)[2][2] = i; }
279 
283 #define multSM(r,s,mm) multSMInternal((r).m,s,(mm).m)
284 #define multSMInternal(r,s,m) \
285  {(r)[0][0] = (s) * (m)[0][0]; (r)[0][1] = (s) * (m)[0][1]; (r)[0][2] = (s) * (m)[0][2]; \
286  (r)[1][0] = (s) * (m)[1][0]; (r)[1][1] = (s) * (m)[1][1]; (r)[1][2] = (s) * (m)[1][2]; \
287  (r)[2][0] = (s) * (m)[2][0]; (r)[2][1] = (s) * (m)[2][1]; (r)[2][2] = (s) * (m)[2][2]; }
288 
292 #define multSI(r,s) multSIInternal((r).m,s)
293 #define multSIInternal(r,s) \
294  {(r)[0][0] = (s); (r)[0][1] = 0; (r)[0][2] = 0; \
295  (r)[1][0] = 0; (r)[1][1] = (s); (r)[1][2] = 0; \
296  (r)[2][0] = 0; (r)[2][1] = 0; (r)[2][2] = (s); }
297 
298 
302 #define addM(r,a,b) addMInternal((r).m,(a).m,(b).m)
303 #define addMInternal(r,a,b) \
304  {(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]; \
305  (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]; \
306  (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]; }
307 
312 #define multVVt(r,a,b) multVVtInternal((r).m,(a).v,(b).v)
313 #define multVVtInternal(r,a,b) \
314  {(r)[0][0] = (a)[0] * (b)[0]; (r)[0][1] = (a)[0] * (b)[1]; (r)[0][2] = (a)[0] * (b)[2]; \
315  (r)[1][0] = (a)[1] * (b)[0]; (r)[1][1] = (a)[1] * (b)[1]; (r)[1][2] = (a)[1] * (b)[2]; \
316  (r)[2][0] = (a)[2] * (b)[0]; (r)[2][1] = (a)[2] * (b)[1]; (r)[2][2] = (a)[2] * (b)[2]; }
317 
318 
322 #define multMV(r,mm,vv) multMVInternal((r).v,(mm).m,(vv).v)
323 #define multMVInternal(r,m,v) \
324  {(r)[0] = (m)[0][0] * (v)[0] + (m)[0][1] * (v)[1] + (m)[0][2] * (v)[2]; \
325  (r)[1] = (m)[1][0] * (v)[0] + (m)[1][1] * (v)[1] + (m)[1][2] * (v)[2]; \
326  (r)[2] = (m)[2][0] * (v)[0] + (m)[2][1] * (v)[1] + (m)[2][2] * (v)[2]; }
327 
328 
332 #define inv(r,mm) invInternal((r).m,(mm).m)
333 #define invInternal(r,m) \
334  {(r)[0][0] = (m)[1][1] * (m)[2][2] - (m)[1][2] * (m)[2][1]; \
335  (r)[0][1] = - (m)[0][1] * (m)[2][2] + (m)[0][2] * (m)[2][1]; \
336  (r)[0][2] = (m)[0][1] * (m)[1][2] - (m)[0][2] * (m)[1][1]; \
337  (r)[1][0] = - (m)[1][0] * (m)[2][2] + (m)[1][2] * (m)[2][0]; \
338  (r)[1][1] = (m)[0][0] * (m)[2][2] - (m)[0][2] * (m)[2][0]; \
339  (r)[1][2] = - (m)[0][0] * (m)[1][2] + (m)[0][2] * (m)[1][0]; \
340  (r)[2][0] = (m)[1][0] * (m)[2][1] - (m)[1][1] * (m)[2][0]; \
341  (r)[2][1] = - (m)[0][0] * (m)[2][1] + (m)[0][1] * (m)[2][0]; \
342  (r)[2][2] = (m)[0][0] * (m)[1][1] - (m)[0][1] * (m)[1][0]; \
343  multSMInternal(r,1/detInternal(m),r); }
344 
348 #define multMM(r,a,b) multMMInternal((r).m,(a).m,(b).m)
349 #define multMMInternal(r,a,b) \
350  {(r)[0][0] = (a)[0][0] * (b)[0][0] + (a)[0][1] * (b)[1][0] + (a)[0][2] * (b)[2][0]; \
351  (r)[0][1] = (a)[0][0] * (b)[0][1] + (a)[0][1] * (b)[1][1] + (a)[0][2] * (b)[2][1]; \
352  (r)[0][2] = (a)[0][0] * (b)[0][2] + (a)[0][1] * (b)[1][2] + (a)[0][2] * (b)[2][2]; \
353  (r)[1][0] = (a)[1][0] * (b)[0][0] + (a)[1][1] * (b)[1][0] + (a)[1][2] * (b)[2][0]; \
354  (r)[1][1] = (a)[1][0] * (b)[0][1] + (a)[1][1] * (b)[1][1] + (a)[1][2] * (b)[2][1]; \
355  (r)[1][2] = (a)[1][0] * (b)[0][2] + (a)[1][1] * (b)[1][2] + (a)[1][2] * (b)[2][2]; \
356  (r)[2][0] = (a)[2][0] * (b)[0][0] + (a)[2][1] * (b)[1][0] + (a)[2][2] * (b)[2][0]; \
357  (r)[2][1] = (a)[2][0] * (b)[0][1] + (a)[2][1] * (b)[1][1] + (a)[2][2] * (b)[2][1]; \
358  (r)[2][2] = (a)[2][0] * (b)[0][2] + (a)[2][1] * (b)[1][2] + (a)[2][2] * (b)[2][2]; }
359 
360 
369 
370 // Number of iteration of the jacobi iterative method.
371 #define ITER_SJ 10
372 
373  Vector x;
374 
375  // The new approximation.
376  Vector x_1;
377 
378  // Clear the vectors.
379  clearV(x);
380  clearV(x_1);
381 
382  int iter;
383  int i,j;
384 
385  // Iterative method.
386  for(iter=0; iter<ITER_SJ; iter++){
387  for(i=0;i<D;i++){
388 
389  if( A.m[i][i] == 0 ) continue;
390 
391  // calculate the new value.
392  x_1.v[i] = 0.0;
393  for(j=0;j<D;j++){
394  if(i != j) x_1.v[i] += A.m[i][j] * x.v[j];
395  }
396  x_1.v[i] = (b.v[i] - x_1.v[i]) / A.m[i][i];
397 
398  }
399  cpyV(x,x_1);
400  }
401 
402  return x;
403 }
404 
405 
417 
418  Vector x;
419 
420  double detA = det(A);
421 
422  // If the matrix has a determinant, we can get the inverse.
423  if(detA != 0){
424 
425  Matrix InvA;
426  inv(InvA,A);
427  multMV(x,InvA,b);
428  } else {
429  //printf("#Warning: singular matrix\n");
430  x = solve_jacobi(A,b);
431  }
432  return x;
433 }
434 
435 
440 
441 
442  if(hit_Rank != 0){
443  hit_tileDomainShapeAlloc(&graph,Node,global_shape);
444  }
445 
446  hit_tileDomainShapeAlloc(&e,Vector,global_shape);
447  hit_tileDomainShapeAlloc(&new_e,Vector,global_shape);
448 
450  hit_tileDomainShapeAlloc(&hessian,Matrix,incident_shape);
451 }
452 
453 
454 
459 
461  hit_tileFree(e);
465 
466 }
467 
468 
469 
470 
471 
476 void print_help(char * name){
477  printf("%s [-n NEWTON METHOD ITERATIONS] [-j JACOBI METHOD ITERATIONS] FILE \n",name);
478 
479  printf(" -n NEWTON METHOD ITERATIONS number of iterations (default 10)\n");
480  printf(" -j JACOBI METHOD ITERATIONS number of iterations (default 100)\n");
481  printf(" FILE input graph file\n");
482 
483 }
484 
485 
486 
494 char *replace_str(char *str, const char *orig, const char *rep){
495 
496  static char buffer[4096];
497  char *p;
498 
499  if(!(p = strstr(str, orig))) return NULL;
500 
501  strncpy(buffer, str, (size_t) (p-str)); // Copy characters from 'str' start to 'orig' st$
502  buffer[p-str] = '\0';
503 
504  sprintf(buffer+(p-str), "%s%s", rep, p+strlen(orig));
505 
506  return buffer;
507 }
508 
509 
510 
511 
516 void random_coordinates(HitShape shape_rc, HitTile_Node graph_rc){
517 
518  srandom(0);
519  if(hit_Rank != 0) return;
520 
521  int n = hit_sigCard(hit_shapeSig(shape_rc,0));
522 
523  int j,i;
524 
525  // Read the coordinates
526  for(j=0;j<D;j++){
527  for(i=0;i<n;i++){
528 
529  double number = (double) random();
530  hit_tileElemAt(graph_rc,1,i).r.v[j] = number;
531 
532  }
533  }
534 
535 
536  srandom(0);
537 
538  for(i=0;i<n;i++){
539 
540  if ((random() % 100) < P_FIXED){
541  hit_tileElemAt(graph_rc,1,i).fixed = 1;
542  } else {
543  hit_tileElemAt(graph_rc,1,i).fixed = 0;
544  }
545  }
546 
547 }
548 
549 
550 
551 
552 /*
553 void read_coordinates(char * name, HitShape shape_rc, HitTile_Node graph_rc){
554 
555  if(hit_Rank != 0) return;
556 
557 
558  FILE * file = fopen(name,"r");
559 
560 
561  if(file == NULL){
562  printf("File not exists\n");
563  return;
564  }
565 
566 
567 
568  char c = (char) getc(file);
569 
570  // Skip first comments
571  if(c == '%'){
572  while( c != '\n' ) c = (char) getc(file);
573  }
574  ungetc(c,file);
575 
576  int res,n,d,i,j;
577 
578  int n_shape = hit_sigCard(hit_shapeSig(shape_rc,0));
579  res = fscanf(file,"%d %d\n",&n,&d);
580 
581  if(n != n_shape) printf("Coordinates: Number of vertices do not match %d\n",n); n = min(n,n_shape);
582  if(d != D) printf("Coordinates: Number of dimensions do not match %d\n",d); d = min(d,D);
583 
584  srandom(0);
585 
586 
587  // Read the coordinates
588  for(j=0;j<D;j++){
589  for(i=0;i<n;i++){
590  if(j<d){
591  double number;
592  res = fscanf(file,"%lf\n",&number);
593  hit_tileElemAt(graph_rc,1,i).r.v[j] = number;
594  } else {
595  hit_tileElemAt(graph_rc,1,i).r.v[j] = 0.0;
596  }
597  }
598  }
599 
600 
601  for(i=0;i<n;i++){
602 
603  if ((random() % 100) < P_FIXED){
604  hit_tileElemAt(graph_rc,1,i).fixed = 1;
605  } else {
606  hit_tileElemAt(graph_rc,1,i).fixed = 0;
607  }
608 
609  }
610 
611  fclose(file);
612 
613 }
614 */
615 
616 
617 
618 
619 
620 
621 
622 
626 void init_graph(int argc, char ** argv){
627 
628 
629  if(argc > 1){
630 
631  // Argument.
632  char c;
633 
634  // Parse the command-line arguments.
635  while ((c = (char) getopt (argc, argv, "hn:j:")) != -1)
636  switch (c) {
637  case 'n':
638  sscanf(optarg,"%d",&iter1);
639  break;
640  case 'j':
641  sscanf(optarg,"%d",&iter2);
642  break;
643  case 'h':
644  if(hit_Rank == 0) print_help(argv[0]);
645  hit_comFinalize();
646  exit(EXIT_SUCCESS);
647  default:
648  abort();
649  break;
650  }
651 
652  if (hit_Rank != 0) return;
653 
654  // Graph file name.
655  char * graph_file = argv[optind];
656 
657  printfRoot("# Graph: %s\n",graph_file);
658  printfRoot("# Iterations %dx%d\n",iter1,iter2);
659 
660  Nvertices = hit_fileHBVertices(graph_file);
661 
662  printfRoot("# Vertices %d\n",Nvertices);
663 
664 
665  HitSig sig = hit_sig(0,Nvertices-1,1);
666  incident_shape = hit_shape(2,sig,sig);
667 
668  hit_tileDomainShapeAlloc(&incident,int,incident_shape);
669 
670 
671  hit_fileHBReadDense(graph_file,&incident);
672 
673 #ifdef DEBUG
674  {int i,j;
675  for(i=0;i<Nvertices;i++){
676  for(j=0;j<Nvertices;j++){
677  if(hit_tileElemAt2(incident,i,j))
678  printf("X");
679  else
680  printf(" ");
681  }
682  printf("\n");
683  }}
684 #endif
685 
686 
687  global_shape = hit_shape(1,hit_sig(0,Nvertices-1,1));
688  hit_tileDomainShapeAlloc(&graph,Node,global_shape);
689 
690  // Coordinates file
691  char * coord_file = replace_str(graph_file,".rb","_coord.mtx");
692  printfRoot("# Coordinates: %s\n",coord_file);
693 
694  // Random coordinates
695  random_coordinates(global_shape,graph);
696 
697  //read_coordinates(coord_file,global_shape,graph);
698 
699 
700  int i;
701  int nfixed = 0;
702  for(i=0;i<Nvertices;i++){
703 
704  nfixed += hit_tileElemAt(graph,1,i).fixed;
705  }
706 
707  printf("# N fixed: %d\n",nfixed);
708 
709  } else {
710 
711  iter1 = iter2 = 3;
712  Nvertices = 8;
713  HitSig sig = hit_sig(0,Nvertices-1,1);
714  incident_shape = hit_shape(2,sig,sig);
715 
716  if(hit_Rank != 0) return;
717 
718  hit_tileDomainShapeAlloc(&incident,int,incident_shape);
719 
720  int i,j;
721 
722  for(i=0;i<Nvertices;i++)
723  for(j=0;j<Nvertices;j++)
724  hit_tileElemAt2(incident,i,j) = 0;
725 
726  hit_tileElemAt2(incident,0,1) = 1;
727  hit_tileElemAt2(incident,1,2) = 1;
728  hit_tileElemAt2(incident,2,3) = 1;
729  hit_tileElemAt2(incident,4,5) = 1;
730  hit_tileElemAt2(incident,5,6) = 1;
731  hit_tileElemAt2(incident,6,7) = 1;
732  hit_tileElemAt2(incident,0,4) = 1;
733  hit_tileElemAt2(incident,1,5) = 1;
734  hit_tileElemAt2(incident,2,6) = 1;
735  hit_tileElemAt2(incident,3,7) = 1;
736  hit_tileElemAt2(incident,1,4) = 1;
737  hit_tileElemAt2(incident,2,5) = 1;
738 
739  hit_tileElemAt2(incident,1,0) = 1;
740  hit_tileElemAt2(incident,2,1) = 1;
741  hit_tileElemAt2(incident,3,2) = 1;
742  hit_tileElemAt2(incident,5,4) = 1;
743  hit_tileElemAt2(incident,6,5) = 1;
744  hit_tileElemAt2(incident,7,6) = 1;
745  hit_tileElemAt2(incident,4,0) = 1;
746  hit_tileElemAt2(incident,5,1) = 1;
747  hit_tileElemAt2(incident,6,2) = 1;
748  hit_tileElemAt2(incident,7,3) = 1;
749  hit_tileElemAt2(incident,4,1) = 1;
750  hit_tileElemAt2(incident,5,2) = 1;
751 
752  // Include each (i,i) edge.
753  hit_tileElemAt2(incident,0,0) = 1;
754  hit_tileElemAt2(incident,1,1) = 1;
755  hit_tileElemAt2(incident,2,2) = 1;
756  hit_tileElemAt2(incident,3,3) = 1;
757  hit_tileElemAt2(incident,4,4) = 1;
758  hit_tileElemAt2(incident,5,5) = 1;
759  hit_tileElemAt2(incident,6,6) = 1;
760  hit_tileElemAt2(incident,7,7) = 1;
761 
762  // Allocate memory
763  global_shape = hit_shape(1,hit_sig(0,Nvertices-1,1));
764  hit_tileDomainShapeAlloc(&graph,Node,global_shape);
765 
766  // Set initial positions
767  Node n0 = {{{0,0,0}},1};
768  hit_tileElemAt(graph,1,0) = n0;
769  Node n1 = {{{100,0,0}},0};
770  hit_tileElemAt(graph,1,1) = n1;
771  Node n2 = {{{200,0,0}},0};
772  hit_tileElemAt(graph,1,2) = n2;
773  Node n3 = {{{300,0,0}},1};
774  hit_tileElemAt(graph,1,3) = n3;
775  Node n4 = {{{0,100,200}},1};
776  hit_tileElemAt(graph,1,4) = n4;
777  Node n5 = {{{100,100,0}},0};
778  hit_tileElemAt(graph,1,5) = n5;
779  Node n6 = {{{200,100,0}},0};
780  hit_tileElemAt(graph,1,6) = n6;
781  Node n7 = {{{300,100,0}},1};
782  hit_tileElemAt(graph,1,7) = n7;
783 
784  }
785 
786 
787 }
788 
789 
790 
795 double gradient_norm(){
796 
797  HitTile_double norm, sum_norm;
798 
799  hit_tileDomainAlloc(&norm, double, 1,1);
800  hit_tileDomainAlloc(&sum_norm, double, 1,1);
801 
802  hit_tileElemAt(norm,1,0) = 0.0;
803 
804  int vertex;
805 
806  hit_shapeIterator(vertex,shape,0){
807 
808  Node current = hit_tileElemAt(graph,1,vertex);
809 
810  if (!current.fixed){
811 
812  calculate_GH(vertex);
813 
814  Vector Gi = hit_tileElemAt(gradient,1,vertex);
815 
816  int d;
817  for(d=0;d<D;d++){
818  hit_tileElemAt(norm,1,0) += pow(Gi.v[d],2);
819  }
820 
821  //printf("vertex: %d: %.3f - %.3f - %.3f \n",vertex,pow(Gi.v[0],2),pow(Gi.v[1],2),pow(Gi.v[2],2));
822 
823  }
824  }
825 
826 
827  HitOp op_sum;
829 
830 
831  HitRanks root = {{0,0,0,-1}};
832  HitCom com_sum = hit_comReduce(lay, root, &norm, &sum_norm, HIT_DOUBLE, op_sum);
833  hit_comDo(&com_sum);
834  hit_comFree(com_sum);
835 
836  double res = 0;
837 
838  if( hit_Rank == 0 ){
839  res = sqrt(hit_tileElemAt(sum_norm,1,0));
840  }
841 
842  hit_tileFree(norm);
843  hit_tileFree(sum_norm);
844 
845 
846  return res;
847 
848 }
849 
850 
851 
852 
856 int main(int argc, char ** argv) {
857 
858  int i;
859  int vertex;
860  HitClock init_time, comp_time;
861 
862  hit_comInit(&argc,&argv);
863 
864  hit_clockStart(init_time);
865 
866  //if(hit_NProcs == 1){
867  // printf("\nERROR: Lauch with at least 2 processes\n\n");
868  // exit(EXIT_FAILURE);
869  //}
870 
871  init_graph(argc,argv);
872 
873  // Create the topology object.
874  HitTopology topo = hit_topology(plug_topPlain);
875 
876  // Create the layout to send the incident matrix
877  HitLayout lay_all = hit_layout(plug_layBlocks,topo,hit_shape(1,hit_sig(0,hit_NProcs-1,1)));
878 
879  // Root processor.
880  HitRanks root = {{0,-1,-1,-1}};
881 
882  // Send the number of vertices to allocate the tile.
883  HitTile_int vertices;
884  hit_tileDomainAlloc(&vertices, int, 1,1);
885  hit_tileElemAt(vertices,1,0) = Nvertices;
886  HitCom sendNVertices = hit_comBroadcast(lay_all,root,&vertices,HIT_INT);
887  hit_comDo(&sendNVertices);
888  Nvertices = hit_tileElemAt(vertices,1,0);
889 
890  // Allocate the tile for the incident matrix
891  if(hit_Rank != 0){
892  HitSig sig = hit_sig(0,Nvertices-1,1);
893  incident_shape = hit_shape(2,sig,sig);
894  hit_tileDomainShapeAlloc(&incident,int,incident_shape);
895  }
896 
897  // Send the incident matrix
898  HitCom sendIncident = hit_comBroadcast(lay_all,root,&incident,HIT_INT);
899  hit_comDo(&sendIncident);
900 
901 
902  // Create the global shape
903  global_shape = hit_shape(1,hit_sig(0,Nvertices-1,1));
904 
905  // Distribute the graph vertices among the processors.
906  lay = hit_layout(plug_layBlocks,topo,global_shape);
907 
908  // Get shape
909  shape = hit_layShape(lay);
910 
911  // Allocate memory of all the variables
912  init_structures();
913 
914  // Create the HitVector type.
915  HitType HitVector;
916  hit_comTypeStruct(&HitVector,Vector,1, v,3,HIT_DOUBLE);
917 
918  // Create the HitNode type
919  HitType HitNode;
920  hit_comTypeStruct(&HitNode,Node,2, r,1,HitVector,fixed,1,HIT_INT);
921 
922 
923  // Create the communication objects.
924  HitCom comA = hit_comBroadcast(lay_all,root,&graph,HitNode);
925  HitPattern patB = hit_pattern( HIT_PAT_UNORDERED );
926  HitPattern patC = hit_pattern( HIT_PAT_UNORDERED );
927 
928  // @javfres 2015-10-09 Support for more processes than nodes
929  //for(i=0;i<hit_NProcs;i++){
930  if(hit_layImActive( lay ))
931  for(i=0;i<hit_layNumActives(lay);i++){
932 
933  HitRanks ranks = {{i,-1,-1,-1}};
934  HitShape otherShape = hit_layShapeOther(lay, ranks);
935  hit_patternAdd(&patB,
936  hit_comBroadcastSelect(lay,ranks,&graph,otherShape,HIT_COM_ARRAYCOORDS,HitNode));
937  hit_patternAdd(&patC,
938  hit_comBroadcastSelect(lay,ranks,&e,otherShape,HIT_COM_ARRAYCOORDS,HitVector));
939 
940  }
941 
942 
943  hit_clockStop(init_time);
944  hit_clockStart(comp_time);
945 
946  // Send all the data from the root proc to the other procs.
947  hit_comDo(&comA);
948 
949  // Calculate the inital gradient norm.
950  //double ignorm = gradient_norm();
951  //printfRoot("# Initial gradient norm: %lf\n",ignorm);
952 
953 
954  // Main loop for the newton method
955  for(i=0;i<iter1;i++){
956  printfRoot("# iter: %d/%d\n",i,iter1);
957 
958  // Iterate trough all the vertices and
959  // obtain the gradient and the hessian.
960 
961  hit_shapeIterator(vertex,shape,0){
962 
963  Node current = hit_tileElemAt(graph,1,vertex);
964  if (!current.fixed)
965  calculate_GH(vertex);
966  }
967 
968 
969  // The initial position displacement is zero.
970  bzero(e.data,sizeof(Vector) * (size_t) Nvertices);
971 
972  // Loop for the jacobi method to solve the system.
973  int j;
974  for(j=0;j<iter2;j++){
975 
976 
977 
978  // Perform a iteration.
979  hit_shapeIterator(vertex,shape,0){
980  Node current = hit_tileElemAt(graph,1,vertex);
981  if (!current.fixed)
982  solve_system_iter(vertex);
983  }
984 
985 
986 
987  // Update the displacement.
988  memcpy(e.data,new_e.data,sizeof(Vector)* (size_t) Nvertices);
989  hit_patternDo(patC);
990 
991  }
992 
993  // Update the position with the final displacement.
994  hit_shapeIterator(vertex,shape,0){
995 
996  Node * current = &(hit_tileElemAt(graph,1,vertex));
997  if (!current->fixed){
998  int d;
999  for(d=0;d<D;d++)
1000  current->r.v[d] -= hit_tileElemAt(e,1,vertex).v[d];
1001  }
1002  }
1003 
1004  hit_patternDo(patB);
1005  }
1006 
1007 
1008  double gnorm = gradient_norm();
1009  printfRoot("# Final gradient norm: %lf\n",gnorm);
1010 
1011  hit_clockStop(comp_time);
1012  hit_clockWorldReduce(init_time);
1013  hit_clockWorldReduce(comp_time);
1014 
1015  printfRoot("# Init time: %lf\n",hit_clockGetSeconds(init_time));
1016  printfRoot("# Comp time: %lf\n",hit_clockGetSeconds(comp_time));
1017  printfRoot("# Total time: %lf\n",hit_clockGetSeconds(init_time)+hit_clockGetSeconds(comp_time));
1018 
1019  hit_comFree(comA);
1020  hit_patternFree(&patB);
1021  hit_patternFree(&patC);
1022  free_structures();
1023 
1024  hit_layFree( lay );
1025  hit_layFree( lay_all );
1026  hit_topFree( topo );
1027  hit_comFinalize();
1028 
1029  return 0;
1030 }
1031 
1032 
1036 Vector g(int i, int j){
1037 
1038  Vector gij;
1039 
1040  if(i == j){
1041  clearV(gij);
1042  return gij;
1043  }
1044 
1045  // Get the two vectors
1046  Vector r_i = (hit_tileElemAt(graph,1,i)).r;
1047  Vector r_j = (hit_tileElemAt(graph,1,j)).r;
1048 
1049  // Auxiliar variables.
1050  Vector v;
1051 
1052  // v = r_i - r_j
1053  subV(v,r_i,r_j);
1054  // n = ||v||
1055  double n = norm(v);
1056 
1057  if( n == 0){
1058  clearV(gij);
1059  return gij;
1060  }
1061 
1062  // gij = - K * ((L - n) / n) * (v)
1063  double scalar = - K * ((L - n) / n);
1064  multSV(gij,scalar,v);
1065 
1066  return gij;
1067 }
1068 
1069 
1070 
1071 
1075 Matrix W(int i, int j){
1076 
1077  Matrix Wij;
1078 
1079  if(i == j){
1080  clearM(Wij);
1081  return Wij;
1082  }
1083 
1084  // Get the two vectors
1085  Vector r_i = (hit_tileElemAt(graph,1,i)).r;
1086  Vector r_j = (hit_tileElemAt(graph,1,j)).r;
1087 
1088  // Auxiliary variables.
1089  Vector v;
1090  Matrix a, b;
1091  double scalar_a, scalar_b;
1092 
1093  // v = r_i - r_j;
1094  subV(v,r_i,r_j);
1095  // n = ||v||
1096  double n = norm(v);
1097 
1098  if( n == 0){
1099  clearM(Wij);
1100  return Wij;
1101  }
1102 
1103  // h = -((L-n)/n) * eye(D) + (L/(n^3)) * (v * v')
1104  scalar_a = -((L-n)/(n));
1105  multSI(a,scalar_a);
1106  scalar_b = (L/(pow(n,3)));
1107  multVVt(b,v,v);
1108  multSM(b,scalar_b,b);
1109  addM(Wij,a,b);
1110  // h = K * h
1111  multSM(Wij,K,Wij);
1112 
1113  return Wij;
1114 }
1115 
1116 
1121 void calculate_GH(int vertex){
1122 
1123  Vector gs;
1124  Matrix Ws;
1125 
1126  clearV(gs);
1127  clearM(Ws);
1128 
1129  int edge;
1130  for(edge=0;edge<Nvertices;edge++){
1131 
1132 
1133  if(! hit_tileElemAt(incident,2,vertex,edge)){
1134  continue;
1135  }
1136 
1137 
1138  // Get the neighbor.
1139  int nghb_index = edge;
1140 
1141  // Calculate gradient[i]
1142  Vector gij = g(vertex,nghb_index);
1143  addV(gs,gs,gij);
1144 
1145  // Calculate part of the hessian[i][i]
1146  Matrix Wij = W(vertex,nghb_index);
1147  addM(Ws,Ws,Wij);
1148 
1149 
1150  // Calculate hessian[i][j] if j is variable.
1151  Node nj = hit_tileElemAt(graph,1,nghb_index);
1152 
1153  if(!nj.fixed){
1154  Matrix Hij = W(nghb_index,vertex);
1155  multSM(Hij,-1,Hij);
1156 
1157  hit_tileElemAt(hessian,2,vertex,nghb_index) = Hij;
1158 
1159 
1160  }
1161 
1162  }
1163 
1164  hit_tileElemAt(gradient,1,vertex) = gs;
1165  hit_tileElemAt(hessian,2,vertex,vertex) = Ws;
1166 }
1167 
1168 
1169 
1175 void solve_system_iter(int vertex){
1176 
1177  Matrix Hii = hit_tileElemAt(hessian,2,vertex,vertex);
1178  Vector Gi = hit_tileElemAt(gradient,1,vertex);
1179 
1180  Vector new_ei;
1181 
1182  // Aux sum vector.
1183  Vector sum;
1184  clearV(sum);
1185 
1186  int edge;
1187 
1188  for(edge=0;edge<Nvertices;edge++){
1189 
1190 
1191  if(! hit_tileElemAt(incident,2,vertex,edge)){
1192  continue;
1193  }
1194 
1195  // Get the neighbor.
1196  int nghb_index = edge;
1197  if(nghb_index == vertex) continue;
1198 
1199  // If the node is fixed, it is no part of the ecuation system.
1200  Node nj = hit_tileElemAt(graph,1,nghb_index);
1201  if(nj.fixed) continue;
1202 
1203  // Get the j displacement and ij hessian part.
1204  Vector ej = hit_tileElemAt(e,1,nghb_index);
1205  Matrix Hij = hit_tileElemAt(hessian,2,vertex,nghb_index);
1206 
1207  // Aux par vector.
1208  Vector part;
1209 
1210  // Calculate and add the contribution of vertex j.
1211  // SUM_{i not j} ( Hij e_j )
1212  multMV(part,Hij,ej);
1213  addV(sum,sum,part);
1214 
1215  }
1216 
1217 
1218  // Gi - SUM_{i not j} ( Hij e_j )
1219  subV(sum,Gi,sum);
1220 
1221  // e_i = H_ii^(-1) * (Gi - SUM_{i not j} ( Hij e_j ) )
1222  new_ei = solve(Hii,sum);
1223 
1224  hit_tileElemAt(new_e,1,vertex) = new_ei;
1225 
1226 
1227 }
#define printfRoot(...)
Definition: spring_dense.c:58
#define D
Definition: spring_dense.c:62
Vector g(int i, int j)
#define L
Definition: spring_dense.c:67
#define hit_shape(nd,...)
Definition: hit_sshape.h:175
MPI_Op HitOp
Definition: hit_com.h:118
void print_help(char *name)
Definition: heat.c:234
#define hit_layShape(lay)
Definition: hit_layout.h:650
#define hit_tileNewType(baseType)
Definition: hit_tile.h:163
void hit_comOpSumDouble(void *, void *, int *, HitType *)
Definition: hit_com.c:2665
#define A
Definition: mg.c:64
HitOp op_sum
Definition: mg.c:176
double v[D]
Definition: spring_bitmap.c:86
#define HIT_INT
Definition: hit_com.h:74
int Nvertices
Definition: spring_dense.c:121
HitCom hit_comBroadcastSelect(HitLayout lay, HitRanks root, const void *tile, HitShape selection, int modeSelect, HitType baseType)
Definition: hit_com.c:731
#define hit_Rank
Definition: hit_com.h:140
int iter2
Definition: spring_bitmap.c:80
int iter
Definition: mmult_bit.c:68
#define addM(r, a, b)
Definition: spring_dense.c:302
#define hit_tileElemAt(var, ndims,...)
Definition: hit_tile.h:519
#define HIT_LAYOUT_NULL_STATIC
Definition: hit_layout.h:300
#define multSV(r, s, vv)
Definition: spring_dense.c:239
void hit_comFree(HitCom issue)
Definition: hit_com.c:1995
#define norm(vv)
Definition: spring_dense.c:232
char * replace_str(char *str, const char *orig, const char *rep)
Matrix W(int i, int j)
void hit_comDo(HitCom *issue)
Definition: hit_com.c:2408
HitTile_Matrix hessian
#define hit_tileDomainAlloc(newVarP, baseType, numDims,...)
Definition: hit_tile.h:336
void solve_system_iter(int vertex)
#define hit_topology(name,...)
Definition: hit_topology.h:308
#define hit_layShapeOther(lay, ranks)
Definition: hit_layout.h:732
#define hit_sigCard(sig)
Definition: hit_sig.h:162
#define det(mm)
Definition: spring_dense.c:262
Definition: hit_sig.h:79
#define hit_fileHBVertices(hbfile)
Definition: hit_file.h:157
#define cpyV(r, a)
Definition: spring_dense.c:210
#define multVVt(r, a, b)
Definition: spring_dense.c:312
HitTile_Vector graph
#define K
Definition: spring_dense.c:65
#define multSI(r, s)
Definition: spring_dense.c:292
#define r(a, b, c)
#define HIT_COM_ARRAYCOORDS
Definition: hit_com.h:345
#define hit_fileHBReadDense(hbfile, tileP)
Definition: hit_file.h:183
#define HIT_PAT_UNORDERED
Definition: hit_pattern.h:81
void hit_patternAdd(HitPattern *pattern, HitCom comm)
Definition: hit_pattern.c:61
Vector solve(Matrix A, Vector b)
void free_structures()
#define hit_NProcs
Definition: hit_com.h:142
#define P_FIXED
Definition: spring_dense.c:70
#define hit_tileDomainShapeAlloc(newVarP, baseType, shape)
Definition: hit_tile.h:354
HitTile_int fixed
#define x
int fixed
Definition: spring_dense.c:101
#define subV(r, a, b)
Definition: spring_dense.c:225
#define addV(r, a, b)
Definition: spring_dense.c:217
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)
Definition: spring_dense.c:322
Vector solve_jacobi(Matrix A, Vector b)
#define hit_clockWorldReduce(c)
Definition: hit_utils.h:156
void hit_patternFree(HitPattern *pattern)
Definition: hit_pattern.c:94
void hit_comInit(int *pargc, char **pargv[])
Definition: hit_com.c:111
HitShape shape
#define hit_layImActive(lay)
Definition: hit_layout.h:797
HitShape incident_shape
Definition: spring_dense.c:128
void random_coordinates()
#define clearM(mm)
Definition: spring_dense.c:246
void hit_layFree(HitLayout lay)
Definition: hit_layout.c:2251
#define hit_clockStop(c)
Definition: hit_utils.h:109
#define clearV(vv)
Definition: spring_dense.c:203
HitShape global_shape
int main(int argc, char *argv[])
Definition: cannonAsync.c:62
#define inv(r, mm)
Definition: spring_dense.c:332
#define m(a, b, c)
HitTile_Vector new_e
Vector r
Definition: spring_dense.c:100
#define hit_comReduce(lay, root, tilePSend, tilePRecv, baseType, operation)
Definition: hit_com.h:725
#define hit_comTypeStruct(new_type, Nstruct, n,...)
Definition: hit_com.h:208
void init_graph(HitTile_double graph, HitShape shape_global)
Definition: heat.c:244
HitTile_Vector gradient
void init_structures()
HitTile_Vector e
#define hit_shapeIterator(var, shape, dim)
Definition: hit_sshape.h:603
#define ITER2
Definition: spring_dense.c:75
int hit_layNumActives(HitLayout lay)
Definition: hit_layout.c:2213
void calculate_GH(int vertex)
#define v(a, b, c)
#define hit_comBroadcast(lay, root, tile, baseType)
Definition: hit_com.h:673
#define multSM(r, s, mm)
Definition: spring_dense.c:283
HitTile_int incident
Definition: spring_dense.c:122
#define hit_layout(name, topo,...)
Definition: hit_layout.h:415
double gradient_norm()
HitLayout lay
Definition: heat.c:107
#define hit_shapeSig(shape, dim)
Definition: hit_sshape.h:400
#define ITER_SJ
#define hit_tileFree(var)
Definition: hit_tile.h:369
#define ITER1
Definition: spring_dense.c:73
#define HIT_DOUBLE
Definition: hit_com.h:78
#define hit_patternDo(pattern)
Definition: hit_pattern.h:157
void hit_comFinalize()
Definition: hit_com.c:159
#define hit_clockGetSeconds(c)
Definition: hit_utils.h:190
int iter1
Definition: spring_bitmap.c:78
double ** buffer
Definition: refMPIluBack.c:103
#define hit_clockStart(c)
Definition: hit_utils.h:87
double m[D][D]
Definition: spring_bitmap.c:93