petsc_solver.hpp 20.8 KB
Newer Older
Pietro Incardona's avatar
Pietro Incardona committed
1 2 3 4 5 6 7 8 9 10
/*
 * petsc_solver.hpp
 *
 *  Created on: Apr 26, 2016
 *      Author: i-bird
 */

#ifndef OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_HPP_
#define OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_HPP_

Pietro Incardona's avatar
Pietro Incardona committed
11 12
#ifdef HAVE_PETSC

Pietro Incardona's avatar
Pietro Incardona committed
13 14 15 16 17
#include "Vector/Vector.hpp"
#include "Eigen/UmfPackSupport"
#include <Eigen/SparseLU>
#include <petscksp.h>
#include <petsctime.h>
18
#include "Plot/GoogleChart.hpp"
Pietro Incardona's avatar
Pietro Incardona committed
19 20 21

#define UMFPACK_NONE 0

22

Pietro Incardona's avatar
Pietro Incardona committed
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
template<typename T>
class petsc_solver
{
public:

	template<typename impl> static Vector<T> solve(const SparseMatrix<T,impl> & A, const Vector<T> & b)
	{
		std::cerr << "Error Petsc only suppor double precision" << "/n";
	}
};

#define SOLVER_NOOPTION 0
#define SOLVER_PRINT_RESIDUAL_NORM_INFINITY 1
#define SOLVER_PRINT_DETERMINANT 2

38 39 40 41 42 43 44 45
size_t cpu_rank;

/*! \brief This class is able to do Matrix inversion in parallel with PETSC solvers
 *
 * #Example of how to solve a lid driven cavity 3D with a PETSC solver in paralle
 * \snippet
 *
 */
Pietro Incardona's avatar
Pietro Incardona committed
46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61
template<>
class petsc_solver<double>
{
	// It contain statistic of the error of the calculated solution
	struct solError
	{
		// infinity norm of the error
		PetscReal err_inf;

		// L1 norm of the error
		PetscReal err_norm;

		// Number of iterations
		PetscInt its;
	};

62 63 64 65 66 67 68 69 70 71 72 73 74
	//It contain statistic of the error at each iteration
	struct itError
	{
		PetscInt it;
		PetscReal err_inf;
		PetscReal err_norm;

		itError(const solError & e)
		{
			it = e.its;
			err_inf = e.err_inf;
			err_norm = e.err_norm;
		}
Pietro Incardona's avatar
Pietro Incardona committed
75 76 77

		// Default constructor
		itError()	{}
78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98
	};

	// It contain the benchmark
	struct solv_bench_info
	{
		// Method name
		std::string method;

		//Method name (short form)
		std::string smethod;

		// time to converge in milliseconds
		double time;

		// Solution error
		solError err;

		// Convergence per iteration
		openfpm::vector<itError> res;
	};

Pietro Incardona's avatar
Pietro Incardona committed
99
	// KSP Relative tolerance
Pietro Incardona's avatar
Pietro Incardona committed
100
//	PetscReal rtol;
Pietro Incardona's avatar
Pietro Incardona committed
101 102

	// KSP Absolute tolerance
Pietro Incardona's avatar
Pietro Incardona committed
103
//	PetscReal abstol;
Pietro Incardona's avatar
Pietro Incardona committed
104 105

	// KSP dtol tolerance to determine that the method diverge
Pietro Incardona's avatar
Pietro Incardona committed
106
//	PetscReal dtol;
Pietro Incardona's avatar
Pietro Incardona committed
107 108 109 110

	// KSP Maximum number of iterations
	PetscInt maxits;

111
	// Main parallel solver
Pietro Incardona's avatar
Pietro Incardona committed
112 113
	KSP ksp;

114 115 116
	// Temporal variable used for calculation of static members
	size_t tmp;

Pietro Incardona's avatar
Pietro Incardona committed
117
	// Solver used
Pietro Incardona's avatar
Pietro Incardona committed
118
//	char s_type[32];
Pietro Incardona's avatar
Pietro Incardona committed
119 120 121 122 123 124 125 126 127 128 129 130 131 132

	// Tollerance set by the solver
	PetscScalar tol;

	// The full set of solvers
	openfpm::vector<std::string> solvs;

	// The full set of preconditionses for simple parallel solvers

	// PCType pre-conditioner type
	PCType pc;

