petsc_solver.hpp 21.9 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_

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"
19 20
#include "Matrix/SparseMatrix.hpp"
#include "Vector/Vector.hpp"
Pietro Incardona's avatar
Pietro Incardona committed
21 22 23

#define UMFPACK_NONE 0

24

Pietro Incardona's avatar
Pietro Incardona committed
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
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

40 41 42 43 44 45 46 47
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
48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
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;
	};

64 65 66 67 68 69 70 71 72 73 74 75 76
	//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
77 78 79

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

	// 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
101 102 103
	// KSP Maximum number of iterations
	PetscInt maxits;

104
	// Main parallel solver
Pietro Incardona's avatar
Pietro Incardona committed
105 106
	KSP ksp;

107 108 109
	// Temporal variable used for calculation of static members
	size_t tmp;

Pietro Incardona's avatar
Pietro Incardona committed
110 111 112 113 114
	// The full set of solvers
	openfpm::vector<std::string> solvs;

	bool try_solve = false;

115 116
	// It contain the solver benchmark results
	openfpm::vector<solv_bench_info> bench;
Pietro Incardona's avatar
Pietro Incardona committed
117

118
	/*! \brief Calculate the residual error at time t for one method
Pietro Incardona's avatar
Pietro Incardona committed
119
	 *
120 121 122 123
	 * \param t time
	 * \param slv solver
	 *
	 * \return the residual
Pietro Incardona's avatar
Pietro Incardona committed
124 125
	 *
	 */
126
	static double calculate_it(double t, solv_bench_info & slv)
Pietro Incardona's avatar
Pietro Incardona committed
127
	{
128
		double s_int = slv.time / slv.res.size();
Pietro Incardona's avatar
Pietro Incardona committed
129

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

133 134 135 136
		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
137 138
	}

139
	/*! \brief Here we write the benchmark report
Pietro Incardona's avatar
Pietro Incardona committed
140 141 142
	 *
	 *
	 */
143
	void write_bench_report()
Pietro Incardona's avatar
Pietro Incardona committed
144
	{
145 146
		if (create_vcluster().getProcessUnitID() != 0)
			return;
Pietro Incardona's avatar
Pietro Incardona committed
147

148 149 150 151
		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
152

153 154
		for (size_t i = 0 ; i < bench.size() ; i++)
			x.add(bench.get(i).smethod);
Pietro Incardona's avatar
Pietro Incardona committed
155

156 157 158 159 160
		// 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
161

162
		// Each colums can have multiple data-set
Pietro Incardona's avatar
Pietro Incardona committed
163

164 165
		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
166

167 168
		// Google charts options
		GCoptions options;
Pietro Incardona's avatar
Pietro Incardona committed
169

170 171 172 173 174
		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
175

176
		GoogleChart cg;
Pietro Incardona's avatar
Pietro Incardona committed
177

178
		std::stringstream ss;
Pietro Incardona's avatar
Pietro Incardona committed
179

180 181 182
		cg.addHTML("<h2>Robustness</h2>");
		cg.AddHistGraph(x,y,yn,options);
		cg.addHTML("<h2>Speed</h2>");
Pietro Incardona's avatar
Pietro Incardona committed
183

184 185
		y.clear();
		yn.clear();
Pietro Incardona's avatar
Pietro Incardona committed
186

187 188
		yn.add("Time in millseconds");
		options.title = std::string("Time");
Pietro Incardona's avatar
Pietro Incardona committed
189

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

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

195
		// Convergence in time
Pietro Incardona's avatar
Pietro Incardona committed
196

197 198 199 200
		x.clear();
		y.clear();
		yn.clear();
		xd.clear();
Pietro Incardona's avatar
Pietro Incardona committed
201

202
		// Get the maximum in time across all the solvers
Pietro Incardona's avatar
Pietro Incardona committed
203

204
		double max_time = 0.0;
Pietro Incardona's avatar
Pietro Incardona committed
205

206
		for (size_t i = 0 ; i < bench.size() ; i++)
Pietro Incardona's avatar
Pietro Incardona committed
207
		{
208 209
			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
210 211
		}

212
		size_t n_int = maxits;
Pietro Incardona's avatar
Pietro Incardona committed
213

214
		// calculate dt
Pietro Incardona's avatar
Pietro Incardona committed
215

216
		double dt = max_time / n_int;
Pietro Incardona's avatar
Pietro Incardona committed
217

218
		// Add
Pietro Incardona's avatar
Pietro Incardona committed
219

220
		// For each solver we have a convergence plot
Pietro Incardona's avatar
Pietro Incardona committed
221

222 223 224 225 226 227 228
		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
229

230 231
		std::stringstream ss2;
		ss2 << "<h2>Convergence with time</h2><br>" << std::endl;
Pietro Incardona's avatar
Pietro Incardona committed
232

233 234 235 236 237
		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
238

239 240
		cg.addHTML(ss2.str());
		cg.AddLinesGraph(xd,y,yn,options);
Pietro Incardona's avatar
Pietro Incardona committed
241

242 243 244
		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
245

246
		cg.addHTML(ss.str());
Pietro Incardona's avatar
Pietro Incardona committed
247

248 249
		cg.write("gc_solver.html");
	}
