First version of ILU0 subdomains solver.

Reorder header list and made it more diff friendly.

[[Imported from SVN: r1251]]

dune/istl/Makefile.am

@@ -3,13 +3,45 @@
 SUBDIRS = . tutorial test paamg
 
 istldir = $(includedir)/dune/istl
-istl_HEADERS = allocator.hh basearray.hh bcrsmatrix.hh bvector.hh gsetc.hh io.hh istlexception.hh vbvector.hh \
-	ilu.hh operators.hh preconditioners.hh solvers.hh indexset.hh communicator.hh remoteindices.hh \
-	interface.hh indicessyncer.hh matrixindexset.hh selection.hh plocalindex.hh localindex.hh scalarproducts.hh \
-	solvercategory.hh bdmatrix.hh matrixutils.hh schwarz.hh owneroverlapcopy.hh matrix.hh supermatrix.hh superlu.hh \
-	pardiso.hh overlappingschwarz.hh scaledidmatrix.hh diagonalmatrix.hh btdmatrix.hh \
-	matrixmatrix.hh repartition.hh matrixredistribute.hh \
-	matrixmarket.hh
+istl_HEADERS = allocator.hh \
+	basearray.hh \
+	bcrsmatrix.hh \
+	bdmatrix.hh \
+	btdmatrix.hh \
+	bvector.hh \
+	communicator.hh \
+	diagonalmatrix.hh \
+	gsetc.hh \
+	ilu.hh \
+	indexset.hh \
+	indicessyncer.hh \
+	interface.hh \
+	io.hh \
+	istlexception.hh \
+	localindex.hh \
+	matrix.hh \
+	matrixindexset.hh \
+	matrixmarket.hh\
+	matrixmatrix.hh \
+	matrixredistribute.hh \
+	matrixutils.hh \
+	operators.hh \
+	overlappingschwarz.hh \
+	owneroverlapcopy.hh \
+	pardiso.hh \
+	plocalindex.hh \
+	preconditioners.hh \
+	remoteindices.hh \
+	repartition.hh \
+	selection.hh \
+	scalarproducts.hh \
+	scaledidmatrix.hh \
+	schwarz.hh \
+	solvercategory.hh \
+	solvers.hh \
+	supermatrix.hh \
+	superlu.hh \
+	vbvector.hh
 
 
 include $(top_srcdir)/am/global-rules

dune/istl/ilusubdomainsolver.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_ISTL_ILUSUBDOMAIN_HH
+#define DUNE_ISTL_ILUSUBDOMAIN_HH
+
+#include <map>
+#include <dune/common/typetraits.hh>
+
+namespace Dune {
+
+  template<class M, class X, class Y>
+  class ILU0SubdomainSolver  {
+  public:
+    //! \brief The matrix type the preconditioner is for.
+    typedef typename Dune::remove_const<M>::type matrix_type;
+    //! \brief The domain type of the preconditioner.
+    typedef X domain_type;
+    //! \brief The range type of the preconditioner.
+    typedef Y range_type;
+    typedef typename X::field_type field_type;
+
+    void apply (X& v, const Y& d)
+    {
+      bilu_backsolve(ILU,v,d);
+    }
+
+    template<class S>
+    void setMatrix(const M& A, S& rowset);
+
+  private:
+    //! \brief The relaxation factor to use.
+    field_type _w;
+    //! \brief The ILU0 decomposition of the matrix.
+    matrix_type ILU;
+  };
+
+  template<class M, class X, class Y>
+  template<class S>
+  void ILU0SubdomainSolver<M,X,Y>::setMatrix(const M& A, S& rowSet)
+  {
+    // Calculate consecutive indices for local problem
+    // while perserving the ordering
+    typedef typename M::size_type size_type;
+    typedef std::map<typename S::key_type,size_type> IndexMap;
+    typedef typename IndexMap::iterator IMIter;
+    IndexMap indexMap;
+    IMIter guess = indexMap.begin();
+    size_type localIndex=0;
+
+    typedef typename S::const_iterator SIter;
+    for(SIter rowIdx = rowSet.begin(), rowEnd=rowSet.end();
+        rowIdx!= rowEnd; ++rowIdx, ++localIndex)
+      guess = indexMap.insert(std::make_pair(*rowIdx,localIndex),
+                              guess);
+
+    // Build Matrix for local subproblem
+    ILU.setSize(rowSet.size(),rowSet.size());
+    ILU.setBuildMode(matrix_type::row_wise);
+
+    // Create sparsity pattern
+    typedef typename matrix_type::CreateIterator CIter;
+    CIter rowCreator = ILU.createbegin();
+    typedef typename S::const_iterator RIter;
+    for(SIter rowIdx = rowSet.begin(), rowEnd=rowSet.end();
+        rowIdx!= rowEnd; ++rowIdx, ++rowCreator) {
+      // See wich row entries are in our subset and add them to
+      // the sparsity pattern
+      guess = indexMap.begin();
+      for(typename matrix_type::ConstColIterator& col=A[rowIdx].begin(),
+          endcol=A[rowIdx].end(); col != endcol; ++col) {
+        // search for the entry in the row set
+        IMIter oldguess = guess;
+        guess = indexMap.find(col.index(), guess);
+        if(guess!=indexMap.end())
+          // add local index to row
+          rowCreator.insert(guess->second);
+        else
+          guess = oldguess;
+      }
+    }
+
+    // Insert the matrix values for the local problem
+    typename matrix_type::iterator iluRow=ILU.begin();
+
+    for(SIter rowIdx = rowSet.begin(), rowEnd=rowSet.end();
+        rowIdx!= rowEnd; ++rowIdx, ++iluRow) {
+      // See wich row entries are in our subset and add them to
+      // the sparsity pattern
+      guess = indexMap.begin();
+      for(typename matrix_type::ConstColIterator& col=A[rowIdx].begin(),
+          endcol=A[rowIdx].end(), localCol=iluRow.begin(); col != endcol; ++col) {
+        // search for the entry in the row set
+        IMIter oldguess = guess;
+        guess = indexMap.find(col.index(), guess);
+        if(guess!=indexMap.end()) {
+          // set local value
+          *localCol*col;
+          ++localCol;
+        }
+      }
+    }
+    bilu0_decomposition(ILU);
+  }
+
+} // end name space DUNE
+
+
+#endif