	bool try_solve = false;

133 134
	// It contain the solver benchmark results
	openfpm::vector<solv_bench_info> bench;
Pietro Incardona's avatar
Pietro Incardona committed
135

136
	/*! \brief Calculate the residual error at time t for one method
Pietro Incardona's avatar
Pietro Incardona committed
137
	 *
138 139 140 141
	 * \param t time
	 * \param slv solver
	 *
	 * \return the residual
Pietro Incardona's avatar
Pietro Incardona committed
142 143
	 *
	 */
144
	static double calculate_it(double t, solv_bench_info & slv)
Pietro Incardona's avatar
Pietro Incardona committed
145
	{
146
		double s_int = slv.time / slv.res.size();
Pietro Incardona's avatar
Pietro Incardona committed
147

148 149
		// Calculate the discrete point in time
		size_t pt = std::floor(t / s_int);
Pietro Incardona's avatar
Pietro Incardona committed
150

151 152 153 154
		if (pt < slv.res.size())
			return slv.res.get(pt).err_inf;

		return slv.res.last().err_inf;
Pietro Incardona's avatar
Pietro Incardona committed
155 156
	}

157
	/*! \brief Here we write the benchmark report
Pietro Incardona's avatar
Pietro Incardona committed
158 159 160
	 *
	 *
	 */
161
	void write_bench_report()
Pietro Incardona's avatar
Pietro Incardona committed
162
	{
163 164
		if (create_vcluster().getProcessUnitID() != 0)
			return;
Pietro Incardona's avatar
Pietro Incardona committed
165

166 167 168 169
		openfpm::vector<std::string> x;
		openfpm::vector<openfpm::vector<double>> y;
		openfpm::vector<std::string> yn;
		openfpm::vector<double> xd;
Pietro Incardona's avatar
Pietro Incardona committed
170

171 172
		for (size_t i = 0 ; i < bench.size() ; i++)
			x.add(bench.get(i).smethod);
Pietro Incardona's avatar
Pietro Incardona committed
173

174 175 176 177 178
		// Each colum can have multiple data set (in this case 4 dataset)
		// Each dataset can have a name
		yn.add("Norm infinity");
		yn.add("Norm average");
		yn.add("Number of iterations");
Pietro Incardona's avatar
Pietro Incardona committed
179

180
		// Each colums can have multiple data-set
Pietro Incardona's avatar
Pietro Incardona committed
181

182 183
		for (size_t i = 0 ; i < bench.size() ; i++)
			y.add({bench.get(i).err.err_inf,bench.get(i).err.err_norm,(double)bench.get(i).err.its});
Pietro Incardona's avatar
Pietro Incardona committed
184

185 186
		// Google charts options
		GCoptions options;
Pietro Incardona's avatar
Pietro Incardona committed
187

188 189 190 191 192
		options.title = std::string("Residual error");
		options.yAxis = std::string("Error");
		options.xAxis = std::string("Method");
		options.stype = std::string("bars");
		options.more = std::string("vAxis: {scaleType: 'log'}");
Pietro Incardona's avatar
Pietro Incardona committed
193

194
		GoogleChart cg;
Pietro Incardona's avatar
Pietro Incardona committed
195

196
		std::stringstream ss;
Pietro Incardona's avatar
Pietro Incardona committed
197

198 199 200
		cg.addHTML("<h2>Robustness</h2>");
		cg.AddHistGraph(x,y,yn,options);
		cg.addHTML("<h2>Speed</h2>");
Pietro Incardona's avatar
Pietro Incardona committed
201

202 203
		y.clear();
		yn.clear();
Pietro Incardona's avatar
Pietro Incardona committed
204

205 206
		yn.add("Time in millseconds");
		options.title = std::string("Time");
Pietro Incardona's avatar
Pietro Incardona committed
207

208 209
		for (size_t i = 0 ; i < bench.size() ; i++)
			y.add({bench.get(i).time});
Pietro Incardona's avatar
Pietro Incardona committed
210

211
		cg.AddHistGraph(x,y,yn,options);
Pietro Incardona's avatar
Pietro Incardona committed
212

213
		// Convergence in time
Pietro Incardona's avatar
Pietro Incardona committed
214

215 216 217 218
		x.clear();
		y.clear();
		yn.clear();
		xd.clear();
Pietro Incardona's avatar
Pietro Incardona committed
219

220
		// Get the maximum in time across all the solvers
Pietro Incardona's avatar
Pietro Incardona committed
221

222
		double max_time = 0.0;
Pietro Incardona's avatar
Pietro Incardona committed
223

224
		for (size_t i = 0 ; i < bench.size() ; i++)
Pietro Incardona's avatar
Pietro Incardona committed
225
		{
226 227
			if (bench.get(i).time > max_time)	{max_time = bench.get(i).time;}
			yn.add(bench.get(i).smethod);
Pietro Incardona's avatar
Pietro Incardona committed
228 229
		}

230
		size_t n_int = 300;
Pietro Incardona's avatar
Pietro Incardona committed
231

232
		// calculate dt
Pietro Incardona's avatar
Pietro Incardona committed
233

234
		double dt = max_time / n_int;
Pietro Incardona's avatar
Pietro Incardona committed
235

236
		// Add
Pietro Incardona's avatar
Pietro Incardona committed
237

238
		// For each solver we have a convergence plot
Pietro Incardona's avatar
Pietro Incardona committed
239

240 241 242 243 244 245 246
		for (double t = dt ; t <= max_time + 0.05 * max_time ; t += dt)
		{
			y.add();
			xd.add(t);
			for (size_t j = 0 ; j < bench.size() ; j++)
				y.last().add(calculate_it(t,bench.get(j)));
		}
Pietro Incardona's avatar
Pietro Incardona committed
247

248 249
		std::stringstream ss2;
		ss2 << "<h2>Convergence with time</h2><br>" << std::endl;
Pietro Incardona's avatar
Pietro Incardona committed
250

251 252 253 254 255
		options.title = std::string("Residual norm infinity");
		options.yAxis = std::string("residual");
		options.xAxis = std::string("Time in milliseconds");
		options.lineWidth = 1.0;
		options.more = std::string("vAxis: {scaleType: 'log'},hAxis: {scaleType: 'log'}");
Pietro Incardona's avatar
Pietro Incardona committed
256

257 258
		cg.addHTML(ss2.str());
		cg.AddLinesGraph(xd,y,yn,options);
Pietro Incardona's avatar
Pietro Incardona committed
259

260 261 262
		ss << "<h2>Solvers Tested</h2><br><br>" << std::endl;
		for (size_t i = 0 ; i < bench.size() ; i++)
			ss << bench.get(i).smethod << " = " << bench.get(i).method << "<br>" << std::endl;
Pietro Incardona's avatar
Pietro Incardona committed
263

264
		cg.addHTML(ss.str());
Pietro Incardona's avatar
Pietro Incardona committed
265

266 267
		cg.write("gc_solver.html");
	}
Pietro Incardona's avatar
Pietro Incardona committed
268

269 270 271 272 273 274 275
	/*! \brief Print a progress bar on standard out
	 *
	 *
	 */
	void print_progress_bar()
	{
		auto & v_cl = create_vcluster();
Pietro Incardona's avatar
Pietro Incardona committed
276

277 278 279
		if (v_cl.getProcessUnitID() == 0)
			std::cout << "-----------------------25-----------------------50-----------------------75----------------------100" << std::endl;
	}
Pietro Incardona's avatar
Pietro Incardona committed
280

281 282 283 284 285 286 287 288 289 290 291
	/*! \brief procedure to collect the absolute residue at each iteration
	 *
	 * \param ksp Solver
	 * \param it Iteration number
	 * \param res resudual
	 * \param custom pointer to data
	 *
	 */
	static PetscErrorCode monitor(KSP ksp,PetscInt it,PetscReal res,void* data)
	{
		petsc_solver<double> * pts = static_cast<petsc_solver<double> *>(data);
Pietro Incardona's avatar
Pietro Incardona committed
292

293 294 295 296
		Mat A;
		Mat P;
		Vec x;
		Vec b;
Pietro Incardona's avatar
Pietro Incardona committed
297

298 299 300
		PETSC_SAFE_CALL(KSPGetOperators(ksp,&A,&P));
		PETSC_SAFE_CALL(KSPGetSolution(ksp,&x));
		PETSC_SAFE_CALL(KSPGetRhs(ksp,&b));
Pietro Incardona's avatar
Pietro Incardona committed
301

302 303 304
		solError err = statSolutionError(A,b,x,ksp);
		itError erri(err);
		erri.it = it;
Pietro Incardona's avatar
Pietro Incardona committed
305

Pietro Incardona's avatar
Pietro Incardona committed
306 307 308 309 310 311 312 313 314 315
		//
		size_t old_size = pts->bench.last().res.size();
		pts->bench.last().res.resize(it+1);

		if (old_size > 0)
		{
			for (size_t i = old_size ; i < it-1 ; i++)
				pts->bench.last().res.get(i) = pts->bench.last().res.get(i-1);
		}

316 317
		// Add the error per iteration
		pts->bench.last().res.add(erri);
Pietro Incardona's avatar
Pietro Incardona committed
318

319
		pts->progress(it);
Pietro Incardona's avatar
Pietro Incardona committed
320

321
		return 0;
Pietro Incardona's avatar
Pietro Incardona committed
322 323
	}

324 325 326 327 328 329 330 331 332
	/*! \brief procedure print the progress of the solver in benchmark mode
	 *
	 * \param ksp Solver
	 * \param it Iteration number
	 * \param res resudual
	 * \param custom pointer to data
	 *
	 */
	static PetscErrorCode monitor_progress(KSP ksp,PetscInt it,PetscReal res,void* data)
Pietro Incardona's avatar
Pietro Incardona committed
333
	{
334
		petsc_solver<double> * pts = (petsc_solver *)data;
Pietro Incardona's avatar
Pietro Incardona committed
335

336
		pts->progress(it);
Pietro Incardona's avatar
Pietro Incardona committed
337

338 339
		return 0;
	}
Pietro Incardona's avatar
Pietro Incardona committed
340

341 342 343 344 345 346 347
	/*! \brief This function print an "*" showing the progress of the solvers
	 *
	 *
	 */
	void progress(PetscInt it)
	{
		PetscInt n_max_it;
Pietro Incardona's avatar
Pietro Incardona committed
348

349
		PETSC_SAFE_CALL(KSPGetTolerances(ksp,NULL,NULL,NULL,&n_max_it));
Pietro Incardona's avatar
Pietro Incardona committed
350

351
		auto & v_cl = create_vcluster();
Pietro Incardona's avatar
Pietro Incardona committed
352

353 354
		if (std::floor(it * 100.0 / n_max_it) > tmp)
		{
Pietro Incardona's avatar
Pietro Incardona committed
355
			size_t diff = it * 100.0 / n_max_it - tmp;
356 357
			tmp = it * 100.0 / n_max_it;
			if (v_cl.getProcessUnitID() == 0)
Pietro Incardona's avatar
Pietro Incardona committed
358 359 360 361
			{
				for (size_t k = 0 ; k <  diff ; k++)	{std::cout << "*";}
				std::cout << std::flush;
			}
362 363
		}
	}
Pietro Incardona's avatar
Pietro Incardona committed
364

365 366 367 368 369 370 371 372
	/*! \brief Print statistic about the solution error and method used
	 *
	 * \param err structure that contain the solution errors
	 *
	 */
	void print_stat(solError & err)
	{
		auto & v_cl = create_vcluster();
Pietro Incardona's avatar
Pietro Incardona committed
373

374 375
		if (v_cl.getProcessUnitID() == 0)
		{
Pietro Incardona's avatar
Pietro Incardona committed
376 377 378 379
			KSPType typ;
			KSPGetType(ksp,&typ);

			std::cout << "Method: " << typ << " " << " pre-conditoner: " << PCJACOBI << "  iterations: " << err.its << std::endl;
380 381 382
			std::cout << "Norm of error: " << err.err_norm << "   Norm infinity: " << err.err_inf << std::endl;
		}
	}
Pietro Incardona's avatar
Pietro Incardona committed
383

384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429
	/*! \brief It convert the KSP type into a human read-able string
	 *
	 *
	 */
	std::string to_string_method(const std::string & solv)
	{
		if (solv == std::string(KSPBCGS))
		{
			return std::string("BiCGStab (Stabilized version of BiConjugate Gradient Squared)");
		}
		else if (solv == std::string(KSPIBCGS))
		{
			return std::string("IBiCGStab (Improved Stabilized version of BiConjugate Gradient Squared)");
		}
		else if (solv == std::string(KSPFBCGS))
		{
			return std::string("FBiCGStab (Flexible Stabilized version of BiConjugate Gradient Squared)");
		}
		else if (solv == std::string(KSPFBCGSR))
		{
			return std::string("FBiCGStab-R (Modified Flexible Stabilized version of BiConjugate Gradient Squared)");
		}
		else if (solv == std::string(KSPBCGSL))
		{
			return std::string("BiCGStab(L) (Stabilized version of BiConjugate Gradient Squared with L search directions)");
		}
		else if (solv == std::string(KSPGMRES))
		{
			return std::string("GMRES(M) (Generalized Minimal Residual method with restart)");
		}
		else if (solv == std::string(KSPFGMRES))
		{
			return std::string("FGMRES(M) (Flexible Generalized Minimal Residual with restart)");
		}
		else if (solv == std::string(KSPLGMRES))
		{
			return std::string("LGMRES(M) (Augment Generalized Minimal Residual method with restart)");
		}
		else if (solv == std::string(KSPPGMRES))
		{
			return std::string("PGMRES(M) (Pipelined Generalized Minimal Residual method)");
		}
		else if (solv == std::string(KSPGCR))
		{
			return std::string("GCR (Generalized Conjugate Residual)");
		}
Pietro Incardona's avatar
Pietro Incardona committed
430

431 432
		return std::string("UNKNOWN");
	}
Pietro Incardona's avatar
Pietro Incardona committed
433

434 435 436 437 438 439 440 441 442 443 444 445
	/*! \brief Allocate a new benchmark slot for a method
	 *
	 * \param str Method name
	 *
	 */
	void new_bench(const std::string & str)
	{
		// Add a new benchmark result
		bench.add();
		bench.last().method = to_string_method(str);
		bench.last().smethod = std::string(str);
	}
Pietro Incardona's avatar
Pietro Incardona committed
446

447 448 449 450 451 452 453 454 455 456 457 458 459 460 461
	/*! \brief Copy the solution if better
	 *
	 * \param res Residual of the solution
	 * \param sol solution
	 * \param best_res residual of the best solution
	 * \param best_sol best solution
	 *
	 */
	void copy_if_better(double res, Vec & sol, double & best_res, Vec & best_sol)
	{
		if (res < best_res)
		{
			VecCopy(sol,best_sol);
			best_res = res;
		}
Pietro Incardona's avatar
Pietro Incardona committed
462 463
	}

464 465 466 467 468 469 470 471
	/*! \brief Try to solve the system x=inv(A)*b using all the Krylov solvers with simple Jacobi pre-conditioner
	 *
	 *   It try to solve the system using JACOBI pre-conditioner and all the Krylov solvers available at the end it write
	 *   a report
	 *
	 * \param A_ Matrix
	 * \param b_ vector of coefficents
	 * \param x solution
Pietro Incardona's avatar
Pietro Incardona committed
472 473
	 *
	 */
474
	void try_solve_simple(Mat & A_, const Vec & b_, Vec & x_)
Pietro Incardona's avatar
Pietro Incardona committed
475
	{
476 477
		Vec best_sol;
		PETSC_SAFE_CALL(VecDuplicate(x_,&best_sol));
Pietro Incardona's avatar
Pietro Incardona committed
478

479 480
		// Best residual
		double best_res = std::numeric_limits<double>::max();
Pietro Incardona's avatar
Pietro Incardona committed
481

482 483
		// Create a new VCluster
		auto & v_cl = create_vcluster();
Pietro Incardona's avatar
Pietro Incardona committed
484

485
		destroyKSP();
Pietro Incardona's avatar
Pietro Incardona committed
486

487 488 489 490 491
		// for each solver
		for (size_t i = 0 ; i < solvs.size() ; i++)
		{
			initKSPForTest();
			setSolver(solvs.get(i).c_str());
Pietro Incardona's avatar
Pietro Incardona committed
492

493 494 495 496 497 498 499
			// Setup for BCGSL, GMRES
			if (solvs.get(i) == std::string(KSPBCGSL))
			{
				// we try from 2 to 6 as search direction
				for (size_t j = 2 ; j < 6 ; j++)
				{
					new_bench(solvs.get(i));
Pietro Incardona's avatar
Pietro Incardona committed
500

501 502 503 504 505
					if (v_cl.getProcessUnitID() == 0)
						std::cout << "L = " << j << std::endl;
					bench.last().smethod += std::string("(") + std::to_string(j) + std::string(")");
					searchDirections(j);
					bench_solve_simple(A_,b_,x_,bench.last());
Pietro Incardona's avatar
Pietro Incardona committed
506

507 508 509 510 511 512 513 514 515 516 517 518
					copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
				}
			}
			else if (solvs.get(i) == std::string(KSPGMRES) ||
					 solvs.get(i) == std::string(std::string(KSPFGMRES)) ||
					 solvs.get(i) == std::string(std::string(KSPLGMRES)) )
			{
				// we try from 2 to 6 as search direction
				for (size_t j = 50 ; j < 300 ; j += 50)
				{
					// Add a new benchmark result
					new_bench(solvs.get(i));
Pietro Incardona's avatar
Pietro Incardona committed
519

520 521 522 523 524
					if (v_cl.getProcessUnitID() == 0)
						std::cout << "Restart = " << j << std::endl;
					bench.last().smethod += std::string("(") + std::to_string(j) + std::string(")");
					setRestart(j);
					bench_solve_simple(A_,b_,x_,bench.last());
Pietro Incardona's avatar
Pietro Incardona committed
525

526 527 528 529 530 531 532
					copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
				}
			}
			else
			{
				// Add a new benchmark result
				new_bench(solvs.get(i));
Pietro Incardona's avatar
Pietro Incardona committed
533

534
				bench_solve_simple(A_,b_,x_,bench.last());
Pietro Incardona's avatar
Pietro Incardona committed
535

536
				copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
Pietro Incardona's avatar
Pietro Incardona committed
537
			}
538 539 540 541 542 543 544 545

			destroyKSP();
		}

		write_bench_report();

		// Copy the best solution to x
		PETSC_SAFE_CALL(VecCopy(best_sol,x_));
Pietro Incardona's avatar
Pietro Incardona committed
546
		PETSC_SAFE_CALL(VecDestroy(&best_sol));
547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591
	}