Pietro Incardona's avatar
Pietro Incardona committed
250

251 252 253 254 255 256 257
	/*! \brief Print a progress bar on standard out
	 *
	 *
	 */
	void print_progress_bar()
	{
		auto & v_cl = create_vcluster();
Pietro Incardona's avatar
Pietro Incardona committed
258

259 260 261
		if (v_cl.getProcessUnitID() == 0)
			std::cout << "-----------------------25-----------------------50-----------------------75----------------------100" << std::endl;
	}
Pietro Incardona's avatar
Pietro Incardona committed
262

263 264 265 266 267 268 269 270 271 272 273
	/*! \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
274

275 276 277 278
		Mat A;
		Mat P;
		Vec x;
		Vec b;
Pietro Incardona's avatar
Pietro Incardona committed
279

280 281 282
		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
283

284 285 286
		solError err = statSolutionError(A,b,x,ksp);
		itError erri(err);
		erri.it = it;
Pietro Incardona's avatar
Pietro Incardona committed
287

288
        itError err_fill;
Pietro Incardona's avatar
Pietro Incardona committed
289

290 291 292 293 294 295 296 297 298 299 300 301 302
        size_t old_size = pts->bench.last().res.size();
        pts->bench.last().res.resize(it+1);

        if (old_size > 0)
                err_fill = pts->bench.last().res.get(old_size-1);
        else
                err_fill = erri;

        for (long int i = old_size ; i < (long int)it ; i++)
                pts->bench.last().res.get(i) = err_fill;

        // Add the error per iteration
        pts->bench.last().res.get(it) = erri;
Pietro Incardona's avatar
Pietro Incardona committed
303

Pietro Incardona's avatar
Pietro Incardona committed
304

305
		pts->progress(it);
Pietro Incardona's avatar
Pietro Incardona committed
306

307
		return 0;
Pietro Incardona's avatar
Pietro Incardona committed
308 309
	}

310 311 312 313 314 315 316 317 318
	/*! \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
319
	{
320
		petsc_solver<double> * pts = (petsc_solver *)data;
Pietro Incardona's avatar
Pietro Incardona committed
321

322
		pts->progress(it);
Pietro Incardona's avatar
Pietro Incardona committed
323

324 325
		return 0;
	}
Pietro Incardona's avatar
Pietro Incardona committed
326

327 328 329 330 331 332 333
	/*! \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
334

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

337
		auto & v_cl = create_vcluster();
Pietro Incardona's avatar
Pietro Incardona committed
338

339 340
		if (std::floor(it * 100.0 / n_max_it) > tmp)
		{
Pietro Incardona's avatar
Pietro Incardona committed
341
			size_t diff = it * 100.0 / n_max_it - tmp;
342 343
			tmp = it * 100.0 / n_max_it;
			if (v_cl.getProcessUnitID() == 0)
Pietro Incardona's avatar
Pietro Incardona committed
344 345 346 347
			{
				for (size_t k = 0 ; k <  diff ; k++)	{std::cout << "*";}
				std::cout << std::flush;
			}
348 349
		}
	}
Pietro Incardona's avatar
Pietro Incardona committed
350

351 352 353 354 355 356 357 358
	/*! \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
359

360 361
		if (v_cl.getProcessUnitID() == 0)
		{
Pietro Incardona's avatar
Pietro Incardona committed
362 363 364 365
			KSPType typ;
			KSPGetType(ksp,&typ);

			std::cout << "Method: " << typ << " " << " pre-conditoner: " << PCJACOBI << "  iterations: " << err.its << std::endl;
366 367 368
			std::cout << "Norm of error: " << err.err_norm << "   Norm infinity: " << err.err_inf << std::endl;
		}
	}
Pietro Incardona's avatar
Pietro Incardona committed
369

370 371 372 373 374 375 376 377 378 379 380 381 382 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
	/*! \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
416

417 418
		return std::string("UNKNOWN");
	}
Pietro Incardona's avatar
Pietro Incardona committed
419

420 421 422 423 424 425 426 427 428 429 430 431
	/*! \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
432

433 434 435 436 437 438 439 440 441 442 443 444 445 446 447
	/*! \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
448 449
	}

450 451 452 453 454 455 456 457
	/*! \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
458 459
	 *
	 */
