diff --git a/dune/solvers/iterationsteps/truncatedblockgsstep.hh b/dune/solvers/iterationsteps/truncatedblockgsstep.hh index 4d1c0c522732cf6f72c0d8ffbbb00b858a4de394..67b9accffd9c134d8906cb0b9869c601d146789f 100644 --- a/dune/solvers/iterationsteps/truncatedblockgsstep.hh +++ b/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