	/*! \brief Benchmark solve simple solving x=inv(A)*b
	 *
	 * \param A_ Matrix A
	 * \param b_ vector b
	 * \param x_ solution x
	 * \param bench structure that store the benchmark information
	 *
	 */
	void bench_solve_simple(Mat & A_, const Vec & b_, Vec & x_, solv_bench_info & bench)
	{
		// timer for benchmark
		timer t;
		t.start();
		// Enable progress monitor
		tmp = 0;
		PETSC_SAFE_CALL(KSPMonitorCancel(ksp));
		PETSC_SAFE_CALL(KSPMonitorSet(ksp,monitor_progress,this,NULL));
		print_progress_bar();
		solve_simple(A_,b_,x_);

		// New line
		if (create_vcluster().getProcessUnitID() == 0)
			std::cout << std::endl;

		t.stop();
		bench.time = t.getwct() * 1000.0;

		// calculate error statistic about the solution and print
		solError err = statSolutionError(A_,b_,x_,ksp);
		print_stat(err);

		bench.err = err;

		// Solve getting the residual per iteration
		tmp = 0;
		PETSC_SAFE_CALL(KSPMonitorCancel(ksp));
		PETSC_SAFE_CALL(KSPMonitorSet(ksp,monitor,this,NULL));
		print_progress_bar();
		solve_simple(A_,b_,x_);

		// New line
		if (create_vcluster().getProcessUnitID() == 0)
			std::cout << std::endl;
Pietro Incardona's avatar
Pietro Incardona committed
592 593 594 595 596 597 598
	}

