From 4f66a8f2e3dba3e431e3efadec356ada77fc680d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jakub=20Klinkovsk=C3=BD?= <klinkjak@fjfi.cvut.cz>
Date: Mon, 22 Feb 2016 21:33:44 +0100
Subject: [PATCH] Adding preconditioner support into the TFQMR solver.
---
src/solvers/linear/krylov/tnlTFQMRSolver.h | 4 +-
.../linear/krylov/tnlTFQMRSolver_impl.h | 64 +++++++++++++------
2 files changed, 49 insertions(+), 19 deletions(-)
diff --git a/src/solvers/linear/krylov/tnlTFQMRSolver.h b/src/solvers/linear/krylov/tnlTFQMRSolver.h
index ab25f3e349..32d250f52a 100644
--- a/src/solvers/linear/krylov/tnlTFQMRSolver.h
+++ b/src/solvers/linear/krylov/tnlTFQMRSolver.h
@@ -75,7 +75,9 @@ class tnlTFQMRSolver : public tnlObject,
bool setSize( IndexType size );
- tnlVector< RealType, Device, IndexType > d, r, w, u, v, r_ast, Au;
+ tnlVector< RealType, Device, IndexType > d, r, w, u, v, r_ast, Au, M_tmp;
+
+ IndexType size;
const MatrixType* matrix;
const PreconditionerType* preconditioner;
diff --git a/src/solvers/linear/krylov/tnlTFQMRSolver_impl.h b/src/solvers/linear/krylov/tnlTFQMRSolver_impl.h
index b5e72e6456..5ddb54a5ec 100644
--- a/src/solvers/linear/krylov/tnlTFQMRSolver_impl.h
+++ b/src/solvers/linear/krylov/tnlTFQMRSolver_impl.h
@@ -18,17 +18,12 @@
#ifndef tnlTFQMRSolver_implH
#define tnlTFQMRSolver_implH
-template< typename RealType,
- typename Vector >
-RealType computeTFQMRNewP( Vector& p,
- const Vector&r,
- const RealType& beta,
- const RealType& omega,
- const Vector& Ap );
-
template< typename Matrix,
typename Preconditioner >
tnlTFQMRSolver< Matrix, Preconditioner > :: tnlTFQMRSolver()
+: size( 0 ),
+ matrix( 0 ),
+ preconditioner( 0 )
{
}
@@ -83,15 +78,29 @@ bool tnlTFQMRSolver< Matrix, Preconditioner > :: solve( const Vector& b, Vector&
dbgFunctionName( "tnlTFQMRSolver", "Solve" );
if( ! this -> setSize( matrix -> getRows() ) ) return false;
- RealType tau, theta, eta, rho, alpha, w_norm;
- RealType b_norm = b. lpNorm( 2.0 );
- if( b_norm == 0.0 )
- b_norm = 1.0;
+ RealType tau, theta, eta, rho, alpha, b_norm, w_norm;
- this -> matrix -> vectorProduct( x, r );
- r. addVector( b, 1.0, -1.0 );
+ if( preconditioner ) {
+ preconditioner -> solve( b, M_tmp );
+ b_norm = M_tmp. lpNorm( ( RealType ) 2.0 );
+
+ matrix -> vectorProduct( x, M_tmp );
+ M_tmp.addVector( b, 1.0, -1.0 );
+ preconditioner -> solve( M_tmp, r );
+ }
+ else {
+ b_norm = b. lpNorm( 2.0 );
+ matrix -> vectorProduct( x, r );
+ r.addVector( b, 1.0, -1.0 );
+ }
w = u = r;
- matrix -> vectorProduct( u, Au );
+ if( preconditioner ) {
+ matrix -> vectorProduct( u, M_tmp );
+ preconditioner -> solve( M_tmp, Au );
+ }
+ else {
+ matrix -> vectorProduct( u, Au );
+ }
v = Au;
d. setValue( 0.0 );
tau = r. lpNorm( 2.0 );
@@ -101,6 +110,9 @@ bool tnlTFQMRSolver< Matrix, Preconditioner > :: solve( const Vector& b, Vector&
// only to avoid compiler warning; alpha is initialized inside the loop
alpha = 0.0;
+ if( b_norm == 0.0 )
+ b_norm = 1.0;
+
this->resetIterations();
this->setResidue( tau / b_norm );
@@ -114,7 +126,13 @@ bool tnlTFQMRSolver< Matrix, Preconditioner > :: solve( const Vector& b, Vector&
else {
// not necessary in odd iter since the previous iteration
// already computed v_{m+1} = A*u_{m+1}
- matrix -> vectorProduct( u, Au );
+ if( preconditioner ) {
+ matrix -> vectorProduct( u, M_tmp );
+ preconditioner -> solve( M_tmp, Au );
+ }
+ else {
+ matrix -> vectorProduct( u, Au );
+ }
}
w.addVector( Au, -alpha );
d.addVector( u, 1.0, theta * theta * eta / alpha );
@@ -137,7 +155,13 @@ bool tnlTFQMRSolver< Matrix, Preconditioner > :: solve( const Vector& b, Vector&
u.addVector( w, 1.0, beta );
v.addVector( Au, beta, beta * beta );
- matrix -> vectorProduct( u, Au );
+ if( preconditioner ) {
+ matrix -> vectorProduct( u, M_tmp );
+ preconditioner -> solve( M_tmp, Au );
+ }
+ else {
+ matrix -> vectorProduct( u, Au );
+ }
v.addVector( Au, 1.0 );
}
else {
@@ -165,13 +189,17 @@ template< typename Matrix,
typename Preconditioner >
bool tnlTFQMRSolver< Matrix, Preconditioner > :: setSize( IndexType size )
{
+ if( this->size == size )
+ return true;
+ this->size = size;
if( ! d. setSize( size ) ||
! r. setSize( size ) ||
! w. setSize( size ) ||
! u. setSize( size ) ||
! v. setSize( size ) ||
! r_ast. setSize( size ) ||
- ! Au. setSize( size ) )
+ ! Au. setSize( size ) ||
+ ! M_tmp. setSize( size ) )
{
cerr << "I am not able to allocate all supporting vectors for the TFQMR solver." << endl;
return false;
--
GitLab