Commit 63dda575 authored by Jakub Klinkovský's avatar Jakub Klinkovský
Browse files

Refactoring ILU preconditioners: added TriangularSolve.h

parent dc4b2498

src/TNL/Solvers/Linear/Preconditioners/CMakeLists.txt

+1 −0
Original line number Diff line number Diff line
@@ -5,6 +5,7 @@ SET( headers Preconditioner.h 
             ILU0_impl.h

             ILUT.h

             ILUT_impl.h

             TriangularSolve.h

   )



INSTALL( FILES ${headers} DESTINATION ${TNL_TARGET_INCLUDE_DIRECTORY}/Solvers/Linear/Preconditioners )

src/TNL/Solvers/Linear/Preconditioners/ILU0.h

+2 −2
Original line number Diff line number Diff line
@@ -62,8 +62,8 @@ public: 
   virtual void solve( ConstVectorViewType b, VectorViewType x ) const override;



protected:

   Matrices::CSR< RealType, DeviceType, IndexType > L;

   Matrices::CSR< RealType, DeviceType, IndexType > U;

   // The factors L and U are stored separately and the rows of U are reversed.

   Matrices::CSR< RealType, DeviceType, IndexType > L, U;

};



template< typename Matrix >

src/TNL/Solvers/Linear/Preconditioners/ILU0_impl.h

+16 −50
Original line number Diff line number Diff line
@@ -13,6 +13,7 @@ 
#pragma once



#include "ILU0.h"

#include "TriangularSolve.h"



#include <TNL/Exceptions/CudaSupportMissing.h>
@@ -107,46 +108,11 @@ void 
ILU0_impl< Matrix, Real, Devices::Host, Index >::

solve( ConstVectorViewType b, VectorViewType x ) const

{

   TNL_ASSERT_EQ( b.getSize(), L.getRows(), "wrong size of the right hand side" );

   TNL_ASSERT_EQ( x.getSize(), L.getRows(), "wrong size of the solution vector" );



   const IndexType N = x.getSize();



   // Step 1: solve y from Ly = b

   for( IndexType i = 0; i < N; i++ ) {

      x[ i ] = b[ i ];



      const auto L_entries = L.getRowLength( i );



      // this condition is to avoid segfaults on empty L.getRow( i )

      if( L_entries > 0 ) {

         const auto L_i = L.getRow( i );



         // loop for j = 0, ..., i - 1; but only over the non-zero entries

         for( IndexType c_j = 0; c_j < L_entries; c_j++ ) {

            const auto j = L_i.getElementColumn( c_j );

            x[ i ] -= L_i.getElementValue( c_j ) * x[ j ];

         }

      }

   }

   triangularSolveLower< true >( L, x, b );



   // Step 2: solve x from Ux = y

   for( IndexType i = N - 1; i >= 0; i-- ) {

      const IndexType U_idx = N - 1 - i;



      const auto U_entries = U.getRowLength( U_idx );

      const auto U_i = U.getRow( U_idx );



      const auto U_ii = U_i.getElementValue( 0 );



      // loop for j = i+1, ..., N-1; but only over the non-zero entries

      for( IndexType c_j = 1; c_j < U_entries ; c_j++ ) {

         const auto j = U_i.getElementColumn( c_j );

         x[ i ] -= U_i.getElementValue( c_j ) * x[ j ];

      }



      x[ i ] /= U_ii;

   }

   triangularSolveUpper< true, true >( U, x, x );

}

src/TNL/Solvers/Linear/Preconditioners/ILUT.h

+1 −2
Original line number Diff line number Diff line
@@ -71,8 +71,7 @@ protected: 
   Real tau = 1e-4;



   // The factors L and U are stored separately and the rows of U are reversed.

   Matrices::CSR< RealType, DeviceType, IndexType > L;

   Matrices::CSR< RealType, DeviceType, IndexType > U;

   Matrices::CSR< RealType, DeviceType, IndexType > L, U;

};



template< typename Matrix, typename Real, typename Index >

src/TNL/Solvers/Linear/Preconditioners/ILUT_impl.h

+4 −41
Original line number Diff line number Diff line
@@ -15,6 +15,8 @@ 
#include <vector>



#include "ILUT.h"

#include "TriangularSolve.h"



#include <TNL/Timer.h>



namespace TNL {
@@ -245,50 +247,11 @@ void 
ILUT_impl< Matrix, Real, Devices::Host, Index >::

solve( ConstVectorViewType b, VectorViewType x ) const

{

   TNL_ASSERT_EQ( b.getSize(), L.getRows(), "wrong size of the right hand side" );

   TNL_ASSERT_EQ( x.getSize(), L.getRows(), "wrong size of the solution vector" );



   const IndexType N = x.getSize();



   // Step 1: solve y from Ly = b

   for( IndexType i = 0; i < N; i++ ) {

      x[ i ] = b[ i ];



      const auto L_entries = L.getRowLength( i );



      // this condition is to avoid segfaults on empty L.getRow( i )

      if( L_entries > 0 ) {

         const auto L_i = L.getRow( i );



         // loop for j = 0, ..., i - 1; but only over the non-zero entries

         for( IndexType c_j = 0; c_j < L_entries; c_j++ ) {

            const auto j = L_i.getElementColumn( c_j );

            // we might have allocated more space than was actually needed due to dropping

            if( j >= N ) break;

            x[ i ] -= L_i.getElementValue( c_j ) * x[ j ];

         }

      }

   }

   triangularSolveLower< false >( L, x, b );



   // Step 2: solve x from Ux = y

   for( IndexType i = N - 1; i >= 0; i-- ) {

      const IndexType U_idx = N - 1 - i;



      const auto U_entries = U.getRowLength( U_idx );

      const auto U_i = U.getRow( U_idx );



      const auto U_ii = U_i.getElementValue( 0 );



      // loop for j = i+1, ..., N-1; but only over the non-zero entries

      for( IndexType c_j = 1; c_j < U_entries ; c_j++ ) {

         const auto j = U_i.getElementColumn( c_j );

         // we might have allocated more space than was actually needed due to dropping

         if( j >= N ) break;

         x[ i ] -= U_i.getElementValue( c_j ) * x[ j ];

      }



      x[ i ] /= U_ii;

   }

   triangularSolveUpper< true, false >( U, x, x );

}



} // namespace Preconditioners