	/*! \brief solve simple use a Krylov solver + Simple selected Parallel Pre-conditioner
	 *
	 */
	void solve_simple(Mat & A_, const Vec & b_, Vec & x_)
	{
Pietro Incardona's avatar
Pietro Incardona committed
599
//		PETSC_SAFE_CALL(KSPSetType(ksp,s_type));
Pietro Incardona's avatar
Pietro Incardona committed
600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626

		// We set the size of x according to the Matrix A
		PetscInt row;
		PetscInt col;
		PetscInt row_loc;
		PetscInt col_loc;

		PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
		PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));

		PC pc;

		// We set the Matrix operators
		PETSC_SAFE_CALL(KSPSetOperators(ksp,A_,A_));

		// We set the pre-conditioner
		PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
		PETSC_SAFE_CALL(PCSetType(pc,PCJACOBI));

		// if we are on on best solve set-up a monitor function

		if (try_solve == true)
		{
			// Disable convergence check
			PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedSkip,NULL,NULL));
		}

627
		PETSC_SAFE_CALL(KSPSetFromOptions(ksp));
Pietro Incardona's avatar
Pietro Incardona committed
628 629 630 631 632 633 634 635 636 637 638 639

		// Solve the system
		PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
	}

	/*! \brief Calculate statistic on the error solution
	 *
	 * \brief A Matrix of the system
	 * \brief b Right hand side of the matrix
	 * \brief x Solution
	 *
	 */