460
	void try_solve_simple(Mat & A_, const Vec & b_, Vec & x_)
Pietro Incardona's avatar
Pietro Incardona committed
461
	{
462 463
		Vec best_sol;
		PETSC_SAFE_CALL(VecDuplicate(x_,&best_sol));
Pietro Incardona's avatar
Pietro Incardona committed
464

465 466
		// Best residual
		double best_res = std::numeric_limits<double>::max();
Pietro Incardona's avatar
Pietro Incardona committed
467

468 469
		// Create a new VCluster
		auto & v_cl = create_vcluster();
Pietro Incardona's avatar
Pietro Incardona committed
470

471
		destroyKSP();
Pietro Incardona's avatar
Pietro Incardona committed
472

473 474 475 476 477
		// 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
478

479 480 481 482 483 484 485
			// 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
486

487 488 489 490 491
					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
492

493 494 495 496 497 498 499 500 501 502 503 504
					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
505

506 507 508 509 510
					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
511

512 513 514 515 516 517 518
					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
519

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

522
				copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
Pietro Incardona's avatar
Pietro Incardona committed
523
			}
524 525 526 527 528 529 530 531

			destroyKSP();
		}

		write_bench_report();

		// Copy the best solution to x
		PETSC_SAFE_CALL(VecCopy(best_sol,x_));
Pietro Incardona's avatar
Pietro Incardona committed
532
		PETSC_SAFE_CALL(VecDestroy(&best_sol));
533 534 535 536 537 538 539 540 541 542 543 544 545 546 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
	}

	/*! \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
578 579 580 581 582 583 584
	}

	/*! \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
585
//		PETSC_SAFE_CALL(KSPSetType(ksp,s_type));
Pietro Incardona's avatar
Pietro Incardona committed
586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612

		// 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));
		}

613
		PETSC_SAFE_CALL(KSPSetFromOptions(ksp));
Pietro Incardona's avatar
Pietro Incardona committed
614 615 616 617 618 619 620 621 622 623 624 625

		// 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
	 *
	 */
626
	static solError statSolutionError(Mat & A_, const Vec & b_, Vec & x_, KSP ksp)
Pietro Incardona's avatar
Pietro Incardona committed
627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657
	{
		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;
658
		err.err_norm = norm / col;
Pietro Incardona's avatar
Pietro Incardona committed
659 660 661 662 663 664
		err.err_inf = norm_inf;
		err.its = its;

		return err;
	}

