Hitmap 1.3
 All Data Structures Namespaces Files Functions Variables Typedefs Friends Macros Groups Pages
spring_csr_sec.c
Go to the documentation of this file.
1 
13 /*
14  * <license>
15  *
16  * Hitmap v1.2
17  *
18  * This software is provided to enhance knowledge and encourage progress in the scientific
19  * community. It should be used only for research and educational purposes. Any reproduction
20  * or use for commercial purpose, public redistribution, in source or binary forms, with or
21  * without modifications, is NOT ALLOWED without the previous authorization of the copyright
22  * holder. The origin of this software must not be misrepresented; you must not claim that you
23  * wrote the original software. If you use this software for any purpose (e.g. publication),
24  * a reference to the software package and the authors must be included.
25  *
26  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS "AS IS" AND ANY
27  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
28  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
29  * THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
31  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
32  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
33  * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
34  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
35  *
36  * Copyright (c) 2007-2015, Trasgo Group, Universidad de Valladolid.
37  * All rights reserved.
38  *
39  * More information on http://trasgo.infor.uva.es/
40  *
41  * </license>
42 */
43 
44 #include <stdio.h>
45 #include <stdlib.h>
46 #include <hitmap.h>
47 #include <hit_com.h>
48 #include <unistd.h>
49 #include <string.h>
50 
51 
53 #define D (3)
54 
56 #define K (1)
57 
58 #define L (0.001)
59 
61 #define P_FIXED 10
62 
64 #define ITER1 10
65 
66 #define ITER2 100
67 
69 int iter1 = ITER1;
71 int iter2 = ITER2;
72 
73 
77 typedef struct {
78  double v[D];
79 } Vector;
80 
84 typedef struct {
85  double m[D][D];
86 } Matrix;
87 
88 
89 // Define the tile_Vector type.
91 
92 // Define the tile_Matrix type.
94 
95 // Define the tile_int type.
96 hit_tileNewType(int);
97 
99 HitTile_Vector gradient;
100 
102 HitTile_Matrix hessian;
103 
105 HitTile_Vector e;
106 
108 HitTile_Vector new_e;
109 
112 
114 HitTile_Vector graph;
115 
117 HitTile_int fixed;
118 
119 /*
120  * Function to calculate the gradient and hessian for a given node.
121  * @param vertex The vertex index.
122  */
123 void calculate_GH(int vertex);
124 
129 void solve_system_iter(int vertex);
130 
131 
132 #ifdef DEBUG
133 
137 void print_matrix(Matrix m){
138 
139  int i,j;
140  for(i=0;i<D;i++){
141  for(j=0;j<D;j++){
142  printf("%.2lf\t",m.m[i][j]);
143  }
144  printf("\n");
145  }
146 }
147 
152 void print_vector(Vector v){
153 
154  int i;
155  for(i=0;i<D;i++){
156  printf("%.2lf\t",v.v[i]);
157  }
158  printf("\n");
159 }
160 #endif
161 
165 #define clearV(vv) clearVInternal((vv).v)
166 #define clearVInternal(v) \
167  ((v)[0] = (v)[1] = (v)[2] = 0.0 )
168 
172 #define cpyV(r,a) cpyVInternal((r).v,(a).v)
173 #define cpyVInternal(r,a) \
174  {(r)[0] = (a)[0]; (r)[1] = (a)[1]; (r)[2] = (a)[2]; }
175 
179 #define addV(r,a,b) addVInternal((r).v,(a).v,(b).v)
180 #define addVInternal(r,a,b) \
181  {(r)[0] = (a)[0] + (b)[0]; (r)[1] = (a)[1] + (b)[1]; (r)[2] = (a)[2] + (b)[2]; }
182 
183 
187 #define subV(r,a,b) subVInternal((r).v,(a).v,(b).v)
188 #define subVInternal(r,a,b) \
189  {(r)[0] = (a)[0] - (b)[0]; (r)[1] = (a)[1] - (b)[1]; (r)[2] = (a)[2] - (b)[2]; }
190 
194 #define norm(vv) normInternal((vv).v)
195 #define normInternal(v) \
196  (sqrt( pow((v)[0],2) + pow((v)[1],2) + pow((v)[2],2) ))
197 
201 #define multSV(r,s,vv) multSVInternal((r).v,s,(vv).v)
202 #define multSVInternal(r,s,v) \
203  {(r)[0] = (s) * (v)[0]; (r)[1] = (s) * (v)[1]; (r)[2] = (s) * (v)[2]; }
204 
208 #define clearM(mm) clearMInternal((mm).m)
209 #define clearMInternal(m) \
210  ((m)[0][0] = (m)[0][1] = (m)[0][2] = \
211  (m)[1][0] = (m)[1][1] = (m)[1][2] = \
212  (m)[2][0] = (m)[2][1] = (m)[2][2] = 0.0 )
213 
214 
224 #define det(mm) detInternal((mm).m)
225 #define detInternal(m) \
226  ( (m)[0][0] * (m)[1][1] * (m)[2][2] \
227  + (m)[0][1] * (m)[1][2] * (m)[2][0] \
228  + (m)[0][2] * (m)[1][0] * (m)[2][1] \
229  - (m)[0][0] * (m)[1][2] * (m)[2][1] \
230  - (m)[0][1] * (m)[1][0] * (m)[2][2] \
231  - (m)[0][2] * (m)[1][1] * (m)[2][0] )
232 
236 #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)
237 #define init_matrixInternal(m,a,b,c,d,e,f,g,h,i) \
238  {(m)[0][0] = a; (m)[0][1] = b; (m)[0][2] = c; \
239  (m)[1][0] = d; (m)[1][1] = e; (m)[1][2] = f; \
240  (m)[2][0] = g; (m)[2][1] = h; (m)[2][2] = i; }
241 
245 #define multSM(r,s,mm) multSMInternal((r).m,s,(mm).m)
246 #define multSMInternal(r,s,m) \
247  {(r)[0][0] = (s) * (m)[0][0]; (r)[0][1] = (s) * (m)[0][1]; (r)[0][2] = (s) * (m)[0][2]; \
248  (r)[1][0] = (s) * (m)[1][0]; (r)[1][1] = (s) * (m)[1][1]; (r)[1][2] = (s) * (m)[1][2]; \
249  (r)[2][0] = (s) * (m)[2][0]; (r)[2][1] = (s) * (m)[2][1]; (r)[2][2] = (s) * (m)[2][2]; }
250 
254 #define multSI(r,s) multSIInternal((r).m,s)
255 #define multSIInternal(r,s) \
256  {(r)[0][0] = (s); (r)[0][1] = 0; (r)[0][2] = 0; \
257  (r)[1][0] = 0; (r)[1][1] = (s); (r)[1][2] = 0; \
258  (r)[2][0] = 0; (r)[2][1] = 0; (r)[2][2] = (s); }
259 
260 
264 #define addM(r,a,b) addMInternal((r).m,(a).m,(b).m)
265 #define addMInternal(r,a,b) \
266  {(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]; \
267  (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]; \
268  (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]; }
269 
274 #define multVVt(r,a,b) multVVtInternal((r).m,(a).v,(b).v)
275 #define multVVtInternal(r,a,b) \
276  {(r)[0][0] = (a)[0] * (b)[0]; (r)[0][1] = (a)[0] * (b)[1]; (r)[0][2] = (a)[0] * (b)[2]; \
277  (r)[1][0] = (a)[1] * (b)[0]; (r)[1][1] = (a)[1] * (b)[1]; (r)[1][2] = (a)[1] * (b)[2]; \
278  (r)[2][0] = (a)[2] * (b)[0]; (r)[2][1] = (a)[2] * (b)[1]; (r)[2][2] = (a)[2] * (b)[2]; }
279 
280 
284 #define multMV(r,mm,vv) multMVInternal((r).v,(mm).m,(vv).v)
285 #define multMVInternal(r,m,v) \
286  {(r)[0] = (m)[0][0] * (v)[0] + (m)[0][1] * (v)[1] + (m)[0][2] * (v)[2]; \
287  (r)[1] = (m)[1][0] * (v)[0] + (m)[1][1] * (v)[1] + (m)[1][2] * (v)[2]; \
288  (r)[2] = (m)[2][0] * (v)[0] + (m)[2][1] * (v)[1] + (m)[2][2] * (v)[2]; }
289 
290 
294 #define inv(r,mm) invInternal((r).m,(mm).m)
295 #define invInternal(r,m) \
296  {(r)[0][0] = (m)[1][1] * (m)[2][2] - (m)[1][2] * (m)[2][1]; \
297  (r)[0][1] = - (m)[0][1] * (m)[2][2] + (m)[0][2] * (m)[2][1]; \
298  (r)[0][2] = (m)[0][1] * (m)[1][2] - (m)[0][2] * (m)[1][1]; \
299  (r)[1][0] = - (m)[1][0] * (m)[2][2] + (m)[1][2] * (m)[2][0]; \
300  (r)[1][1] = (m)[0][0] * (m)[2][2] - (m)[0][2] * (m)[2][0]; \
301  (r)[1][2] = - (m)[0][0] * (m)[1][2] + (m)[0][2] * (m)[1][0]; \
302  (r)[2][0] = (m)[1][0] * (m)[2][1] - (m)[1][1] * (m)[2][0]; \
303  (r)[2][1] = - (m)[0][0] * (m)[2][1] + (m)[0][1] * (m)[2][0]; \
304  (r)[2][2] = (m)[0][0] * (m)[1][1] - (m)[0][1] * (m)[1][0]; \
305  multSMInternal(r,1/detInternal(m),r); }
306 
310 #define multMM(r,a,b) multMMInternal((r).m,(a).m,(b).m)
311 #define multMMInternal(r,a,b) \
312  {(r)[0][0] = (a)[0][0] * (b)[0][0] + (a)[0][1] * (b)[1][0] + (a)[0][2] * (b)[2][0]; \
313  (r)[0][1] = (a)[0][0] * (b)[0][1] + (a)[0][1] * (b)[1][1] + (a)[0][2] * (b)[2][1]; \
314  (r)[0][2] = (a)[0][0] * (b)[0][2] + (a)[0][1] * (b)[1][2] + (a)[0][2] * (b)[2][2]; \
315  (r)[1][0] = (a)[1][0] * (b)[0][0] + (a)[1][1] * (b)[1][0] + (a)[1][2] * (b)[2][0]; \
316  (r)[1][1] = (a)[1][0] * (b)[0][1] + (a)[1][1] * (b)[1][1] + (a)[1][2] * (b)[2][1]; \
317  (r)[1][2] = (a)[1][0] * (b)[0][2] + (a)[1][1] * (b)[1][2] + (a)[1][2] * (b)[2][2]; \
318  (r)[2][0] = (a)[2][0] * (b)[0][0] + (a)[2][1] * (b)[1][0] + (a)[2][2] * (b)[2][0]; \
319  (r)[2][1] = (a)[2][0] * (b)[0][1] + (a)[2][1] * (b)[1][1] + (a)[2][2] * (b)[2][1]; \
320  (r)[2][2] = (a)[2][0] * (b)[0][2] + (a)[2][1] * (b)[1][2] + (a)[2][2] * (b)[2][2]; }
321 
322 
331 
332 // Number of iteration of the jacobi iterative method.
333 #define ITER_SJ 10
334 
335  Vector x;
336 
337  // The new approximation.
338  Vector x_1;
339 
340  // Clear the vectors.
341  clearV(x);
342  clearV(x_1);
343 
344  int iter;
345  int i,j;
346 
347  // Iterative method.
348  for(iter=0; iter<ITER_SJ; iter++){
349  for(i=0;i<D;i++){
350 
351  if( A.m[i][i] == 0 ) continue;
352 
353  // calculate the new value.
354  x_1.v[i] = 0.0;
355  for(j=0;j<D;j++){
356  if(i != j) x_1.v[i] += A.m[i][j] * x.v[j];
357  }
358  x_1.v[i] = (b.v[i] - x_1.v[i]) / A.m[i][i];
359 
360  }
361  cpyV(x,x_1);
362  }
363 
364  return x;
365 }
366 
367 
368 
380 
381  Vector x;
382 
383  double detA = det(A);
384 
385  // If the matrix has a determinant, we can get the inverse.
386  if(detA != 0){
387 
388  Matrix InvA;
389  inv(InvA,A);
390  multMV(x,InvA,b);
391  } else {
392  //printf("#Warning: singular matrix\n");
393  x = solve_jacobi(A,b);
394  }
395  return x;
396 }
397 
402 
407 }
408 
413 
416  hit_tileFree(e);
420 
421  hit_shapeFree(shape);
422 }
423 
424 
429 void print_help(char * name){
430  printf("%s [-n NEWTON METHOD ITERATIONS] [-j JACOBI METHOD ITERATIONS] FILE \n",name);
431 
432  printf(" -n NEWTON METHOD ITERATIONS number of iterations (default 10)\n");
433  printf(" -j JACOBI METHOD ITERATIONS number of iterations (default 100)\n");
434  printf(" FILE input graph file\n");
435 
436 }
437 
438 
444 
445  srandom(0);
446  if(hit_Rank != 0) return;
447 
448  int n = hit_cShapeNvertices(shape);
449 
450  int j,i;
451 
452  // Set the coordinates
453  for(j=0;j<D;j++){
454  for(i=0;i<n;i++){
455  double number = (double) random();
456  hit_gcTileVertexAt(graph,i).v[j] = number;
457  }
458  }
459 
460  srandom(0);
461  for(i=0;i<n;i++){
462 
463  if ((random() % 100) < P_FIXED){
464  hit_gcTileVertexAt(fixed,i) = 1;
465  } else {
466  hit_gcTileVertexAt(fixed,i) = 0;
467  }
468  }
469 
470 
471 }
472 
473 
474 
475 
476 
477 
481 void init_graph(int argc, char ** argv){
482 
483  // Define the spring structure.
484  shape = HIT_CSR_SHAPE_NULL;
485 
486  // Argument.
487  char c;
488 
489  // Parse the command-line arguments.
490  while ((c = (char) getopt (argc, argv, "hn:j:")) != -1)
491  switch (c) {
492  case 'n':
493  sscanf(optarg,"%d",&iter1);
494  break;
495  case 'j':
496  sscanf(optarg,"%d",&iter2);
497  break;
498  case 'h':
499  print_help(argv[0]);
500  exit(EXIT_SUCCESS);
501  default:
502  abort ();
503  break;
504  }
505 
506 
507  if(optind == (argc-1)){
508 
509  // Graph file name.
510  char * graph_file = argv[optind];
511 
512  printf("# Graph: %s\n",graph_file);
513  printf("# iterations %d x %d\n",iter1,iter2);
514 
515  shape = hit_fileHBRead(graph_file);
516 
517  printf("# Vertices %d\n",hit_cShapeNvertices(shape));
518 
519  // Allocate memory
522 
523 
524  // Random coordinates
526 
527  int i;
528  int nfixed = 0;
529  for(i=0;i<hit_cShapeNvertices(shape);i++){
530 
531  nfixed += hit_gcTileVertexAt(fixed,i);
532  }
533 
534  printf("# N fixed: %d\n",nfixed);
535 
536 
537 
538  } else {
539 
540  iter1 = iter2 = 3;
541 
542  hit_cShapeAddEdge2(&shape,0,1);
543  hit_cShapeAddEdge2(&shape,1,2);
544  hit_cShapeAddEdge2(&shape,2,3);
545  hit_cShapeAddEdge2(&shape,4,5);
546  hit_cShapeAddEdge2(&shape,5,6);
547  hit_cShapeAddEdge2(&shape,6,7);
548  hit_cShapeAddEdge2(&shape,0,4);
549  hit_cShapeAddEdge2(&shape,1,5);
550  hit_cShapeAddEdge2(&shape,2,6);
551  hit_cShapeAddEdge2(&shape,3,7);
552  hit_cShapeAddEdge2(&shape,1,4);
553  hit_cShapeAddEdge2(&shape,2,5);
554 
555  // Include each (i,i) edge.
556  hit_cShapeAddEdge2(&shape,0,0);
557  hit_cShapeAddEdge2(&shape,1,1);
558  hit_cShapeAddEdge2(&shape,2,2);
559  hit_cShapeAddEdge2(&shape,3,3);
560  hit_cShapeAddEdge2(&shape,4,4);
561  hit_cShapeAddEdge2(&shape,5,5);
562  hit_cShapeAddEdge2(&shape,6,6);
563  hit_cShapeAddEdge2(&shape,7,7);
564 
565  // Allocate memory
568 
569  // Set initial positions
570  Vector v0 = {{0,0,0}};
571  hit_gcTileVertexAt(graph,0) = v0;
572  Vector v1 = {{100,0,0}};
573  hit_gcTileVertexAt(graph,1) = v1;
574  Vector v2 = {{200,0,0}};
575  hit_gcTileVertexAt(graph,2) = v2;
576  Vector v3 = {{300,0,0}};
577  hit_gcTileVertexAt(graph,3) = v3;
578  Vector v4 = {{0,100,200}};
579  hit_gcTileVertexAt(graph,4) = v4;
580  Vector v5 = {{100,100,0}};
581  hit_gcTileVertexAt(graph,5) = v5;
582  Vector v6 = {{200,100,0}};
583  hit_gcTileVertexAt(graph,6) = v6;
584  Vector v7 = {{300,100,0}};
585  hit_gcTileVertexAt(graph,7) = v7;
586 
587  // Set the fix/free status
588  hit_gcTileVertexAt(fixed,0) = 1;
589  hit_gcTileVertexAt(fixed,1) = 0;
590  hit_gcTileVertexAt(fixed,2) = 0;
591  hit_gcTileVertexAt(fixed,3) = 1;
592  hit_gcTileVertexAt(fixed,4) = 1;
593  hit_gcTileVertexAt(fixed,5) = 0;
594  hit_gcTileVertexAt(fixed,6) = 0;
595  hit_gcTileVertexAt(fixed,7) = 1;
596 
597  }
598 
599 }
600 
605 double gradient_norm(){
606 
607  double norm = 0.0;
608 
609  int vertex;
610 
611  hit_cShapeVertexIterator(vertex,shape) {
612  if (!hit_gcTileVertexAt(fixed,vertex)){
613 
614  calculate_GH(vertex);
615 
616  Vector Gi = hit_gcTileVertexAt(gradient,vertex);
617 
618  int d;
619  for(d=0;d<D;d++){
620  norm += pow(Gi.v[d],2);
621  }
622  }
623  }
624 
625  return sqrt(norm);
626 }
627 
628 
629 
630 
634 int main(int argc, char ** argv) {
635 
636  int i;
637  int vertex;
638  HitClock init_time, comp_time;
639 
640  hit_comInit(&argc,&argv);
641 
642  hit_clockStart(init_time);
643  init_graph(argc,argv);
644 
645  // Allocate memory of all the variables
646  init_structures();
647 
648  // printf("# Initial gradient norm: %lf\n",gradient_norm());
649  hit_clockStop(init_time);
650  hit_clockStart(comp_time);
651 
652  // Main loop for the newton method
653  for(i=0;i<iter1;i++){
654  printf("# iter: %d/%d\n",i,iter1);
655 
656  // Iterate trough all the vertices and
657  // obtain the gradient and the hessian.
658  hit_cShapeVertexIterator(vertex,shape) {
659  if (!hit_gcTileVertexAt(fixed,vertex))
660  calculate_GH(vertex);
661  }
662 
663 
664  // The initial position displacement is zero.
666 
667  // Loop for the jacobi method to solve the system.
668  int j;
669  for(j=0;j<iter2;j++){
670 
671  // Perform a iteration.
672  hit_cShapeVertexIterator(vertex,shape) {
673  if (!hit_gcTileVertexAt(fixed,vertex))
674  solve_system_iter(vertex);
675  }
676 
677  // Update the displacement.
679 
680  }
681 
682  // Update the position with the final displacement.
683  hit_cShapeVertexIterator(vertex,shape) {
684  if (!hit_gcTileVertexAt(fixed,vertex)){
685  Vector * current = &(hit_gcTileVertexAt(graph,vertex));
686  subV(*current,*current,hit_gcTileVertexAt(e,vertex));
687  }
688  }
689 
690  }
691 
692  hit_clockStop(comp_time);
693 
694  printf("# Final gradient norm: %lf\n",gradient_norm());
695  printf("# Init time: %lf\n",hit_clockGetSeconds(init_time));
696  printf("# Comp time: %lf\n",hit_clockGetSeconds(comp_time));
697  printf("# Total time: %lf\n",hit_clockGetSeconds(init_time)+hit_clockGetSeconds(comp_time));
698 
699  free_structures();
700  hit_comFinalize();
701 
702  return 0;
703 }
704 
705 
706 
710 Vector g(int i, int j){
711 
712  Vector gij;
713 
714  if(i == j){
715  clearV(gij);
716  return gij;
717  }
718 
719  // Get the two vectors
720  Vector r_i = (hit_gcTileVertexAt(graph,i));
721  Vector r_j = (hit_gcTileVertexAt(graph,j));
722 
723  // Auxiliar variables.
724  Vector v;
725 
726  // v = r_i - r_j
727  subV(v,r_i,r_j);
728  // n = ||v||
729  double n = norm(v);
730 
731  if( n == 0){
732  clearV(gij);
733  return gij;
734  }
735 
736  // gij = - K * ((L - n) / n) * (v)
737  double scalar = - K * ((L - n) / n);
738  multSV(gij,scalar,v);
739 
740  return gij;
741 }
742 
743 
747 Matrix W(int i, int j){
748 
749  Matrix Wij;
750 
751  if(i == j){
752  clearM(Wij);
753  return Wij;
754  }
755 
756  // Get the two vectors
757  Vector r_i = (hit_gcTileVertexAt(graph,i));
758  Vector r_j = (hit_gcTileVertexAt(graph,j));
759 
760  // Auxiliary variables.
761  Vector v;
762  Matrix a, b;
763  double scalar_a, scalar_b;
764 
765  // v = r_i - r_j;
766  subV(v,r_i,r_j);
767  // n = ||v||
768  double n = norm(v);
769 
770  if( n == 0){
771  clearM(Wij);
772  return Wij;
773  }
774 
775  // h = -((L-n)/n) * eye(D) + (L/(n^3)) * (v * v')
776  scalar_a = -((L-n)/(n));
777  multSI(a,scalar_a);
778  scalar_b = (L/(pow(n,3)));
779  multVVt(b,v,v);
780  multSM(b,scalar_b,b);
781  addM(Wij,a,b);
782  // h = K * h
783  multSM(Wij,K,Wij);
784 
785  return Wij;
786 }
787 
788 
794 void calculate_GH(int vertex){
795 
796  Vector gs;
797  Matrix Ws;
798 
799  clearV(gs);
800  clearM(Ws);
801 
802  int edge;
803  hit_cShapeEdgeIterator(edge,shape,vertex){
804 
805  // Get the neighbor.
806  int nghb_index = hit_cShapeEdgeTarget(shape,edge);
807 
808  // Calculate gradient[i]
809  Vector gij = g(vertex,nghb_index);
810  addV(gs,gs,gij);
811 
812  // Calculate part of the hessian[i][i]
813  Matrix Wij = W(vertex,nghb_index);
814  addM(Ws,Ws,Wij);
815 
816  // Calculate hessian[i][j] if j is variable.
817  //Node nj = hit_gcTileVertexAt(graph,nghb_index);
818  if(!hit_gcTileVertexAt(fixed,nghb_index)){
819 
820  Matrix Hij = W(nghb_index,vertex);
821  multSM(Hij,-1,Hij);
822 
823  hit_gcTileEdgeIteratorAt(hessian,vertex,edge) = Hij;
824  }
825  }
826 
827  hit_gcTileVertexAt(gradient,vertex) = gs;
828  hit_gcTileEdgeAt(hessian,vertex,vertex) = Ws;
829 
830 }
831 
837 void solve_system_iter(int vertex){
838 
839  Matrix Hii = hit_gcTileEdgeAt(hessian,vertex,vertex);
840  Vector Gi = hit_gcTileVertexAt(gradient,vertex);
841 
842  Vector new_ei;
843 
844  // Aux sum vector.
845  Vector sum;
846  clearV(sum);
847 
848  int edge;
849  hit_cShapeEdgeIterator(edge,shape,vertex){
850 
851  // Get the neighbor.
852  int nghb_index = hit_cShapeEdgeTarget(shape,edge);
853  if(nghb_index == vertex) continue;
854 
855  // If the node is fixed, it is no part of the equation system.
856  //Node nj = hit_gcTileVertexAt(graph,nghb_index);
857  if(hit_gcTileVertexAt(fixed,nghb_index)) continue;
858 
859  // Get the j displacement and ij hessian part.
860  Vector ej = hit_gcTileVertexAt(e,nghb_index);
861  Matrix Hij = hit_gcTileEdgeIteratorAt(hessian,vertex,edge);
862 
863  // Aux par vector.
864  Vector part;
865 
866  // Calculate and add the contribution of vertex j.
867  // SUM_{i not j} ( Hij e_j )
868  multMV(part,Hij,ej);
869  addV(sum,sum,part);
870 
871  }
872 
873  // Gi - SUM_{i not j} ( Hij e_j )
874  subV(sum,Gi,sum);
875 
876  // e_i = H_ii^(-1) * (Gi - SUM_{i not j} ( Hij e_j ) )
877  new_ei = solve(Hii,sum);
878  hit_gcTileVertexAt(new_e,vertex) = new_ei;
879 }
880 
881 
#define subV(r, a, b)
Vector g(int i, int j)
void hit_gcTileCopyVertices(void *destP, void *srcP)
Definition: hit_gctile.c:145
void print_help(char *name)
Definition: heat.c:234
#define hit_tileNewType(baseType)
Definition: hit_tile.h:163
#define hit_cShapeEdgeTarget(s, edge)
#define A
Definition: mg.c:64
double v[D]
Definition: spring_bitmap.c:86
#define multSV(r, s, vv)
#define hit_Rank
Definition: hit_com.h:140
#define D
int iter2
Definition: spring_bitmap.c:80
#define hit_cShapeVertexIterator(var, shape)
Definition: hit_cshape.h:392
int iter
Definition: mmult_bit.c:68
Matrix W(int i, int j)
#define multSI(r, s)
HitTile_Matrix hessian
#define ITER2
#define addV(r, a, b)
#define L
void solve_system_iter(int vertex)
#define hit_cShapeEdgeIterator(var, shape, vertex)
Definition: hit_cshape.h:416
#define hit_gcTileEdgeIteratorAt(var, vertex, edge_index)
HitTile_Vector graph
#define HIT_EDGES
Definition: hit_tile.h:214
void hit_gcTileClearVertices(void *varP)
Definition: hit_gctile.c:132
Vector solve(Matrix A, Vector b)
void free_structures()
#define hit_gcTileDomainShapeAlloc(var, baseType, shape, allocOpts)
Definition: hit_gctile.h:99
#define det(mm)
HitTile_int fixed
#define x
#define hit_gcTileVertexAt(var, vertex)
Definition: hit_gctile.h:157
#define hit_cShapeNvertices(shape)
#define ITER_SJ
#define norm(vv)
#define addM(r, a, b)
Vector solve_jacobi(Matrix A, Vector b)
#define cpyV(r, a)
void hit_comInit(int *pargc, char **pargv[])
Definition: hit_com.c:111
HitShape shape
#define clearM(mm)
#define multVVt(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
#define hit_clockStop(c)
Definition: hit_utils.h:109
int main(int argc, char *argv[])
Definition: cannonAsync.c:62
#define m(a, b, c)
HitTile_Vector new_e
#define inv(r, mm)
#define clearV(vv)
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)
double gradient_norm()
#define P_FIXED
#define hit_tileFree(var)
Definition: hit_tile.h:369
void hit_shapeFree(HitShape s)
Definition: hit_shape.c:84
#define hit_gcTileEdgeAt(var, pos1, pos2)
Definition: hit_gctile.h:167
#define multMV(r, mm, vv)
#define K
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 HIT_VERTICES
Definition: hit_tile.h:208
#define ITER1
#define hit_clockStart(c)
Definition: hit_utils.h:87
double m[D][D]
Definition: spring_bitmap.c:93
#define multSM(r, s, mm)