640
	static solError statSolutionError(Mat & A_, const Vec & b_, Vec & x_, KSP ksp)
Pietro Incardona's avatar
Pietro Incardona committed
641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671
	{
		PetscScalar neg_one = -1.0;

		// We set the size of x according to the Matrix A
		PetscInt row;
		PetscInt col;
		PetscInt row_loc;
		PetscInt col_loc;
		PetscInt its;

		PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
		PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));

		/*
		      Here we calculate the residual error
	    */
		PetscReal norm;
		PetscReal norm_inf;

		// Get a vector r for the residual
		Vector<double,PETSC_BASE> r(row,row_loc);
		Vec & r_ = r.getVec();

		PETSC_SAFE_CALL(MatMult(A_,x_,r_));
		PETSC_SAFE_CALL(VecAXPY(r_,neg_one,b_));

		PETSC_SAFE_CALL(VecNorm(r_,NORM_1,&norm));
		PETSC_SAFE_CALL(VecNorm(r_,NORM_INFINITY,&norm_inf));
		PETSC_SAFE_CALL(KSPGetIterationNumber(ksp,&its));

		solError err;
672
		err.err_norm = norm / col;
Pietro Incardona's avatar
Pietro Incardona committed
673 674 675 676 677 678
		err.err_inf = norm_inf;
		err.its = its;

		return err;
	}

