* added some docu

* made the local solver exchangeable

[[Imported from SVN: r7058]]

dune/solvers/iterationsteps/truncatedblockgsstep.hh

@@ -117,6 +117,14 @@ struct RecursiveGSStep<1>
 
 } // end namespace TruncatedBlockGSStepNS
 
+/**  \brief nonrecursive specialization of the TruncatedBlockGSStep
+  *
+  *  \tparam SparseMatrixType the type of the sparse matrix in the SumOperator
+  *  \tparam LowRankMatrixType the type of the lowrank matrix in the SumOperator
+  *  \tparam VectorType the type of the solution vector
+  *  \tparam BitVectorType the type of the truncation field
+  *  \tparam localSolver static switch for local solver, see specializations of method localSolve
+  */
 template<class SparseMatrixType, class LowRankMatrixType, class VectorType, class BitVectorType>
 class TruncatedBlockGSStep<SumOperator<SparseMatrixType, LowRankMatrixType>, VectorType, BitVectorType>
     : public LinearIterationStep<SumOperator<SparseMatrixType, LowRankMatrixType>, VectorType, BitVectorType>
@@ -125,14 +133,22 @@ class TruncatedBlockGSStep<SumOperator<SparseMatrixType, LowRankMatrixType>, Vec
     typedef typename VectorType::block_type VectorBlock;
     typedef typename VectorType::field_type field_type;
 
+    typedef typename SparseMatrixType::block_type SpBlock;
+    typedef typename LowRankMatrixType::block_type LRBlock;
+    typedef typename StaticMatrix::Promote<SpBlock, LRBlock>::Type MBlock;
+    typedef typename VectorType::block_type VBlock;
+
 public:
 
     //! Default constructor.  Doesn't init anything
-    TruncatedBlockGSStep() {}
+    TruncatedBlockGSStep(int localSolver=0):
+        localSolver(localSolver)
+    {}
 
     //! Constructor with a linear problem
-    TruncatedBlockGSStep(MatrixType& mat, VectorType& x, VectorType& rhs)
-        : LinearIterationStep<MatrixType,VectorType>(mat, x, rhs)
+    TruncatedBlockGSStep(MatrixType& mat, VectorType& x, VectorType& rhs, int localSolver=0):
+        LinearIterationStep<MatrixType,VectorType>(mat, x, rhs),
+        localSolver(localSolver)
     {}
 
     virtual VectorType getSol()
@@ -145,10 +161,6 @@ public:
     {
         typedef typename SparseMatrixType::row_type::ConstIterator SparseMatrixColumnIterator;
         typedef typename LowRankMatrixType::LowRankFactorType::row_type::ConstIterator LowRankFactorColumnIterator;
-        typedef typename SparseMatrixType::block_type SpBlock;
-        typedef typename LowRankMatrixType::block_type LRBlock;
-        typedef typename StaticMatrix::Promote<SpBlock, LRBlock>::Type MBlock;
-        typedef typename VectorType::block_type VBlock;
 
         const MatrixType& mat = *this->mat_;
         const SparseMatrixType& sparse_mat = mat.sparseMatrix();
@@ -200,11 +212,12 @@ public:
 
             // "solve" Local system by one GS step
 //            if (Aii!=0)
+#if 0
             {
                 for(typename MBlock::size_type i=0; i<Aii.N(); ++i)
                 {
                     const typename MBlock::field_type& aii = Aii[i][i];
-                    if (aii>0)
+                    if (fabs(aii)>1e-13)
                     {
                         if (not(ignore[row][i]))
                         {
@@ -219,6 +232,12 @@ public:
                     }
                 }
             }
+#endif
+            if (localSolver==1)
+                localExactSolve(Aii, r, x[row], ignore[row]);
+            else
+                localGSSolve(Aii, r, x[row], ignore[row]);
+
 
             // update the mx
             for (typename LowRankMatrixType::size_type k=0; k<lowrank_mat.lowRankFactor().N(); ++k)
@@ -228,6 +247,63 @@ public:
             }
         }
     }
+
+    private:
+
+        int localSolver;
+
+        /**  \brief solve method for local problems
+          *
+          *  This version "solves" by one GS step.
+          *  This is the "default" behaviour
+          */
+        void localGSSolve(MBlock& A, VBlock& b, VBlock& x, typename BitVectorType::const_reference ignore)
+        {
+            for(typename MBlock::size_type i=0; i<A.N(); ++i)
+            {
+                const typename MBlock::field_type& aii = A[i][i];
+                if (fabs(aii)>1e-13)
+                {
+                    if (not(ignore[i]))
+                    {
+                        x[i] = b[i];
+                        typename MBlock::row_type::Iterator inner_it = A[i].begin();
+                        typename MBlock::row_type::Iterator inner_end = A[i].end();
+                        for(; inner_it!=inner_end; ++inner_it)
+                            if (inner_it.index()!=i)
+                                x[i] -= (*inner_it) * x[inner_it.index()];
+                        x[i] /= aii;
+                    }
+                }
+            }
+        }
+
+        /**  \brief solve method for local problems
+          *
+          *  uses the solve method of FieldMatrix
+          */
+        void localExactSolve(MBlock& A, VBlock& b, VBlock& x, typename BitVectorType::const_reference ignore)
+        {
+            for(typename MBlock::size_type i=0; i<A.N(); ++i)
+            {
+                if (ignore[i])
+                {
+                    A[i] = 0.0;
+                }
+
+                typename MBlock::row_type::Iterator inner_it = A[i].begin();
+                typename MBlock::row_type::Iterator inner_end = A[i].end();
+                for(; inner_it!=inner_end; ++inner_it)
+                    if (inner_it.index()==i and *inner_it==0.0)
+                    {
+                        *inner_it = 1.0;
+                        b[i] = x[i];
+                    }
+            }
+
+            A.solve(x,b);
+
+        }
 };
 
 #endif