LoAdSG
cg_method.h
1//
2// Created by scherner on 26.07.21.
3//
4
5#ifndef GRUN_CG_METHOD_H
6#define GRUN_CG_METHOD_H
7
8
9#include <sys/time.h>
10#include <unordered_set>
11#include "../MatrixVectorMultiplication/MatrixVectorHomogen.h"
12#include "../MatrixVectorMultiplication/MatrixVectorInhomogen.h"
13
14class CG {
15public:
16 static bool solveNeumann(int maxIteration,
17 double eps,
18 VectorSparseG &x,
19
20 VectorSparseG &f,
21
22
23
24 int *iterations,
25 MatrixVectorInhomogen& matrix);
26
27 template<class Stencil>
28 static bool solveHomogen(double eps,
29 VectorSparseG &x,
30 VectorSparseG &f,
31 int *iterations,
32 VectorSparseG &usol,
33 MatrixVectorHomogen& matrix,Stencil stencil,double* time_precon);
34
35 static const int maxIteration = 500;
36 static double error[MaximumDepth][maxIteration];
37 static int count;
38
39
40 static void printError() {
41 ofstream Datei;
42 Datei.open("../results/error_CG.gnu", std::ios::out);
43
44 for (int i = 0; i < MaximumDepth; i++) {
45 if (error[i][0] != 0) {
46 Datei << "\t" << error[i][0];
47 cout << "\t" << error[i][0];
48 }
49 }
50 Datei << endl;
51 cout << endl;
52 for (int i = 1; i < maxIteration; i++) {
53 Datei << i;
54 cout << i;
55 for (int j = 0; j < MaximumDepth; j++) {
56 if (error[j][0] != 0) {
57 Datei << "\t" << error[j][i];
58 cout << "\t " << error[j][i];
59 }
60
61
62 }
63
64 Datei << "\n";
65 cout << "\n";
66
67 }
68 Datei << endl;
69 Datei.close();
70 }
71};
72
73
74class OperatorT {
75public:
76 static void apply(VectorSparseG &input, VectorSparseG &output);
77
78 static void applyTranspose(VectorSparseG &input, VectorSparseG &output);
79};
80
81template<class StencilClass>
82class Preconditioning{
83public:
84 Preconditioning(VectorSparseG& z_, StencilClass stencilClass) : z(z_), boundary(false) {
85
86 if(stencilClass.getTypeMatrixVectorMultiplication()==StencilOnTheFly) precond3(stencilClass);
87 else{
88 precondLocalStiffnessMatrix(stencilClass);
89 }
90
91 MPI_Barrier(MPI_COMM_WORLD);
92 }
93
94
95 void precondMV(StencilClass stencilClass){
96 int world_rank;
97 int world_size;
98#ifdef MY_MPI_ON
99
100 MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
101
102 MPI_Comm_size(MPI_COMM_WORLD, &world_size);
103#endif
104
105
106
107 VectorSparseG v(*z.getSparseGrid());
108 VectorSparseG Av(*z.getSparseGrid());
109 MultiLevelAdaptiveSparseGrid mgrid(z.getSparseGrid());
110 MatrixVectorHomogen matrixVectorHomogen(*z.getSparseGrid(), mgrid, 1, 0);
111
112 for (size_t i = 0; i < z.getSparseGrid()->getMaximalOccupiedSecondTable(); i++) {
113 if (z.getSparseGrid()->getActiveTable()[i]) {
114
115 IndexDimension Index = z.getSparseGrid()->getIndexOfTable(i);
116
117
118 v.setValue(Index,1);
119 matrixVectorHomogen.multiplication(v,Av,stencilClass);
120
121
122
123
124 v = 0.0;
125 double coeff=Av.getValue(i);
126 if(world_rank==0) {
127 Index.PrintCoord();
128 cout << " " << coeff << endl;
129 }
130 z.setValue(i, coeff);
131
132 }
133
134
135 }
136
137
138
139 }
140
141 void precond1(StencilClass stencilClass){
142#ifdef MY_MPI_ON
143 int world_rank;
144 MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
145 int world_size;
146 MPI_Comm_size(MPI_COMM_WORLD, &world_size);
147#endif
148
149
150
151 VectorSparseG v(*z.getSparseGrid());
152
153 for (size_t i = 0; i < z.getSparseGrid()->getMaximalOccupiedSecondTable(); i++) {
154 if (z.getSparseGrid()->getActiveTable()[i]) {
155#ifdef MY_MPI_ON
156 if (world_rank == i % world_size ) {
157
158#endif
159 IndexDimension Index = z.getSparseGrid()->getIndexOfTable(i);
160
161 MultiDimFiveCompass mc_outer;
162 bool todo[mc_outer.getMaxShift()];
163
164 IndexDimension index_mc[mc_outer.getMaxShift()];
165 double basis_coeff_mc[mc_outer.getMaxShift()];
166
167 for (int i = 0; i < mc_outer.getMaxShift(); i++) {
168 todo[i] = false;
169 }
170 Depth T(Index);
171 stencilClass.initialize(T);
172
173
174 for (MultiDimFiveCompass mc; mc.goon(); ++mc) {
175 double basis_coeff = 1.0;
176 bool next = true;
177 IndexDimension J;
178 J = Index.nextFive(&mc, T, &basis_coeff, &next);
179
180
181 if (next) {
182 todo[mc.getShiftNumber()] = true;
183 index_mc[mc.getShiftNumber()] = J;
184 basis_coeff_mc[mc.getShiftNumber()] = basis_coeff;
185
186 v.setValue(J, basis_coeff);
187
188
189 }
190
191
192 }
193
194 double coeff = 0.0;
195
196 #pragma omp parallel for schedule(dynamic)
197 for (int j = 0; j < mc_outer.getMaxShift(); j++) {
198
199
200 if (todo[j]) {
201
202 IndexDimension J = index_mc[j];
203 double stencil = applyStencilLocal(J, v, T, stencilClass);
204
205
206#pragma omp critical
207 stencil *= basis_coeff_mc[j];
208
209#pragma omp critical
210 coeff += stencil;
211
212
213 }
214
215 }
216
217
218 v = 0.0;
219
220 z.setValue(i, coeff);
221#ifdef MY_MPI_ON
222 }
223#endif
224 }
225
226
227 }
228
229#ifdef MY_MPI_ON
230 MPI_Allreduce(MPI_IN_PLACE,z.getDatatableVector(),z.getSparseGrid()->getMaximalOccupiedSecondTable(),MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
231#endif
232
233 }
234
235
236 void precond2(StencilClass stencilClass){
237
238#pragma omp parallel
239 {
240 VectorSparseG v(*z.getSparseGrid());
241 #pragma omp for schedule(runtime)
242 for (size_t i = 0; i < z.getSparseGrid()->getMaximalOccupiedSecondTable(); i++) {
243 if (z.getSparseGrid()->getActiveTable()[i]) {
244
245 IndexDimension Index = z.getSparseGrid()->getIndexOfTable(i);
246
247 MultiDimFiveCompass mc_outer;
248 bool todo[mc_outer.getMaxShift()];
249
250
251 IndexDimension index_mc[mc_outer.getMaxShift()];
252 double basis_coeff_mc[mc_outer.getMaxShift()];
253
254 for (int j = 0; j < mc_outer.getMaxShift(); j++) {
255 todo[j] = false;
256 }
257 Depth T(Index);
258
259
260 for (MultiDimFiveCompass mc; mc.goon(); ++mc) {
261 double basis_coeff = 1.0;
262 bool next = true;
263 IndexDimension J;
264 J = Index.nextFive(&mc, T, &basis_coeff, &next);
265
266
267 if (next) {
268
269 todo[mc.getShiftNumber()] = true;
270 index_mc[mc.getShiftNumber()] = J;
271 basis_coeff_mc[mc.getShiftNumber()] = basis_coeff;
272
273 v.setValue(J, basis_coeff);
274
275
276 }
277
278
279 }
280
281 double coeff = 0.0;
282
283
284 for (int j = 0; j < mc_outer.getMaxShift(); j++) {
285
286
287 if (todo[j]) {
288
289 IndexDimension J = index_mc[j];
290 double stencil = applyStencilLocal(J, v, T, stencilClass);
291
292
293 stencil *= basis_coeff_mc[j];
294
295
296 coeff += stencil;
297
298
299
300 }
301
302 }
303
304
305 v = 0.0;
306
307 z.setValue(i, coeff);
308
309 }
310 }
311
312
313 }
314
315
316
317
318
319 }
320
321
322 void precond3(StencilClass stencilClass){
323#ifdef MY_MPI_ON
324 int world_rank;
325 MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
326 int world_size;
327 MPI_Comm_size(MPI_COMM_WORLD, &world_size);
328#endif
329
330
331
332 int iter_i=0;
333 DepthList depthList(*z.getSparseGrid());
334
335
336
337
338 for (auto it = depthList.begin_all(); it != depthList.end_all(); ++it) {
339
340
341 Depth T = *it;
342
343 SingleDepthHashGrid &depthGrid = z.getSparseGrid()->getMultiDepthHashGrid()->getGridForDepth(T);
344 const auto &mapping = depthGrid._mapPosToGridPos;
345 const auto mapping_size = mapping.size();
346
347
348 if (depthGrid.getNumberOfEntries() > 0) {
349
350
351#ifdef MY_MPI_ON
352 if (world_rank == iter_i % world_size ) {
353#endif
354
355 stencilClass.initialize(T);
356#pragma omp parallel
357 {
358 VectorSparseG v(*z.getSparseGrid());
359 #pragma omp for schedule(dynamic)
360 for (size_t i = 0; i < mapping_size; i++) {
361
362
363
364 if (z.getSparseGrid()->getActiveTable()[mapping[i]]) {
365 IndexDimension Index = depthGrid._map.getIndexOfTable(i);
366
367 MultiDimFiveCompass mc_outer;
368 bool todo[mc_outer.getMaxShift()];
369
370 IndexDimension index_mc[mc_outer.getMaxShift()];
371 double basis_coeff_mc[mc_outer.getMaxShift()];
372
373 for (int j = 0; j < mc_outer.getMaxShift(); j++) {
374 todo[j] = false;
375 }
376
377
378 for (MultiDimFiveCompass mc; mc.goon(); ++mc) {
379 double basis_coeff = 1.0;
380 bool next = true;
381 IndexDimension J;
382 if (boundary) {
383 J = Index.nextFive_Neumann(&mc, T, &basis_coeff, &next);
384 } else
385 J = Index.nextFive(&mc, T, &basis_coeff, &next);
386
387
388 if (next) {
389 todo[mc.getShiftNumber()] = true;
390 index_mc[mc.getShiftNumber()] = J;
391 basis_coeff_mc[mc.getShiftNumber()] = basis_coeff;
392
393 v.setValue(J, basis_coeff);
394
395
396 }
397
398
399 }
400
401 double coeff = 0.0;
402
403
404 for (int j = 0; j < mc_outer.getMaxShift(); j++) {
405
406
407 if (todo[j]) {
408
409 IndexDimension J = index_mc[j];
410 double stencil = applyStencilLocal(J, v, T, stencilClass);
411 stencil *= basis_coeff_mc[j];
412
413
414 coeff += stencil;
415
416
417 }
418
419 }
420
421
422 v = 0.0;
423
424 z.setValue(mapping[i], coeff);
425
426 }
427
428
429
430 }
431
432
433 }
434
435
436#ifdef MY_MPI_ON
437 }
438#endif
439 }
440
441 iter_i++;
442 }
443
444#ifdef MY_MPI_ON
445
446 MPI_Allreduce(MPI_IN_PLACE,z.getDatatableVector(),z.getSparseGrid()->getMaximalOccupiedSecondTable(),MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
447#endif
448
449
450
451 }
452
453
454 void precond4(StencilClass stencilClass){
455#pragma omp parallel
456 {
457 VectorSparseG v(z.getSparseGrid());
458 VectorSparseG u(z.getSparseGrid());
459#pragma omp for schedule(dynamic)
460 for (unsigned long k = 0; k < z.getSparseGrid()->getMaximalOccupiedSecondTable(); k++) {
461
462 if (z.getSparseGrid()->getActiveTable()[k]) {
463
464
465 IndexDimension Index = z.getSparseGrid()->getIndexOfTable(k);
466 Depth Tneu(Index);
467 stencilClass.initialize(Tneu);
468
469
470 for (MultiDimFiveCompass mc; mc.goon(); ++mc) {
471 double basis_coeff = 1.0;
472 bool next = true;
473 IndexDimension J;
474 if (boundary) {
475 J = Index.nextFive_Neumann(&mc, Tneu, &basis_coeff, &next);
476 } else
477 J = Index.nextFive(&mc, Tneu, &basis_coeff, &next);
478
479
480 if (next) {
481 v.setValue(J, basis_coeff);
482 }
483
484
485
486
487 }
488
489
490 for (MultiDimFiveCompass mc; mc.goon(); ++mc) {
491 double basis_coeff = 1.0;
492 bool next = true;
493 IndexDimension J;
494 if (boundary) {
495 J = Index.nextFive_Neumann(&mc, Tneu, &basis_coeff, &next);
496 } else
497 J = Index.nextFive(&mc, Tneu, &basis_coeff, &next);
498
499
500 if (next) {
501 basis_coeff = 0.0;
502
503 // basis_coeff = CalcStencilValue_Boundary(J, Tneu, v, Stiffness, stencil);
504 basis_coeff = applyStencilLocal(J,v,Tneu,stencilClass);
505
506 u.setValue(J, basis_coeff);
507 }
508
509 }
510
511
512 double coeff = 0.0;
513 for (MultiDimFiveCompass mc; mc.goon(); ++mc) {
514 double basis_coeff = 1.0;
515 bool next = true;
516 IndexDimension J;
517 if (boundary) {
518 J = Index.nextFive_Neumann(&mc, Tneu, &basis_coeff, &next);
519 } else
520 J = Index.nextFive(&mc, Tneu, &basis_coeff, &next);
521
522 if (next) {
523 coeff = coeff + u.getValue(J) * basis_coeff;
524
525 }
526 v.setValue(J, 0.0);
527 u.setValue(J, 0.0);
528
529 }
530 v=0.0;
531 u=0.0;
532
533 z.setValue(k, coeff);
534
535 }
536 }
537
538 }
539
540 }
541
542
543
544 void precondLocalStiffnessMatrix(StencilClass stencilClass){
545 int world_rank;
546#ifdef MY_MPI_ON
547
548 MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
549 int world_size;
550 MPI_Comm_size(MPI_COMM_WORLD, &world_size);
551#endif
552
553
554
555 int iter_i=0;
556 DepthList depthList(*z.getSparseGrid());
557
558
559
560
561 for (auto it = depthList.begin_all(); it != depthList.end_all(); ++it) {
562
563
564 Depth T = *it;
565 int node = stencilClass.getNodeForDepth(T);
566 if(world_rank==node){
567
568 SingleDepthHashGrid &depthGrid = z.getSparseGrid()->getMultiDepthHashGrid()->getGridForDepth(T);
569 const auto &mapping = depthGrid._mapPosToGridPos;
570 const auto mapping_size = mapping.size();
571
572
573 if (depthGrid.getNumberOfEntries() > 0) {
574
575
576#pragma omp parallel
577 {
578 VectorSparseG v(*z.getSparseGrid());
579 VectorSparseG v_output(*z.getSparseGrid());
580
581#pragma omp for schedule(dynamic)
582 for (size_t i = 0; i < mapping_size; i++) {
583
584
585 if (z.getSparseGrid()->getActiveTable()[mapping[i]]) {
586 IndexDimension Index = depthGrid._map.getIndexOfTable(i);
587
588
589 MultiDimFiveCompass mc_outer;
590 bool todo[mc_outer.getMaxShift()];
591
592 IndexDimension index_mc[mc_outer.getMaxShift()];
593 double basis_coeff_mc[mc_outer.getMaxShift()];
594
595 for (int j = 0; j < mc_outer.getMaxShift(); j++) {
596 todo[j] = false;
597 }
598
599 std::vector<IndexDimension> IndicesVector;
600 for (MultiDimFiveCompass mc; mc.goon(); ++mc) {
601 double basis_coeff = 1.0;
602 bool next = true;
603 IndexDimension J;
604 if (boundary) {
605 J = Index.nextFive_Neumann(&mc, T, &basis_coeff, &next);
606 } else
607 J = Index.nextFive(&mc, T, &basis_coeff, &next);
608
609
610 if (next) {
611 todo[mc.getShiftNumber()] = true;
612 index_mc[mc.getShiftNumber()] = J;
613 basis_coeff_mc[mc.getShiftNumber()] = basis_coeff;
614
615 v.setValue(J, basis_coeff);
616 IndicesVector.push_back(J);
617
618
619 }
620
621
622 }
623
624 stencilClass.applyLocalStiffnessMatricesOnIndices_onNode(v, v_output, T, IndicesVector);
625
626 double coeff = 0.0;
627
628
629 for (int j = 0; j < mc_outer.getMaxShift(); j++) {
630
631
632 if (todo[j]) {
633
634 IndexDimension J = index_mc[j];
635 double stencil;
636
637
638 stencil = v_output.getValue(J);
639
640
641 stencil *= basis_coeff_mc[j];
642
643
644 coeff += stencil;
645
646
647 }
648
649 }
650
651
652 v = 0.0;
653 v_output = 0.0;
654
655
656 z.setValue(mapping[i], coeff);
657
658 }
659
660
661
662 }
663
664
665 }
666
667
668 }
669
670 }
671
672 iter_i++;
673 }
674
675#ifdef MY_MPI_ON
676
677 MPI_Allreduce(MPI_IN_PLACE,z.getDatatableVector(),z.getSparseGrid()->getMaximalOccupiedSecondTable(),MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
678#endif
679
680 }
681
682
683
684
685
686
687
688
689
690
691 double applyStencilLocal(IndexDimension &Index, const VectorSparseG &input,Depth &T, StencilClass stencilClass) {
692 double val = 0.0;
693 IndexDimension NextIndex;
694 //stencilClass.initialize(T);
695 for (MultiDimCompass mc; mc.goon(); ++mc) {
696 NextIndex = Index.nextThree2(&mc, T.returnTiefen());
697
698 if(mc.goon()) {
699 double value = 0.0;
700 value = stencilClass.returnValue(Index, mc);
701
702 val = val + value * input.getValue(NextIndex);
703
704
705 }
706
707 }
708 return val;
709 };
710
711
712
713 void apply(VectorSparseG* u){
714 *u = ((*u)*z);
715 }
716
717 void apply_inverse(VectorSparseG* u){
718 for (unsigned long k = 0; k < z.getSparseGrid()->getMaximalOccupiedSecondTable(); k++) {
719 if(z.getSparseGrid()->getActiveTable()[k]) {
720 double val1 = u->getValue(k);
721 double val2 = z.getValue(k);
722 u->setValue(k,double(val1/val2));
723
724
725 }
726 }
727 }
728
729 void apply_NEU(VectorSparseG* u){
730 for (unsigned long k = 0; k < z.getSparseGrid()->getMaximalOccupiedSecondTable(); k++) {
731 if(z.getSparseGrid()->getActiveTable()[k]) {
732 double val1 = u->getValue(k);
733 double val2 = z.getValue(k);
734 u->setValue(k, double(val1/ sqrt(val2)));
735 }
736 }
737 }
738
739 double getValue(unsigned long k){
740 return z.getValue(k);
741 }
742
743 VectorSparseG* getZ(){
744 return &z;
745 }
746
747
748private:
749 VectorSparseG& z;
750 bool boundary;
751
752
753
754};
755
756
757
758//INLINE TEMPLATE FUNCTIONS
759template<class Stencil>
760bool
761CG::solveHomogen(double eps, VectorSparseG &x,VectorSparseG &f,int *iterations,
762 VectorSparseG &usol, MatrixVectorHomogen& matrix,Stencil stencil,double* time_precon) {
763
764
765 double tau;
766 double delta;
767 double delta_prime=1.0;
768 double beta;
769
770
771 AdaptiveSparseGrid* grid = (x.getSparseGrid());
772 ListOfDepthOrderedSubgrids list(*grid);
773
774 if(!(x.getSparseGrid()->getKey()==f.getSparseGrid()->getKey()))exit(1);
775
776
777
778 VectorSparseG z(grid);
779 VectorSparseG r(grid);
780 VectorSparseG g(grid);
781 VectorSparseG d(grid);
782
783
784
785
786 VectorSparseG h(grid);
787 double j = 1.0;
788
789 // solves Ax = f;
790
791
792 // r = A*x - f;
793
794 // Start measuring time
795 struct timeval begin_all, end_all;
796 gettimeofday(&begin_all, 0);
797
798
799
800
801 // Start measuring time
802 struct timeval begin, end;
803 gettimeofday(&begin, 0);
804 //dirichlet_grid.completeDirichletGrid();
805 Preconditioning<Stencil> P(z,stencil);
806
807 // Stop measuring time and calculate the elapsed time
808 gettimeofday(&end, 0);
809 long seconds = end.tv_sec - begin.tv_sec;
810 long microseconds = end.tv_usec - begin.tv_usec;
811 double duration_def = seconds + microseconds*1e-6;
812 *time_precon = duration_def;
813
814
815 matrix.multiplication<Stencil>(x, r, stencil);
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830 // r = r - f = Ax - f;
831 r = f - r;
832
833 h = r;
834
835
836
837 P.apply_inverse(&h);
838
839
840 d = -1.0*h;
841
842
843
844
845 delta = product(r, h);
846
847
848
849
850 int k = 0;
851
852 for (int i = 1; i <= maxIteration && delta > eps; ++i) {
853 grid->WorkOnHangingNodes = false;
854 k = i;
855
856 g = 0.0;
857
858
859 matrix.multiplication<Stencil>(d, g,stencil);
860
861
862 tau = double(delta / product(d, g));
863
864 r = r + (tau * g);
865 x = x - (tau * d);
866
867 h = j * r;
868
869 P.apply_inverse(&h);
870 delta_prime = product(r, h);
871 if (delta_prime < eps*eps) break;
872
873
874 beta = delta_prime / delta;
875 delta = delta_prime;
876 d = -1.0*h + beta * d;
877
878
879 }
880
881
882 //cout << "CG finished after " << k << " iterations and residuum delta = " << sqrt(delta_prime) << endl;
883 *iterations = k;
884 count++;
885 return true;
886}
887
888#endif //GRUN_CG_METHOD_H
AdaptiveSparseGrid
Definition sparseGrid.h:277
ListOfDepthOrderedSubgrids
Definition ListOfDepthOrderedGrids.h:115
MatrixVectorHomogen
Definition MatrixVectorHomogen.h:32
MatrixVectorHomogen::multiplication
void multiplication(VectorSparseG &prew, VectorSparseG &Ax, StencilClass &stencilClass)
MultiDimCompass
Definition index.h:356
MultiDimFiveCompass
Definition index.h:436