679 680 681 682 683 684 685
	/*! \brief initialize the KSP object
	 *
	 *
	 */
	void initKSP()
	{
		PETSC_SAFE_CALL(KSPCreate(PETSC_COMM_WORLD,&ksp));
Pietro Incardona's avatar
Pietro Incardona committed
686
//		PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
687 688 689 690 691 692 693 694 695
	}

	/*! \brief initialize the KSP object for solver testing
	 *
	 *
	 */
	void initKSPForTest()
	{
		PETSC_SAFE_CALL(KSPCreate(PETSC_COMM_WORLD,&ksp));
Pietro Incardona's avatar
Pietro Incardona committed
696 697 698 699
//		PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
//		PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol,dtol,300));

		setMaxIter(300);
700 701 702 703 704 705 706 707 708 709 710 711 712

		// Disable convergence check
		PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedSkip,NULL,NULL));
	}

	/*! \brief Destroy the KSP object
	 *
	 *
	 */
	void destroyKSP()
	{
		KSPDestroy(&ksp);
	}
Pietro Incardona's avatar
Pietro Incardona committed
713 714 715

public:

Pietro Incardona's avatar
Pietro Incardona committed
716 717 718 719
	~petsc_solver()
	{
	}

Pietro Incardona's avatar
Pietro Incardona committed
720 721
	petsc_solver()
	{
722
		initKSP();
Pietro Incardona's avatar
Pietro Incardona committed
723 724 725

		// Add the solvers

726 727 728 729 730
		solvs.add(std::string(KSPBCGS));
//		solvs.add(std::string(KSPIBCGS)); <--- Produce invalid argument
		solvs.add(std::string(KSPFBCGS));
//		solvs.add(std::string(KSPFBCGSR)); <--- Nan production problems
		solvs.add(std::string(KSPBCGSL));
Pietro Incardona's avatar
Pietro Incardona committed
731
		solvs.add(std::string(KSPGMRES));
732 733 734
		solvs.add(std::string(KSPFGMRES));
		solvs.add(std::string(KSPLGMRES));
//		solvs.add(std::string(KSPPGMRES)); <--- Take forever
Pietro Incardona's avatar
Pietro Incardona committed
735
//		solvs.add(std::string(KSPGCR));
Pietro Incardona's avatar
Pietro Incardona committed
736 737

		setSolver(KSPGMRES);
Pietro Incardona's avatar
Pietro Incardona committed
738 739
	}