665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680
	/*! \brief initialize the KSP object
	 *
	 *
	 */
	void initKSP()
	{
		PETSC_SAFE_CALL(KSPCreate(PETSC_COMM_WORLD,&ksp));
	}

	/*! \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
681

682
		setMaxIter(maxits);
683 684 685 686 687 688 689 690 691 692 693 694 695

		// 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
696 697 698

public:

incardon's avatar
incardon committed
699 700
	typedef Vector<double,PETSC_BASE> return_type;

Pietro Incardona's avatar
Pietro Incardona committed
701 702
	~petsc_solver()
	{
incardon's avatar
incardon committed
703
		PETSC_SAFE_CALL(KSPDestroy(&ksp));
Pietro Incardona's avatar
Pietro Incardona committed
704 705
	}

Pietro Incardona's avatar
Pietro Incardona committed
706
	petsc_solver()
707
	:maxits(300)
Pietro Incardona's avatar
Pietro Incardona committed
708
	{
709
		initKSP();
Pietro Incardona's avatar
Pietro Incardona committed
710 711 712

		// Add the solvers

713 714
		solvs.add(std::string(KSPBCGS));
//		solvs.add(std::string(KSPIBCGS)); <--- Produce invalid argument
715
//		solvs.add(std::string(KSPFBCGS));
716 717
//		solvs.add(std::string(KSPFBCGSR)); <--- Nan production problems
		solvs.add(std::string(KSPBCGSL));
Pietro Incardona's avatar
Pietro Incardona committed
718
		solvs.add(std::string(KSPGMRES));
719 720 721
		solvs.add(std::string(KSPFGMRES));
		solvs.add(std::string(KSPLGMRES));
//		solvs.add(std::string(KSPPGMRES)); <--- Take forever
Pietro Incardona's avatar
Pietro Incardona committed
722
//		solvs.add(std::string(KSPGCR));
Pietro Incardona's avatar
Pietro Incardona committed
723 724

		setSolver(KSPGMRES);
Pietro Incardona's avatar
Pietro Incardona committed
725 726
	}

727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754
	/*! \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
755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774
	/*! \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
775
		PetscOptionsSetValue("-ksp_type",type);
Pietro Incardona's avatar
Pietro Incardona committed
776 777 778 779 780 781 782 783 784 785 786
	}

	/*! \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
787
		PetscOptionsSetValue("-ksp_rtol",std::to_string(rtol_).c_str());
Pietro Incardona's avatar
Pietro Incardona committed
788 789 790 791 792 793 794 795 796 797 798
	}

	/*! \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
799
		PetscOptionsSetValue("-ksp_atol",std::to_string(abstol_).c_str());
Pietro Incardona's avatar
Pietro Incardona committed
800 801 802 803 804 805 806 807 808 809 810
	}

	/*! \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
811 812 813 814 815 816 817 818 819 820
		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());
821
		maxits = n;
Pietro Incardona's avatar
Pietro Incardona committed
822 823 824 825 826 827 828 829
	}

	/*! 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
	 *
830
	 * \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
831 832 833 834
	 *
	 */
	void searchDirections(PetscInt l)
	{
835 836 837 838 839 840 841 842 843 844 845
		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
846 847 848 849 850 851 852 853 854 855
	}

	/*! \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
	 *
incardon's avatar
incardon committed
856 857 858 859 860 861
	 * \param A sparse matrix
	 * \param b vector
	 * \param initial_guess true if x has the initial guess
	 *
	 * \return the solution
	 *
Pietro Incardona's avatar
Pietro Incardona committed
862
	 */
incardon's avatar
incardon committed
863
	Vector<double,PETSC_BASE> solve(SparseMatrix<double,int,PETSC_BASE> & A, const Vector<double,PETSC_BASE> & b, bool initial_guess = false)
Pietro Incardona's avatar
Pietro Incardona committed
864 865 866 867 868 869 870 871 872 873 874 875
	{
		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));
incardon's avatar
incardon committed
876
		PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_FALSE));
Pietro Incardona's avatar
Pietro Incardona committed
877 878 879 880 881 882 883 884 885 886 887 888

		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;
	}
incardon's avatar
incardon committed
889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931

	/*! \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
	 *
	 * \param A sparse matrix
	 * \param b vector
	 * \param x solution and initial guess
	 * \param initial_guess true if x has the initial guess
	 *
	 * \return true if succeed
	 *
	 */
	bool solve(SparseMatrix<double,int,PETSC_BASE> & A, Vector<double,PETSC_BASE> & x, 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));
		PETSC_SAFE_CALL(KSPSetInitialGuessNonzero(ksp,PETSC_TRUE));


		Vec & x_ = x.getVec();

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

		x.update();
		return true;
	}
Pietro Incardona's avatar
Pietro Incardona committed
932 933
};

934
#endif
Pietro Incardona's avatar
Pietro Incardona committed
935 936

#endif /* OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_HPP_ */