Hitmap 1.3
 All Data Structures Namespaces Files Functions Variables Typedefs Friends Macros Groups Pages
spring_bitmap_sec.c
Go to the documentation of this file.
1 
17 /*
18  * <license>
19  *
20  * Hitmap v1.2
21  *
22  * This software is provided to enhance knowledge and encourage progress in the scientific
23  * community. It should be used only for research and educational purposes. Any reproduction
24  * or use for commercial purpose, public redistribution, in source or binary forms, with or
25  * without modifications, is NOT ALLOWED without the previous authorization of the copyright
26  * holder. The origin of this software must not be misrepresented; you must not claim that you
27  * wrote the original software. If you use this software for any purpose (e.g. publication),
28  * a reference to the software package and the authors must be included.
29  *
30  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS "AS IS" AND ANY
31  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
32  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
33  * THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
34  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
35  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
36  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
37  * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
38  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
39  *
40  * Copyright (c) 2007-2015, Trasgo Group, Universidad de Valladolid.
41  * All rights reserved.
42  *
43  * More information on http://trasgo.infor.uva.es/
44  *
45  * </license>
46 */
47 
48 #include <stdio.h>
49 #include <stdlib.h>
50 #include <hitmap.h>
51 #include <hit_com.h>
52 #include <unistd.h>
53 #include <string.h>
54 
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 
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 
103 HitTile_Vector gradient;
104 
106 HitTile_Matrix hessian;
107 
109 HitTile_Vector e;
110 
112 HitTile_Vector new_e;
113 
116 
118 HitTile_Vector graph;
119 
121 HitTile_int fixed;
122 
123 /*
124  * Function to calculate the gradient and hessian for a given node.
125  * @param vertex The vertex index.
126  */
127 void calculate_GH(int vertex);
128 
133 void solve_system_iter(int vertex);
134 
135 
136 #ifdef DEBUG
137 
141 void print_matrix(Matrix m){
142 
143  int i,j;
144  for(i=0;i<D;i++){
145  for(j=0;j<D;j++){
146  printf("%.2lf\t",m.m[i][j]);
147  }
148  printf("\n");
149  }
150 }
151 
156 void print_vector(Vector v){
157 
158  int i;
159  for(i=0;i<D;i++){
160  printf("%.2lf\t",v.v[i]);
161  }
162  printf("\n");
163 }
164 #endif
165 
169 #define clearV(vv) clearVInternal((vv).v)
170 #define clearVInternal(v) \
171  ((v)[0] = (v)[1] = (v)[2] = 0.0 )
172 
176 #define cpyV(r,a) cpyVInternal((r).v,(a).v)
177 #define cpyVInternal(r,a) \
178  {(r)[0] = (a)[0]; (r)[1] = (a)[1]; (r)[2] = (a)[2]; }
179 
183 #define addV(r,a,b) addVInternal((r).v,(a).v,(b).v)
184 #define addVInternal(r,a,b) \
185  {(r)[0] = (a)[0] + (b)[0]; (r)[1] = (a)[1] + (b)[1]; (r)[2] = (a)[2] + (b)[2]; }
186 
187 
191 #define subV(r,a,b) subVInternal((r).v,(a).v,(b).v)
192 #define subVInternal(r,a,b) \
193  {(r)[0] = (a)[0] - (b)[0]; (r)[1] = (a)[1] - (b)[1]; (r)[2] = (a)[2] - (b)[2]; }
194 
198 #define norm(vv) normInternal((vv).v)
199 #define normInternal(v) \
200  (sqrt( pow((v)[0],2) + pow((v)[1],2) + pow((v)[2],2) ))
201 
205 #define multSV(r,s,vv) multSVInternal((r).v,s,(vv).v)
206 #define multSVInternal(r,s,v) \
207  {(r)[0] = (s) * (v)[0]; (r)[1] = (s) * (v)[1]; (r)[2] = (s) * (v)[2]; }
208 
212 #define clearM(mm) clearMInternal((mm).m)
213 #define clearMInternal(m) \
214  ((m)[0][0] = (m)[0][1] = (m)[0][2] = \
215  (m)[1][0] = (m)[1][1] = (m)[1][2] = \
216  (m)[2][0] = (m)[2][1] = (m)[2][2] = 0.0 )
217 
218 
228 #define det(mm) detInternal((mm).m)
229 #define detInternal(m) \
230  ( (m)[0][0] * (m)[1][1] * (m)[2][2] \
231  + (m)[0][1] * (m)[1][2] * (m)[2][0] \
232  + (m)[0][2] * (m)[1][0] * (m)[2][1] \
233  - (m)[0][0] * (m)[1][2] * (m)[2][1] \
234  - (m)[0][1] * (m)[1][0] * (m)[2][2] \
235  - (m)[0][2] * (m)[1][1] * (m)[2][0] )
236 
240 #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)
241 #define init_matrixInternal(m,a,b,c,d,e,f,g,h,i) \
242  {(m)[0][0] = a; (m)[0][1] = b; (m)[0][2] = c; \
243  (m)[1][0] = d; (m)[1][1] = e; (m)[1][2] = f; \
244  (m)[2][0] = g; (m)[2][1] = h; (m)[2][2] = i; }
245 
249 #define multSM(r,s,mm) multSMInternal((r).m,s,(mm).m)
250 #define multSMInternal(r,s,m) \
251  {(r)[0][0] = (s) * (m)[0][0]; (r)[0][1] = (s) * (m)[0][1]; (r)[0][2] = (s) * (m)[0][2]; \
252  (r)[1][0] = (s) * (m)[1][0]; (r)[1][1] = (s) * (m)[1][1]; (r)[1][2] = (s) * (m)[1][2]; \
253  (r)[2][0] = (s) * (m)[2][0]; (r)[2][1] = (s) * (m)[2][1]; (r)[2][2] = (s) * (m)[2][2]; }
254 
258 #define multSI(r,s) multSIInternal((r).m,s)
259 #define multSIInternal(r,s) \
260  {(r)[0][0] = (s); (r)[0][1] = 0; (r)[0][2] = 0; \
261  (r)[1][0] = 0; (r)[1][1] = (s); (r)[1][2] = 0; \
262  (r)[2][0] = 0; (r)[2][1] = 0; (r)[2][2] = (s); }
263 
264 
268 #define addM(r,a,b) addMInternal((r).m,(a).m,(b).m)
269 #define addMInternal(r,a,b) \
270  {(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]; \
271  (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]; \
272  (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]; }
273 
278 #define multVVt(r,a,b) multVVtInternal((r).m,(a).v,(b).v)
279 #define multVVtInternal(r,a,b) \
280  {(r)[0][0] = (a)[0] * (b)[0]; (r)[0][1] = (a)[0] * (b)[1]; (r)[0][2] = (a)[0] * (b)[2]; \
281  (r)[1][0] = (a)[1] * (b)[0]; (r)[1][1] = (a)[1] * (b)[1]; (r)[1][2] = (a)[1] * (b)[2]; \
282  (r)[2][0] = (a)[2] * (b)[0]; (r)[2][1] = (a)[2] * (b)[1]; (r)[2][2] = (a)[2] * (b)[2]; }
283 
284 
288 #define multMV(r,mm,vv) multMVInternal((r).v,(mm).m,(vv).v)
289 #define multMVInternal(r,m,v) \
290  {(r)[0] = (m)[0][0] * (v)[0] + (m)[0][1] * (v)[1] + (m)[0][2] * (v)[2]; \
291  (r)[1] = (m)[1][0] * (v)[0] + (m)[1][1] * (v)[1] + (m)[1][2] * (v)[2]; \
292  (r)[2] = (m)[2][0] * (v)[0] + (m)[2][1] * (v)[1] + (m)[2][2] * (v)[2]; }
293 
294 
298 #define inv(r,mm) invInternal((r).m,(mm).m)
299 #define invInternal(r,m) \
300  {(r)[0][0] = (m)[1][1] * (m)[2][2] - (m)[1][2] * (m)[2][1]; \
301  (r)[0][1] = - (m)[0][1] * (m)[2][2] + (m)[0][2] * (m)[2][1]; \
302  (r)[0][2] = (m)[0][1] * (m)[1][2] - (m)[0][2] * (m)[1][1]; \
303  (r)[1][0] = - (m)[1][0] * (m)[2][2] + (m)[1][2] * (m)[2][0]; \
304  (r)[1][1] = (m)[0][0] * (m)[2][2] - (m)[0][2] * (m)[2][0]; \
305  (r)[1][2] = - (m)[0][0] * (m)[1][2] + (m)[0][2] * (m)[1][0]; \
306  (r)[2][0] = (m)[1][0] * (m)[2][1] - (m)[1][1] * (m)[2][0]; \
307  (r)[2][1] = - (m)[0][0] * (m)[2][1] + (m)[0][1] * (m)[2][0]; \
308  (r)[2][2] = (m)[0][0] * (m)[1][1] - (m)[0][1] * (m)[1][0]; \
309  multSMInternal(r,1/detInternal(m),r); }
310 
314 #define multMM(r,a,b) multMMInternal((r).m,(a).m,(b).m)
315 #define multMMInternal(r,a,b) \
316  {(r)[0][0] = (a)[0][0] * (b)[0][0] + (a)[0][1] * (b)[1][0] + (a)[0][2] * (b)[2][0]; \
317  (r)[0][1] = (a)[0][0] * (b)[0][1] + (a)[0][1] * (b)[1][1] + (a)[0][2] * (b)[2][1]; \
318  (r)[0][2] = (a)[0][0] * (b)[0][2] + (a)[0][1] * (b)[1][2] + (a)[0][2] * (b)[2][2]; \
319  (r)[1][0] = (a)[1][0] * (b)[0][0] + (a)[1][1] * (b)[1][0] + (a)[1][2] * (b)[2][0]; \
320  (r)[1][1] = (a)[1][0] * (b)[0][1] + (a)[1][1] * (b)[1][1] + (a)[1][2] * (b)[2][1]; \
321  (r)[1][2] = (a)[1][0] * (b)[0][2] + (a)[1][1] * (b)[1][2] + (a)[1][2] * (b)[2][2]; \
322  (r)[2][0] = (a)[2][0] * (b)[0][0] + (a)[2][1] * (b)[1][0] + (a)[2][2] * (b)[2][0]; \
323  (r)[2][1] = (a)[2][0] * (b)[0][1] + (a)[2][1] * (b)[1][1] + (a)[2][2] * (b)[2][1]; \
324  (r)[2][2] = (a)[2][0] * (b)[0][2] + (a)[2][1] * (b)[1][2] + (a)[2][2] * (b)[2][2]; }
325 
326 
335 
336 // Number of iteration of the jacobi iterative method.
337 #define ITER_SJ 10
338 
339  Vector x;
340 
341  // The new approximation.
342  Vector x_1;
343 
344  // Clear the vectors.
345  clearV(x);
346  clearV(x_1);
347 
348  int iter;
349  int i,j;
350 
351  // Iterative method.
352  for(iter=0; iter<ITER_SJ; iter++){
353  for(i=0;i<D;i++){
354 
355  if( A.m[i][i] == 0 ) continue;
356 
357  // calculate the new value.
358  x_1.v[i] = 0.0;
359  for(j=0;j<D;j++){
360  if(i != j) x_1.v[i] += A.m[i][j] * x.v[j];
361  }
362  x_1.v[i] = (b.v[i] - x_1.v[i]) / A.m[i][i];
363 
364  }
365  cpyV(x,x_1);
366  }
367 
368  return x;
369 }
370 
371 
372 
384 
385  Vector x;
386 
387  double detA = det(A);
388 
389  // If the matrix has a determinant, we can get the inverse.
390  if(detA != 0){
391 
392  Matrix InvA;
393  inv(InvA,A);
394  multMV(x,InvA,b);
395  } else {
396  //printf("#Warning: singular matrix\n");
397  x = solve_jacobi(A,b);
398  }
399  return x;
400 }
401 
406 
411 }
412 
417 
420  hit_tileFree(e);
424 
425  hit_shapeFree(shape);
426 }
427 
428 
433 void print_help(char * name){
434  printf("%s [-n NEWTON METHOD ITERATIONS] [-j JACOBI METHOD ITERATIONS] FILE \n",name);
435 
436  printf(" -n NEWTON METHOD ITERATIONS number of iterations (default 10)\n");
437  printf(" -j JACOBI METHOD ITERATIONS number of iterations (default 100)\n");
438  printf(" FILE input graph file\n");
439 
440 }
441 
442 
443 
449 
450  srandom(0);
451  if(hit_Rank != 0) return;
452 
453  int n = hit_bShapeNvertices(shape);
454 
455  int j,i;
456 
457  // Set the coordinates
458  for(j=0;j<D;j++){
459  for(i=0;i<n;i++){
460  double number = (double) random();
461  hit_gbTileVertexAt(graph,i).v[j] = number;
462  }
463  }
464 
465  srandom(0);
466  for(i=0;i<n;i++){
467 
468  if ((random() % 100) < P_FIXED){
469  hit_gbTileVertexAt(fixed,i) = 1;
470  } else {
471  hit_gbTileVertexAt(fixed,i) = 0;
472  }
473  }
474 
475 
476 }
477 
478 
479 
480 
481 
482 
486 void init_graph(int argc, char ** argv){
487 
488  // Define the spring structure.
489  shape = HIT_BITMAP_SHAPE_NULL;
490 
491  // Argument.
492  char c;
493 
494  // Parse the command-line arguments.
495  while ((c = (char) getopt (argc, argv, "hn:j:")) != -1)
496  switch (c) {
497  case 'n':
498  sscanf(optarg,"%d",&iter1);
499  break;
500  case 'j':
501  sscanf(optarg,"%d",&iter2);
502  break;
503  case 'h':
504  print_help(argv[0]);
505  exit(EXIT_SUCCESS);
506  default:
507  abort ();
508  break;
509  }
510 
511 
512  if(optind == (argc-1)){
513 
514  // Graph file name.
515  char * graph_file = argv[optind];
516 
517  printf("# Graph: %s\n",graph_file);
518  printf("# iterations %d x %d\n",iter1,iter2);
519 
520  shape = hit_fileHBReadBitmap(graph_file);
521 
522  printf("# Vertices %d\n",hit_bShapeNvertices(shape));
523 
524  // Allocate memory
527 
528 
529  // Random coordinates
531 
532  int i;
533  int nfixed = 0;
534  for(i=0;i<hit_bShapeNvertices(shape);i++){
535  nfixed += hit_gbTileVertexAt(fixed,i);
536  }
537 
538  printf("# N fixed: %d\n",nfixed);
539 
540 
541 
542  } else {
543 
544  iter1 = iter2 = 3;
545 
546  hit_bShapeAddEdge2(&shape,0,1);
547  hit_bShapeAddEdge2(&shape,1,2);
548  hit_bShapeAddEdge2(&shape,2,3);
549  hit_bShapeAddEdge2(&shape,4,5);
550  hit_bShapeAddEdge2(&shape,5,6);
551  hit_bShapeAddEdge2(&shape,6,7);
552  hit_bShapeAddEdge2(&shape,0,4);
553  hit_bShapeAddEdge2(&shape,1,5);
554  hit_bShapeAddEdge2(&shape,2,6);
555  hit_bShapeAddEdge2(&shape,3,7);
556  hit_bShapeAddEdge2(&shape,1,4);
557  hit_bShapeAddEdge2(&shape,2,5);
558 
559  // Include each (i,i) edge.
560  hit_bShapeAddEdge2(&shape,0,0);
561  hit_bShapeAddEdge2(&shape,1,1);
562  hit_bShapeAddEdge2(&shape,2,2);
563  hit_bShapeAddEdge2(&shape,3,3);
564  hit_bShapeAddEdge2(&shape,4,4);
565  hit_bShapeAddEdge2(&shape,5,5);
566  hit_bShapeAddEdge2(&shape,6,6);
567  hit_bShapeAddEdge2(&shape,7,7);
568 
569  // Allocate memory
572 
573  // Set initial positions
574  Vector v0 = {{0,0,0}};
575  hit_gbTileVertexAt(graph,0) = v0;
576  Vector v1 = {{100,0,0}};
577  hit_gbTileVertexAt(graph,1) = v1;
578  Vector v2 = {{200,0,0}};
579  hit_gbTileVertexAt(graph,2) = v2;
580  Vector v3 = {{300,0,0}};
581  hit_gbTileVertexAt(graph,3) = v3;
582  Vector v4 = {{0,100,200}};
583  hit_gbTileVertexAt(graph,4) = v4;
584  Vector v5 = {{100,100,0}};
585  hit_gbTileVertexAt(graph,5) = v5;
586  Vector v6 = {{200,100,0}};
587  hit_gbTileVertexAt(graph,6) = v6;
588  Vector v7 = {{300,100,0}};
589  hit_gbTileVertexAt(graph,7) = v7;
590 
591  // Set the fix/free status
592  hit_gbTileVertexAt(fixed,0) = 1;
593  hit_gbTileVertexAt(fixed,1) = 0;
594  hit_gbTileVertexAt(fixed,2) = 0;
595  hit_gbTileVertexAt(fixed,3) = 1;
596  hit_gbTileVertexAt(fixed,4) = 1;
597  hit_gbTileVertexAt(fixed,5) = 0;
598  hit_gbTileVertexAt(fixed,6) = 0;
599  hit_gbTileVertexAt(fixed,7) = 1;
600 
601  }
602 
603 }
604 
609 double gradient_norm(){
610 
611  double norm = 0.0;
612 
613  int vertex;
614 
615  hit_bShapeVertexIterator(vertex,shape) {
616  if (!hit_gbTileVertexAt(fixed,vertex)){
617 
618  calculate_GH(vertex);
619 
620  Vector Gi = hit_gbTileVertexAt(gradient,vertex);
621 
622  int d;
623  for(d=0;d<D;d++){
624  norm += pow(Gi.v[d],2);
625  }
626  }
627  }
628 
629  return sqrt(norm);
630 }
631 
632 
633 
634 
638 int main(int argc, char ** argv) {
639 
640  int i;
641  int vertex;
642  HitClock init_time, comp_time;
643 
644  hit_comInit(&argc,&argv);
645 
646  hit_clockStart(init_time);
647  init_graph(argc,argv);
648 
649  // Allocate memory of all the variables
650  init_structures();
651 
652  // printf("# Initial gradient norm: %lf\n",gradient_norm());
653  hit_clockStop(init_time);
654  hit_clockStart(comp_time);
655 
656  // Main loop for the newton method
657  for(i=0;i<iter1;i++){
658  printf("# iter: %d/%d\n",i,iter1);
659 
660  // Iterate trough all the vertices and
661  // obtain the gradient and the hessian.
662  hit_bShapeVertexIterator(vertex,shape) {
663  if (!hit_gbTileVertexAt(fixed,vertex))
664  calculate_GH(vertex);
665  }
666 
667 
668  // The initial position displacement is zero.
670 
671  // Loop for the jacobi method to solve the system.
672  int j;
673  for(j=0;j<iter2;j++){
674 
675  // Perform a iteration.
676  hit_bShapeVertexIterator(vertex,shape) {
677  if (!hit_gbTileVertexAt(fixed,vertex))
678  solve_system_iter(vertex);
679  }
680 
681  // Update the displacement.
683 
684  }
685 
686  // Update the position with the final displacement.
687  hit_bShapeVertexIterator(vertex,shape) {
688  if (!hit_gbTileVertexAt(fixed,vertex)){
689  Vector * current = &(hit_gbTileVertexAt(graph,vertex));
690  subV(*current,*current,hit_gbTileVertexAt(e,vertex));
691  }
692  }
693 
694  }
695 
696  hit_clockStop(comp_time);
697 
698  printf("# Final gradient norm: %lf\n",gradient_norm());
699  printf("# Init time: %lf\n",hit_clockGetSeconds(init_time));
700  printf("# Comp time: %lf\n",hit_clockGetSeconds(comp_time));
701  printf("# Total time: %lf\n",hit_clockGetSeconds(init_time)+hit_clockGetSeconds(comp_time));
702 
703  free_structures();
704  hit_comFinalize();
705 
706  return 0;
707 }
708 
709 
710 
714 Vector g(int i, int j){
715 
716  Vector gij;
717 
718  if(i == j){
719  clearV(gij);
720  return gij;
721  }
722 
723  // Get the two vectors
724  Vector r_i = (hit_gbTileVertexAt(graph,i));
725  Vector r_j = (hit_gbTileVertexAt(graph,j));
726 
727  // Auxiliar variables.
728  Vector v;
729 
730  // v = r_i - r_j
731  subV(v,r_i,r_j);
732  // n = ||v||
733  double n = norm(v);
734 
735  if( n == 0){
736  clearV(gij);
737  return gij;
738  }
739 
740  // gij = - K * ((L - n) / n) * (v)
741  double scalar = - K * ((L - n) / n);
742  multSV(gij,scalar,v);
743 
744  return gij;
745 }
746 
747 
751 Matrix W(int i, int j){
752 
753  Matrix Wij;
754 
755  if(i == j){
756  clearM(Wij);
757  return Wij;
758  }
759 
760  // Get the two vectors
761  Vector r_i = (hit_gbTileVertexAt(graph,i));
762  Vector r_j = (hit_gbTileVertexAt(graph,j));
763 
764  // Auxiliary variables.
765  Vector v;
766  Matrix a, b;
767  double scalar_a, scalar_b;
768 
769  // v = r_i - r_j;
770  subV(v,r_i,r_j);
771  // n = ||v||
772  double n = norm(v);
773 
774  if( n == 0){
775  clearM(Wij);
776  return Wij;
777  }
778 
779  // h = -((L-n)/n) * eye(D) + (L/(n^3)) * (v * v')
780  scalar_a = -((L-n)/(n));
781  multSI(a,scalar_a);
782  scalar_b = (L/(pow(n,3)));
783  multVVt(b,v,v);
784  multSM(b,scalar_b,b);
785  addM(Wij,a,b);
786  // h = K * h
787  multSM(Wij,K,Wij);
788 
789  return Wij;
790 }
791 
792 
798 void calculate_GH(int vertex){
799 
800  Vector gs;
801  Matrix Ws;
802 
803  clearV(gs);
804  clearM(Ws);
805 
806  int edge;
807  hit_bShapeEdgeIterator(edge,shape,vertex){
808 
809  // Get the neighbor.
810  int nghb_index = hit_bShapeEdgeTarget(shape,edge);
811 
812  // Calculate gradient[i]
813  Vector gij = g(vertex,nghb_index);
814  addV(gs,gs,gij);
815 
816  // Calculate part of the hessian[i][i]
817  Matrix Wij = W(vertex,nghb_index);
818  addM(Ws,Ws,Wij);
819 
820  // Calculate hessian[i][j] if j is variable.
821  //Node nj = hit_gbTileVertexAt(graph,nghb_index);
822  if(!hit_gbTileVertexAt(fixed,nghb_index)){
823 
824  Matrix Hij = W(nghb_index,vertex);
825  multSM(Hij,-1,Hij);
826 
827  hit_gbTileEdgeIteratorAt(hessian,vertex,edge) = Hij;
828  }
829  }
830 
831  hit_gbTileVertexAt(gradient,vertex) = gs;
832  hit_gbTileEdgeAt(hessian,vertex,vertex) = Ws;
833 
834 }
835 
841 void solve_system_iter(int vertex){
842 
843  Matrix Hii = hit_gbTileEdgeAt(hessian,vertex,vertex);
844  Vector Gi = hit_gbTileVertexAt(gradient,vertex);
845 
846  Vector new_ei;
847 
848  // Aux sum vector.
849  Vector sum;
850  clearV(sum);
851 
852  int edge;
853  hit_bShapeEdgeIterator(edge,shape,vertex){
854 
855  // Get the neighbor.
856  int nghb_index = hit_bShapeEdgeTarget(shape,edge);
857  if(nghb_index == vertex) continue;
858 
859  // If the node is fixed, it is no part of the equation system.
860  //Node nj = hit_gbTileVertexAt(graph,nghb_index);
861  if(hit_gbTileVertexAt(fixed,nghb_index)) continue;
862 
863  // Get the j displacement and ij hessian part.
864  Vector ej = hit_gbTileVertexAt(e,nghb_index);
865  Matrix Hij = hit_gbTileEdgeIteratorAt(hessian,vertex,edge);
866 
867  // Aux par vector.
868  Vector part;
869 
870  // Calculate and add the contribution of vertex j.
871  // SUM_{i not j} ( Hij e_j )
872  multMV(part,Hij,ej);
873  addV(sum,sum,part);
874 
875  }
876 
877  // Gi - SUM_{i not j} ( Hij e_j )
878  subV(sum,Gi,sum);
879 
880  // e_i = H_ii^(-1) * (Gi - SUM_{i not j} ( Hij e_j ) )
881  new_ei = solve(Hii,sum);
882  hit_gbTileVertexAt(new_e,vertex) = new_ei;
883 }
884 
885 
#define addM(r, a, b)
Vector g(int i, int j)
#define multSM(r, s, mm)
void print_help(char *name)
Definition: heat.c:234
#define hit_tileNewType(baseType)
Definition: hit_tile.h:163
#define hit_bShapeAddEdge2(shape, x, y)
Definition: hit_bshape.h:432
#define A
Definition: mg.c:64
double v[D]
Definition: spring_bitmap.c:86
#define ITER2
#define P_FIXED
#define hit_fileHBReadBitmap(hbfile)
Definition: hit_file.h:138
#define hit_Rank
Definition: hit_com.h:140
int iter2
Definition: spring_bitmap.c:80
int iter
Definition: mmult_bit.c:68
#define hit_bShapeEdgeTarget(s, edge)
#define hit_bShapeEdgeIterator(var, shape, vertex)
Definition: hit_bshape.h:494
#define ITER_SJ
#define D
Matrix W(int i, int j)
void hit_gbTileCopyVertices(void *destP, void *srcP)
Definition: hit_gbtile.c:159
void hit_gbTileClearVertices(void *varP)
Definition: hit_gbtile.c:146
HitTile_Matrix hessian
void solve_system_iter(int vertex)
#define norm(vv)
HitTile_Vector graph
#define HIT_EDGES
Definition: hit_tile.h:214
#define cpyV(r, a)
#define hit_gbTileVertexAt(var, vertex)
Definition: hit_gbtile.h:154
Vector solve(Matrix A, Vector b)
#define multSV(r, s, vv)
void free_structures()
#define hit_bShapeVertexIterator(var, shape)
Definition: hit_bshape.h:455
#define ITER1
HitTile_int fixed
#define clearV(vv)
#define x
#define L
#define det(mm)
Vector solve_jacobi(Matrix A, Vector b)
#define multSI(r, s)
void hit_comInit(int *pargc, char **pargv[])
Definition: hit_com.c:111
HitShape shape
#define K
#define subV(r, a, b)
#define multVVt(r, a, b)
#define hit_bShapeNvertices(shape)
#define inv(r, mm)
void random_coordinates()
#define hit_gbTileEdgeAt(var, pos1, pos2)
Definition: hit_gbtile.h:164
#define hit_clockStop(c)
Definition: hit_utils.h:109
int main(int argc, char *argv[])
Definition: cannonAsync.c:62
#define addV(r, a, b)
#define m(a, b, c)
HitTile_Vector new_e
HitShape HIT_BITMAP_SHAPE_NULL
Definition: hit_shape.c:71
void init_graph(HitTile_double graph, HitShape shape_global)
Definition: heat.c:244
HitTile_Vector gradient
void init_structures()
#define clearM(mm)
HitTile_Vector e
void calculate_GH(int vertex)
#define v(a, b, c)
double gradient_norm()
#define multMV(r, mm, vv)
#define hit_tileFree(var)
Definition: hit_tile.h:369
void hit_shapeFree(HitShape s)
Definition: hit_shape.c:84
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
#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