"Fossies" - the Fresh Open Source Software Archive

Member "dune-istl-2.7.1/dune/istl/ildl.hh" (26 Nov 2020, 6543 Bytes) of package /linux/misc/dune/dune-istl-2.7.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "ildl.hh" see the Fossies "Dox" file reference documentation and the last Fossies "Diffs" side-by-side code changes report: 2.6.0_vs_2.7.0.

    1 #ifndef DUNE_ISTL_ILDL_HH
    2 #define DUNE_ISTL_ILDL_HH
    3 
    4 #include <dune/common/scalarvectorview.hh>
    5 #include <dune/common/scalarmatrixview.hh>
    6 #include "ilu.hh"
    7 
    8 /**
    9  * \file
   10  *
   11  * \brief   Incomplete LDL decomposition
   12  * \author  Martin Nolte
   13  **/
   14 
   15 namespace Dune
   16 {
   17 
   18   // bildl_subtractBCT
   19   // -----------------
   20 
   21   template< class K, int m, int n >
   22   inline static void bildl_subtractBCT ( const FieldMatrix< K, m, n > &B, const FieldMatrix< K, m, n > &CT, FieldMatrix< K, m, n > &A )
   23   {
   24     for( int i = 0; i < m; ++i )
   25     {
   26       for( int j = 0; j < n; ++j )
   27       {
   28         for( int k = 0; k < n; ++k )
   29           A[ i ][ j ] -= B[ i ][ k ] * CT[ j ][ k ];
   30       }
   31     }
   32   }
   33 
   34   template< class K >
   35   inline static void bildl_subtractBCT ( const K &B, const K &CT, K &A,
   36                                          typename std::enable_if_t<Dune::IsNumber<K>::value>* sfinae = nullptr )
   37   {
   38     A -= B * CT;
   39   }
   40 
   41   template< class Matrix >
   42   inline static void bildl_subtractBCT ( const Matrix &B, const Matrix &CT, Matrix &A,
   43                                          typename std::enable_if_t<!Dune::IsNumber<Matrix>::value>* sfinae = nullptr )
   44   {
   45     for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
   46     {
   47       auto &&A_i = *i;
   48       auto &&B_i = B[ i.index() ];
   49       const auto ikend = B_i.end();
   50       for( auto j = A_i.begin(), jend = A_i.end(); j != jend; ++j )
   51       {
   52         auto &&A_ij = *j;
   53         auto &&CT_j = CT[ j.index() ];
   54         const auto jkend = CT_j.end();
   55         for( auto ik = B_i.begin(), jk = CT_j.begin(); (ik != ikend) && (jk != jkend); )
   56         {
   57           if( ik.index() == jk.index() )
   58           {
   59             bildl_subtractBCT( *ik, *jk, A_ij );
   60             ++ik; ++jk;
   61           }
   62           else if( ik.index() < jk.index() )
   63             ++ik;
   64           else
   65             ++jk;
   66         }
   67       }
   68     }
   69   }
   70 
   71 
   72 
   73   // bildl_decompose
   74   // ---------------
   75 
   76   /**
   77    * \brief  compute ILDL decomposition of a symmetric matrix A
   78    * \author Martin Nolte
   79    *
   80    * \param[inout]  A  matrix to decompose
   81    *
   82    * \note A is overwritten by the factorization.
   83    * \note Only the lower half of A is used.
   84    **/
   85   template< class Matrix >
   86   inline void bildl_decompose ( Matrix &A )
   87   {
   88     for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
   89     {
   90       auto &&A_i = *i;
   91 
   92       auto ij = A_i.begin();
   93       for( ; ij.index() < i.index(); ++ij )
   94       {
   95         auto &&A_ij = *ij;
   96         auto &&A_j = A[ ij.index() ];
   97 
   98         // store L_ij Dj in A_ij (note: for k < i: A_kj = L_kj)
   99         // L_ij Dj = A_ij - \sum_{k < j} (L_ik D_k) L_jk^T
  100         auto ik = A_i.begin();
  101         auto jk = A_j.begin();
  102         while( (ik != ij) && (jk.index() < ij.index()) )
  103         {
  104           if( ik.index() == jk.index() )
  105           {
  106             bildl_subtractBCT(*ik, *jk, A_ij);
  107             ++ik; ++jk;
  108           }
  109           else if( ik.index() < jk.index() )
  110             ++ik;
  111           else
  112             ++jk;
  113         }
  114       }
  115 
  116       if( ij.index() != i.index() )
  117         DUNE_THROW( ISTLError, "diagonal entry missing" );
  118 
  119       // update diagonal and multiply A_ij by D_j^{-1}
  120       auto &&A_ii = *ij;
  121       for( auto ik = A_i.begin(); ik != ij; ++ik )
  122       {
  123         auto &&A_ik = *ik;
  124         const auto &A_k = A[ ik.index() ];
  125 
  126         auto B = A_ik;
  127         Impl::asMatrix(A_ik).rightmultiply( Impl::asMatrix(*A_k.find( ik.index() )) );
  128         bildl_subtractBCT( B, A_ik, A_ii );
  129       }
  130       try
  131       {
  132         Impl::asMatrix(A_ii).invert();
  133       }
  134       catch( const Dune::FMatrixError &e )
  135       {
  136         DUNE_THROW( MatrixBlockError, "ILDL failed to invert matrix block A[" << i.index() << "][" << ij.index() << "]" << e.what(); th__ex.r = i.index(); th__ex.c = ij.index() );
  137       }
  138     }
  139   }
  140 
  141 
  142 
  143   // bildl_backsolve
  144   // ---------------
  145 
  146   template< class Matrix, class X, class Y >
  147   inline void bildl_backsolve ( const Matrix &A, X &v, const Y &d, bool isLowerTriangular = false )
  148   {
  149     // solve L v = d, note: Lii = I
  150     for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
  151     {
  152       const auto &A_i = *i;
  153       v[ i.index() ] = d[ i.index() ];
  154       for( auto ij = A_i.begin(); ij.index() < i.index(); ++ij )
  155       {
  156         auto&& vi = Impl::asVector( v[ i.index() ] );
  157         Impl::asMatrix(*ij).mmv(Impl::asVector( v[ ij.index() ] ), vi);
  158       }
  159     }
  160 
  161     // solve D w = v, note: diagonal stores Dii^{-1}
  162     if( isLowerTriangular )
  163     {
  164       // The matrix is lower triangular, so the diagonal entry is the
  165       // last one in each row.
  166       for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
  167       {
  168         const auto &A_i = *i;
  169         const auto ii = A_i.beforeEnd();
  170         assert( ii.index() == i.index() );
  171         // We need to be careful here: Directly using
  172         // auto rhs = Impl::asVector(v[ i.index() ]);
  173         // is not OK in case this is a proxy. Hence
  174         // we first have to copy the value. Notice that
  175         // this is still not OK, if the vector type itself returns
  176         // proxy references.
  177         auto rhsValue = v[ i.index() ];
  178         auto&& rhs = Impl::asVector(rhsValue);
  179         auto&& vi = Impl::asVector( v[ i.index() ] );
  180         Impl::asMatrix(*ii).mv(rhs, vi);
  181       }
  182     }
  183     else
  184     {
  185       // Without assumptions on the sparsity pattern we have to search
  186       // for the diagonal entry in each row.
  187       for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
  188       {
  189         const auto &A_i = *i;
  190         const auto ii = A_i.find( i.index() );
  191         assert( ii.index() == i.index() );
  192         // We need to be careful here: Directly using
  193         // auto rhs = Impl::asVector(v[ i.index() ]);
  194         // is not OK in case this is a proxy. Hence
  195         // we first have to copy the value. Notice that
  196         // this is still not OK, if the vector type itself returns
  197         // proxy references.
  198         auto rhsValue = v[ i.index() ];
  199         auto&& rhs = Impl::asVector(rhsValue);
  200         auto&& vi = Impl::asVector( v[ i.index() ] );
  201         Impl::asMatrix(*ii).mv(rhs, vi);
  202       }
  203     }
  204 
  205     // solve L^T v = w, note: only L is stored
  206     // note: we perform the operation column-wise from right to left
  207     for( auto i = A.beforeEnd(), iend = A.beforeBegin(); i != iend; --i )
  208     {
  209       const auto &A_i = *i;
  210       for( auto ij = A_i.begin(); ij.index() < i.index(); ++ij )
  211       {
  212         auto&& vij = Impl::asVector( v[ ij.index() ] );
  213         Impl::asMatrix(*ij).mmtv(Impl::asVector( v[ i.index() ] ), vij);
  214       }
  215     }
  216   }
  217 
  218 } // namespace Dune
  219 
  220 #endif // #ifndef DUNE_ISTL_ILDL_HH