OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
petsc_solver_AMG_report.hpp
1/*
2 * petsc_solver_report.hpp
3 *
4 * Created on: Jun 22, 2017
5 * Author: i-bird
6 */
7
8#ifndef OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_AMG_REPORT_HPP_
9#define OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_AMG_REPORT_HPP_
10
11#include "config.h"
12
13#ifdef HAVE_PETSC
14
15#include <fstream>
16#include "Solvers/petsc_solver.hpp"
17#include "Vector/map_vector.hpp"
18
24{
26 double time_setup;
27
29 double time_solve;
30
32 double residual;
33
35 double score;
36};
37
43{
45 size_t num_sweep_test = 10;
46
48 double target_accuracy = 1;
49
52
55
58
61
64
67
70
81 petsc_solver<double> & solver,
83 {
84 timer tm_solve;
85 tm_solve.start();
86 auto x_ = solver.solve(A,b);
87 tm_solve.stop();
88
89 solError serr = solver.get_residual_error(A,x_,b);
90
91 // In order to measure the time to solve we solve again the system
92
93 timer tm_solve2;
94 tm_solve2.start();
95 solver.solve(x_,b);
96 tm_solve2.stop();
97
98 double time1 = tm_solve.getwct();
99 double time2 = tm_solve2.getwct();
100
101 Vcluster<> & v_cl = create_vcluster();
102 v_cl.max(time1);
103 v_cl.max(time2);
104 v_cl.execute();
105
106 // Save the result
108 tmp.time_setup = time1 - time2;
109 tmp.time_solve = time2;
110 tmp.residual = serr.err_inf;
111
112 perf_amg.add(tmp);
113
114 if (v_cl.getProcessUnitID() == 0)
115 std::cout << "Time1: " << time1 << " Time2: " << time2 << std::endl;
116 }
117
125 {
126 Vcluster<> & v_cl = create_vcluster();
127
129
131 solver.setSolver(KSPRICHARDSON);
132 solver.setAbsTol(0.01);
133 solver.setMaxIter(1);
134
136
137 solver.setPreconditioner(PCHYPRE_BOOMERAMG);
138 solver.setPreconditionerAMG_nl(6);
139 solver.setPreconditionerAMG_maxit(3);
140 solver.setPreconditionerAMG_relax("SOR/Jacobi");
141 solver.setPreconditionerAMG_cycleType("V",6,6);
142 solver.setPreconditionerAMG_coarsenNodalType(0);
143 solver.log_monitor();
144
145
146 // Reset the Matrix A
147 A.getMatrixTriplets();
148 if (v_cl.getProcessUnitID() == 0) {std::cout << "Benchmarking HMIS coarsener" << std::endl;}
149 solver.setPreconditionerAMG_coarsen("HMIS");
150 benchmark(A,b,solver,perf_amg_coars);
151 method.add("HMIS");
152
153 // Reset the Matrix A
154 A.getMatrixTriplets();
155 if (v_cl.getProcessUnitID() == 0) {std::cout << "Benchmarking PMIS coarsener" << std::endl;}
156 solver.setPreconditionerAMG_coarsen("PMIS");
157 benchmark(A,b,solver,perf_amg_coars);
158 method.add("PMIS");
159
160 // Reset the Matrix A
161 A.getMatrixTriplets();
162 if (v_cl.getProcessUnitID() == 0) {std::cout << "Benchmarking Falgout coarsener" << std::endl;}
163 solver.setPreconditionerAMG_coarsen("Falgout");
164 benchmark(A,b,solver,perf_amg_coars);
165 method.add("Falgout");
166
167 // Reset the Matrix A
168 A.getMatrixTriplets();
169 if (v_cl.getProcessUnitID() == 0) {std::cout << "Benchmarking modifiedRuge-Stueben coarsener" << std::endl;}
170 solver.setPreconditionerAMG_coarsen("modifiedRuge-Stueben");
171 benchmark(A,b,solver,perf_amg_coars);
172 method.add("modifiedRuge-Stueben");
173
174 // Reset the Matrix A
175 A.getMatrixTriplets();
176 if (v_cl.getProcessUnitID() == 0) {std::cout << "Benchmarking Ruge-Stueben coarsener" << std::endl;}
177 solver.setPreconditionerAMG_coarsen("Ruge-Stueben");
178 benchmark(A,b,solver,perf_amg_coars);
179 method.add("Ruge-Stueben");
180
181 // Reset the Matrix A
182 A.getMatrixTriplets();
183 if (v_cl.getProcessUnitID() == 0) {std::cout << "Benchmarking CLJP coarsener" << std::endl;}
184 solver.setPreconditionerAMG_coarsen("CLJP");
185 benchmark(A,b,solver,perf_amg_coars);
186 method.add("CLJP");
187
188 }
189
197 double score_solver(double t_solve, double t_m)
198 {
199 return 1.0/t_solve * std::min(1.0,target_accuracy/t_m);
200 }
201
208 {
209 std::cout << "Composing coarsening report" << std::endl;
210
215
216 for (size_t i = 0 ; i < perf_amg_coars.size() ; i++)
217 x.add(method.get(i));
218
219 // Each colum can have multiple data set (in this case 4 dataset)
220 // Each dataset can have a name
221 yn.add("Norm infinity");
222 yn.add("time to setup");
223 yn.add("time to solve");
224
225 // Each colums can have multiple data-set
226 for (size_t i = 0 ; i < perf_amg_coars.size() ; i++)
227 y.add({perf_amg_coars.get(i).residual,perf_amg_coars.get(i).time_setup,perf_amg_coars.get(i).time_solve});
228
229 // Google charts options
230 GCoptions options;
231
232 options.title = std::string("Overview of different coarsener on a V cycle with 6 relaxation SOR/Jacobi sweeps for each level up and down");
233 options.yAxis = std::string("Error/time");
234 options.xAxis = std::string("Method");
235 options.stype = std::string("bars");
236 options.more = std::string("vAxis: {scaleType: 'log'}");
237
238 std::stringstream ss;
239
240 cg.addHTML("<h2>Coarsener</h2>");
241 cg.AddHistGraph(x,y,yn,options);
242 }
243
251 {
252 cg.addHTML("<h2>V cycle sweep</h2>");
253 for(size_t k = 0 ; k < perf_amg_sweep_sym.size() / num_sweep_test ; k++)
254 {
259
260 x.clear();
261 for (size_t i = 0 ; i < num_sweep_test ; i++)
262 x.add(v_cycle_tested_sym.get(i));
263
264 // Each colum can have multiple data set (in this case 4 dataset)
265 // Each dataset can have a name
266 yn.add("error");
267 yn.add("time");
268 yn.add("score");
269
270 // Each colums can have multiple data-set
271 for (size_t i = 0 ; i < num_sweep_test ; i++)
272 {
273 double score = score_solver(perf_amg_sweep_sym.get(k*num_sweep_test+i).time_solve,perf_amg_sweep_sym.get(k*num_sweep_test+i).residual);
274 y.add({perf_amg_sweep_sym.get(k*num_sweep_test+i).residual,perf_amg_sweep_sym.get(k*num_sweep_test+i).time_solve,score});
275 }
276
277 // Google charts options
278 GCoptions options;
279
280 options.title = std::string("Coarsener: " + coars.get(k)) + std::string(": Overview of V cycle with n sweeps symmetrically up and down");
281 options.yAxis = std::string("error/time/score");
282 options.xAxis = std::string("sweeps");
283 options.stype = std::string("bars");
284 options.more = std::string("vAxis: {scaleType: 'log'}");
285
286 std::stringstream ss;
287
288 cg.AddHistGraph(x,y,yn,options);
289 }
290 }
291
299 {
300 cg.addHTML("<h2>V cycle sweep asymmetric test</h2>");
301 for(size_t k = 0 ; k < coars.size() ; k++)
302 {
307
308 size_t num_sweep_test = asym_tests.get(k+1) - asym_tests.get(k);
309 size_t sweep_start = asym_tests.get(k);
310
311 x.clear();
312 for (size_t i = 0 ; i < num_sweep_test ; i++)
313 {x.add(std::string("(" + std::to_string(v_cycle_tested_asym.get(i+sweep_start).first) + "," + std::to_string(v_cycle_tested_asym.get(i+sweep_start).second) + ")"));}
314
315 // Each colum can have multiple data set (in this case 4 dataset)
316 // Each dataset can have a name
317 yn.add("error");
318 yn.add("time");
319 yn.add("score");
320
321 // Each colums can have multiple data-set
322 for (size_t i = 0 ; i < num_sweep_test ; i++)
323 {
324 double score = score_solver(perf_amg_sweep_asym.get(sweep_start+i).time_solve,perf_amg_sweep_asym.get(sweep_start+i).residual);
325 y.add({perf_amg_sweep_asym.get(sweep_start+i).residual,perf_amg_sweep_asym.get(sweep_start+i).time_solve,score});
326 }
327
328 // Google charts options
329 GCoptions options;
330
331 options.title = std::string("Coarsener: " + coars.get(k)) + std::string(": Overview of V cycle with n sweeps asymmetrically distributed");
332 options.yAxis = std::string("error/time/score");
333 options.xAxis = std::string("sweeps");
334 options.stype = std::string("bars");
335 options.more = std::string("vAxis: {scaleType: 'log'}");
336
337 std::stringstream ss;
338
339 cg.AddHistGraph(x,y,yn,options);
340 }
341 }
342
352 openfpm::vector<std::string> & coarsener_to_test)
353 {
354 Vcluster<> & v_cl = create_vcluster();
355
357
359 solver.setSolver(KSPRICHARDSON);
360 solver.setAbsTol(0.01);
361 solver.setMaxIter(1);
362
365
366 for (size_t k = 0 ; k < coarsener_to_test.size(); k++)
367 {
368 for (size_t i = 1 ; i <= num_sweep_test ; i++)
369 {
370
371 solver.setPreconditioner(PCHYPRE_BOOMERAMG);
372 solver.setPreconditionerAMG_nl(6);
373 solver.setPreconditionerAMG_maxit(3);
374 solver.setPreconditionerAMG_relax("SOR/Jacobi");
375 solver.setPreconditionerAMG_cycleType("V",i,i);
376 solver.setPreconditionerAMG_coarsenNodalType(0);
377 solver.setPreconditionerAMG_coarsen(coarsener_to_test.get(k));
378 solver.log_monitor();
379
380 A.getMatrixTriplets();
381 solver.setPreconditionerAMG_coarsen(coarsener_to_test.get(k).c_str());
382 benchmark(A,b,solver,perf_amg_sweep_sym);
383
384 double score = score_solver(perf_amg_sweep_sym.get(k*num_sweep_test+(i-1)).time_solve,perf_amg_sweep_sym.get(k*num_sweep_test+(i-1)).residual);
385 perf_amg_sweep_sym.last().score = score;
386
387 if (v_cl.getProcessUnitID() == 0)
388 std::cout << "Coarsener: " << coarsener_to_test.get(k) << " Sweep: " << i << " Score: " << score << std::endl;
389
390 v_cycle_tested_sym.add(i);
391 }
392 }
393
394 }
395
407 openfpm::vector<std::string> & coarsener_to_test)
408 {
409 Vcluster<> & v_cl = create_vcluster();
410
412
414 solver.setSolver(KSPRICHARDSON);
415 solver.setAbsTol(0.01);
416 solver.setMaxIter(1);
417
420
421 for (size_t k = 0 ; k < coarsener_to_test.size(); k++)
422 {
424 size_t num_tot_sweep = 2*ids_ts.get(k);
425 for (size_t i = 0 ; i <= num_tot_sweep ; i++)
426 {
427
428 solver.setPreconditioner(PCHYPRE_BOOMERAMG);
429 solver.setPreconditionerAMG_nl(6);
430 solver.setPreconditionerAMG_maxit(3);
431 solver.setPreconditionerAMG_relax("SOR/Jacobi");
432 solver.setPreconditionerAMG_cycleType("V",i,num_tot_sweep-i);
433 solver.setPreconditionerAMG_coarsenNodalType(0);
434 solver.setPreconditionerAMG_coarsen(coarsener_to_test.get(k));
435 solver.log_monitor();
436
437 A.getMatrixTriplets();
438 solver.setPreconditionerAMG_coarsen(coarsener_to_test.get(k).c_str());
439
440 benchmark(A,b,solver,perf_amg_sweep_asym);
441
442 double score = score_solver(perf_amg_sweep_asym.last().time_solve,perf_amg_sweep_asym.last().residual);
443
444 if (v_cl.getProcessUnitID() == 0)
445 std::cout << "Coarsener: " << coarsener_to_test.get(k) << " Sweep: (" << i << "," << num_tot_sweep-i << ") Score: " << score << std::endl;
446
447 v_cycle_tested_asym.add(std::pair<size_t,size_t>(i,num_tot_sweep-i));
448 }
449
450 }
451
453 }
454
461 {
462 if (create_vcluster().getProcessUnitID() != 0)
463 return;
464
465 GoogleChart cg;
467
468 // Write V cycles performances
469 write_report_cycle(cg,coars);
470
471 write_report_cycle_asym(cg,coars);
472
473 cg.write("gc_AMG_preconditioners.html");
474 }
475
482 {
483 int fastest = 0;
484 double best_time = perf_amg_coars.get(0).time_solve;
485 int accurate = 0;
486 double best_accuracy = perf_amg_coars.get(0).residual;
487
488 for (size_t i = 1 ; i < perf_amg_coars.size() ; i++)
489 {
490 if (perf_amg_coars.get(i).time_solve < best_time)
491 {
492 fastest = i;
493 best_time = perf_amg_coars.get(i).time_solve;
494 }
495
496 if (perf_amg_coars.get(i).residual < best_accuracy)
497 {
498 accurate = i;
499 best_accuracy = perf_amg_coars.get(i).residual;
500 }
501 }
502
503 coars.add(method.get(fastest));
504 coars.add(method.get(accurate));
505 }
506
515 openfpm::vector<size_t> & sw_optimal,
516 size_t nm)
517 {
518 size_t page = perf.size() / nm;
519
520 for (size_t k = 0 ; k < nm ; k++)
521 {
522 int score_id = page*k;
523 double best_score = perf.get(score_id).score;
524
525 for (size_t i = page*k ; i < page*k+page ; i++)
526 {
527 if (perf.get(i).score > best_score)
528 {
529 score_id = i;
530 best_score = perf.get(i).score;
531 }
532 }
533
534 sw_optimal.add(v_cycle_tested_sym.get(score_id));
535 }
536 }
537
538public:
539
550 void setTergetAMGAccuracy(double t_a)
551 {
552 target_accuracy = t_a;
553 }
554
562 void setMaxSweep(int max_sw)
563 {
564 num_sweep_test = max_sw;
565 }
566
576 {
577 test_coarsener(A,b);
578
580 best_coarsener(coa_methods);
581
582 test_cycle_type(A,b,coa_methods);
583 openfpm::vector<size_t> sweep_optimal;
584 best_score(perf_amg_sweep_sym,sweep_optimal,coa_methods.size());
585
586 test_cycle_asym(A,b,sweep_optimal,coa_methods);
587
588 write_report(coa_methods);
589 }
590};
591
592#endif
593
594#endif /* OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_AMG_REPORT_HPP_ */
Small class to produce graph with Google chart in HTML.
void write(std::string file)
It write the graphs on file in html format using Google charts.
void AddHistGraph(openfpm::vector< Y > &y)
Add an histogram graph.
void addHTML(const std::string &html)
Add HTML text.
Sparse Matrix implementation.
void execute()
Execute all the requests.
size_t getProcessUnitID()
Get the process unit id.
void max(T &num)
Get the maximum number across all processors (or reduction with infinity norm)
Implementation of VCluster class.
Definition VCluster.hpp:59
Sparse Matrix implementation stub object when OpenFPM is compiled with no linear algebra support.
Definition Vector.hpp:40
Implementation of 1-D std::vector like structure.
size_t size()
Stub size.
Class to test AMG solvers.
void test_cycle_asym(SparseMatrix< double, int, PETSC_BASE > &A, const Vector< double, PETSC_BASE > &b, openfpm::vector< size_t > &ids_ts, openfpm::vector< std::string > &coarsener_to_test)
test best asymmetric combination of sweeps
void test_coarsener(SparseMatrix< double, int, PETSC_BASE > &A, const Vector< double, PETSC_BASE > &b)
test the corasener for this problem
double target_accuracy
Target accuracy for the preconditioner.
void setMaxSweep(int max_sw)
Set the target accuracy to score the AMG solver.
double score_solver(double t_solve, double t_m)
Score the solver.
void benchmark(SparseMatrix< double, int, PETSC_BASE > &A, const Vector< double, PETSC_BASE > &b, petsc_solver< double > &solver, openfpm::vector< AMG_time_err_coars > &perf_amg)
benchmark the selected setting for the preconditioner
openfpm::vector< AMG_time_err_coars > perf_amg_sweep_sym
It contain the performance of several AMG methods sweep configuration.
openfpm::vector< std::string > method
List of methods.
openfpm::vector< std::pair< size_t, size_t > > v_cycle_tested_asym
List of tested cycles in case of asymmetric.
openfpm::vector< size_t > v_cycle_tested_sym
List of tested cycles.
void test_cycle_type(SparseMatrix< double, int, PETSC_BASE > &A, const Vector< double, PETSC_BASE > &b, openfpm::vector< std::string > &coarsener_to_test)
test the corasener for this problem
openfpm::vector< size_t > asym_tests
starting point of each coarsening test for the asymmetric test
void write_report(openfpm::vector< std::string > &coars)
void write_report_cycle(GoogleChart &cg, openfpm::vector< std::string > &coars)
Write the report for coarsening.
size_t num_sweep_test
Number of sweeps for test.
void best_coarsener(openfpm::vector< std::string > &coars)
take the most accurate and the fastest AMG
void best_score(openfpm::vector< AMG_time_err_coars > &perf, openfpm::vector< size_t > &sw_optimal, size_t nm)
Return the best scoring solver.
openfpm::vector< AMG_time_err_coars > perf_amg_coars
It contains the performance of several AMG methods coarsener.
void try_solve(SparseMatrix< double, int, PETSC_BASE > &A, const Vector< double, PETSC_BASE > &b)
Try to use AMG pre-conditioner and check how they they perform.
void setTergetAMGAccuracy(double t_a)
Set the target accuracy to score the AMG solver.
void write_report_cycle_asym(GoogleChart &cg, openfpm::vector< std::string > &coars)
Write the report for coarsening.
void write_report_coars(GoogleChart &cg)
Write the report for coarsening.
openfpm::vector< AMG_time_err_coars > perf_amg_sweep_asym
It contain the performance of several AMG methods sweep configuration.
In case T does not match the PETSC precision compilation create a stub structure.
static Vector< T > solve(const SparseMatrix< T, impl > &A, const Vector< T > &b)
Solve the linear system.
Class for cpu time benchmarking.
Definition timer.hpp:28
void stop()
Stop the timer.
Definition timer.hpp:119
void start()
Start the timer.
Definition timer.hpp:90
double getwct()
Return the elapsed real time.
Definition timer.hpp:130
It contain information about the performance of the AMG.
double score
assign a score to the solver
double time_solve
time to solve
double residual
residual norm
double time_setup
time to setup
Google chart options.
std::string xAxis
X axis name.
std::string more
more
std::string title
Title of the chart.
std::string stype
std::string yAxis
Y axis name.
It contain statistic of the error of the calculated solution.
PetscReal err_inf
infinity norm of the error