740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767
	/*! \brief Add a test solver
	 *
	 * The try solve function use the most robust solvers in PETSC, if you want to add
	 * additionally solver like KSPIBCGS,KSPFBCGSR,KSPPGMRES, use addTestSolver(std::string(KSPIBCGS))
	 *
	 */
	void addTestSolver(std::string & solver)
	{
		solvs.add(solver);
	}

	/*! \brief Remove a test solver
	 *
	 * The try solve function use the most robust solvers in PETSC, if you want to remove
	 * a solver like use removeTestSolver(std::string(KSPIBCGS))
	 *
	 */
	void removeTestSolver(const std::string & solver)
	{
		// search for the test solver and remove it

		for (size_t i = 0 ; i < solvs.size() ; i++)
		{
			if (solvs.get(i) == solver)
				return;
		}
	}

Pietro Incardona's avatar
Pietro Incardona committed
768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787
	/*! \brief Set the solver to test all the solvers and generate a report
	 *
	 * In this mode the system will try different Solvers, Preconditioner and
	 * combination of solvers in order to find the best solver in speed, and
	 * precision. As output it will produce a performance report
	 *
	 *
	 */
	void best_solve()
	{
		try_solve = true;
	}

	/*! \brief Set the Petsc solver
	 *
	 * \see KSPType in PETSC manual for a list of all PETSC solvers
	 *
	 */
	void setSolver(KSPType type)
	{
Pietro Incardona's avatar
Pietro Incardona committed
788
		PetscOptionsSetValue("-ksp_type",type);
Pietro Incardona's avatar
Pietro Incardona committed
789 790 791 792 793 794 795 796 797 798 799
	}

	/*! \brief Set the relative tolerance as stop criteria
	 *
	 * \see PETSC manual KSPSetTolerances for an explanation
	 *
	 * \param rtol Relative tolerance
	 *
	 */
	void setRelTol(PetscReal rtol_)
	{
Pietro Incardona's avatar
Pietro Incardona committed
800
		PetscOptionsSetValue("-ksp_rtol",std::to_string(rtol_).c_str());
Pietro Incardona's avatar
Pietro Incardona committed
801 802 803 804 805 806 807 808 809 810 811
	}

	/*! \brief Set the absolute tolerance as stop criteria
	 *
	 * \see PETSC manual KSPSetTolerances for an explanation
	 *
	 * \param abstol Absolute tolerance
	 *
	 */
	void setAbsTol(PetscReal abstol_)
	{
Pietro Incardona's avatar
Pietro Incardona committed
812
		PetscOptionsSetValue("-ksp_atol",std::to_string(abstol_).c_str());
Pietro Incardona's avatar
Pietro Incardona committed
813 814 815 816 817 818 819 820 821 822 823
	}

	/*! \brief Set the divergence tolerance
	 *
	 * \see PETSC manual KSPSetTolerances for an explanation
	 *
	 * \param divtol
	 *
	 */
	void setDivTol(PetscReal dtol_)
	{
Pietro Incardona's avatar
Pietro Incardona committed
824 825 826 827 828 829 830 831 832 833
		PetscOptionsSetValue("-ksp_dtol",std::to_string(dtol_).c_str());
	}

	/*! \brief Set the maximum number of iteration for Krylov solvers
	 *
	 *
	 */
	void setMaxIter(PetscInt n)
	{
		PetscOptionsSetValue("-ksp_max_it",std::to_string(n).c_str());
Pietro Incardona's avatar
Pietro Incardona committed
834 835 836 837 838 839 840 841
	}

	/*! For the BiCGStab(L) it define the number of search directions
	 *
	 * BiCG Methods can fail for base break-down (or near break-down). BiCGStab(L) try
	 * to avoid this problem. Such method has a parameter L (2 by default) Bigger is L
	 * more the system will try to avoid the breakdown
	 *
842
	 * \param l Increasing L should reduce the probability of failure of the solver because of break-down of the base
Pietro Incardona's avatar
Pietro Incardona committed
843 844 845 846
	 *
	 */
	void searchDirections(PetscInt l)
	{
847 848 849 850 851 852 853 854 855 856 857
		PetscOptionsSetValue("-ksp_bcgsl_ell",std::to_string(l).c_str());
	}

	/*! \brief For GMRES based method,  the number of Krylov directions to orthogonalize against
	 *
	 * \param n number of directions
	 *
	 */
	void setRestart(PetscInt n)
	{
		PetscOptionsSetValue("-ksp_gmres_restart",std::to_string(n).c_str());
Pietro Incardona's avatar
Pietro Incardona committed
858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895
	}

	/*! \brief Here we invert the matrix and solve the system
	 *
	 *  \warning umfpack is not a parallel solver, this function work only with one processor
	 *
	 *  \note if you want to use umfpack in a NON parallel, but on a distributed data, use solve with triplet
	 *
	 *	\tparam impl Implementation of the SparseMatrix
	 *
	 */
	Vector<double,PETSC_BASE> solve(SparseMatrix<double,int,PETSC_BASE> & A, const Vector<double,PETSC_BASE> & b)
	{
		Mat & A_ = A.getMat();
		const Vec & b_ = b.getVec();

		// We set the size of x according to the Matrix A
		PetscInt row;
		PetscInt col;
		PetscInt row_loc;
		PetscInt col_loc;

		PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
		PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));

		Vector<double,PETSC_BASE> x(row,row_loc);
		Vec & x_ = x.getVec();

		if (try_solve == true)
			try_solve_simple(A_,b_,x_);
		else
			solve_simple(A_,b_,x_);

		x.update();
		return x;
	}
};

Pietro Incardona's avatar
Pietro Incardona committed
896
#endif
Pietro Incardona's avatar
Pietro Incardona committed
897 898

#endif /* OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_HPP_ */