// $Id$
// Version: $Name$
//
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2009 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
class PreconditionSSOR : public PreconditionRelaxation<MATRIX>
{
public:
+
+ /**
+ * A typedef to the base class.
+ */
+ typedef PreconditionRelaxation<MATRIX> BaseClass;
+
+
+ /**
+ * Initialize matrix and
+ * relaxation parameter. The
+ * matrix is just stored in the
+ * preconditioner object. The
+ * relaxation parameter should be
+ * larger than zero and smaller
+ * than 2 for numerical
+ * reasons. It defaults to 1.
+ */
+ void initialize (const MATRIX &A,
+ typename BaseClass::AdditionalData parameters = BaseClass::AdditionalData());
+
/**
* Apply preconditioner.
*/
*/
template<class VECTOR>
void Tvmult (VECTOR&, const VECTOR&) const;
+
+ private:
+ /**
+ * An array that stores for each matrix
+ * row where the first position after
+ * the diagonal is located.
+ */
+ std::vector<unsigned int> pos_right_of_diagonal;
};
//---------------------------------------------------------------------------
+template <class MATRIX>
+inline void
+PreconditionSSOR<MATRIX>::initialize (const MATRIX &rA,
+ typename BaseClass::AdditionalData parameters)
+{
+ this->PreconditionRelaxation<MATRIX>::initialize (rA, parameters);
+
+ // calculate the positions first after
+ // the diagonal.
+ const std::size_t * rowstart_ptr =
+ this->A->get_sparsity_pattern().get_rowstart_indices();
+ const unsigned int * const colnums =
+ this->A->get_sparsity_pattern().get_column_numbers();
+ const unsigned int n = this->A->n();
+ pos_right_of_diagonal.resize(n);
+ for (unsigned int row=0; row<n; ++row, ++rowstart_ptr)
+ {
+ // find the first element in this line
+ // which is on the right of the diagonal.
+ // we need to precondition with the
+ // elements on the left only.
+ // note: the first entry in each
+ // line denotes the diagonal element,
+ // which we need not check.
+ pos_right_of_diagonal[row] =
+ std::lower_bound (&colnums[*rowstart_ptr+1],
+ &colnums[*(rowstart_ptr+1)],
+ row)
+ - colnums;
+ }
+}
+
+
template <class MATRIX>
template<class VECTOR>
inline void
PreconditionSSOR<MATRIX>::vmult (VECTOR &dst, const VECTOR &src) const
{
Assert (this->A!=0, ExcNotInitialized());
- this->A->precondition_SSOR (dst, src, this->relaxation);
+ this->A->precondition_SSOR (dst, src, this->relaxation, pos_right_of_diagonal);
}
//---------------------------------------------------------------------------
-// Copyright (C) 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+// Copyright (C) 2002, 2003, 2004, 2005, 2006, 2009 by the deal.II authors
// by the deal.II authors and Stephen "Cheffo" Kolaroff
//
// This file is subject to QPL and may not be distributed
//---------------------------------------------------------------------------
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008, 2009 by the deal.II authors
// by the deal.II authors and Stephen "Cheffo" Kolaroff
//
// This file is subject to QPL and may not be distributed
* @<double@></tt>; others can be generated in application programs (see the
* section on @ref Instantiations in the manual).
*
- * @author Wolfgang Bangerth, 2008; unified interface: Ralf Hartmann
+ * @author Wolfgang Bangerth, 2008, 2009; unified interface: Ralf Hartmann
*/
template <typename number>
class SparseILU : public SparseLUDecomposition<number>
//---------------------------------------------------------------------------
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008, 2009 by the deal.II authors
// by the deal.II authors and Stephen "Cheffo" Kolaroff
//
// This file is subject to QPL and may not be distributed
// $Id$
// Version: $Name: $
//
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
/**
* Return the number of actually
- * nonzero elements of this
- * matrix.
+ * nonzero elements of this matrix. It
+ * is possible to specify the parameter
+ * <tt>threshold</tt> in order to count
+ * only the elements that have absolute
+ * value greater than the threshold.
*
- * Note, that this function does
- * (in contrary to
- * n_nonzero_elements()) not
- * count all entries of the
- * sparsity pattern but only the
- * ones that are nonzero.
- */
- unsigned int n_actually_nonzero_elements () const;
+ * Note, that this function does (in
+ * contrary to n_nonzero_elements())
+ * not count all entries of the
+ * sparsity pattern but only the ones
+ * that are nonzero (or whose absolute
+ * value is greater than threshold).
+ */
+ unsigned int n_actually_nonzero_elements (const double threshold = 0.) const;
/**
* Return a (constant) reference
* <tt>src</tt>.
*/
template <typename somenumber>
- void precondition_SSOR (Vector<somenumber> &dst,
- const Vector<somenumber> &src,
- const number om = 1.) const;
+ void precondition_SSOR (Vector<somenumber> &dst,
+ const Vector<somenumber> &src,
+ const number om = 1.,
+ const std::vector<unsigned int>&pos_right_of_diagonal=std::vector<unsigned int>()) const;
/**
* Apply SOR preconditioning
// $Id$
// Version: $Name$
//
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
template <typename number>
template <typename somenumber>
void
-SparseMatrix<number>::precondition_SSOR (Vector<somenumber> &dst,
- const Vector<somenumber> &src,
- const number om) const
+SparseMatrix<number>::precondition_SSOR (Vector<somenumber> &dst,
+ const Vector<somenumber> &src,
+ const number om,
+ const std::vector<unsigned int> &pos_right_of_diagonal) const
{
// to understand how this function works
// you may want to take a look at the CVS
const std::size_t *rowstart_ptr = &cols->rowstart[0];
somenumber *dst_ptr = &dst(0);
+ // case when we have stored the position
+ // just right of the diagonal (then we
+ // don't have to search for it).
+ if (pos_right_of_diagonal.size() != 0)
+ {
+ Assert (pos_right_of_diagonal.size() == dst.size(),
+ ExcDimensionMismatch (pos_right_of_diagonal.size(), dst.size()));
+
+ // forward sweep
+ for (unsigned int row=0; row<n; ++row, ++dst_ptr, ++rowstart_ptr)
+ {
+ *dst_ptr = src(row);
+ const unsigned int first_right_of_diagonal_index =
+ pos_right_of_diagonal[row];
+ Assert (first_right_of_diagonal_index <= *(rowstart_ptr+1),
+ ExcInternalError());
+ double s = *dst_ptr;
+ for (unsigned int j=(*rowstart_ptr)+1; j<first_right_of_diagonal_index; ++j)
+ s -= val[j] * dst(cols->colnums[j]);
+
+ // divide by diagonal element
+ *dst_ptr = s * om / val[*rowstart_ptr];
+ };
+
+ rowstart_ptr = &cols->rowstart[0];
+ dst_ptr = &dst(0);
+ for (unsigned int row=0; row<n; ++row, ++rowstart_ptr, ++dst_ptr)
+ *dst_ptr *= (2.-om)*val[*rowstart_ptr];
+
+ // backward sweep
+ rowstart_ptr = &cols->rowstart[n-1];
+ dst_ptr = &dst(n-1);
+ for (int row=n-1; row>=0; --row, --rowstart_ptr, --dst_ptr)
+ {
+ const unsigned int end_row = *(rowstart_ptr+1);
+ const unsigned int first_right_of_diagonal_index
+ = pos_right_of_diagonal[row];
+ double s = *dst_ptr;
+ for (unsigned int j=first_right_of_diagonal_index; j<end_row; ++j)
+ s -= val[j] * dst(cols->colnums[j]);
+
+ *dst_ptr = s * om / val[*rowstart_ptr];
+ };
+ return;
+ }
+
+ // case when we need to get the position
+ // of the first element right of the
+ // diagonal manually for each sweep.
// forward sweep
for (unsigned int row=0; row<n; ++row, ++dst_ptr, ++rowstart_ptr)
{
// $Id: sparse_matrix_matrix.in.h 15011 2007-08-22 16:59:41Z kanschat $
// Version: $Name$
//
-// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2006, 2007 by the deal.II authors
+// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2006, 2007, 2009 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
template void SparseMatrix<S1>::
precondition_SSOR<S2> (Vector<S2> &,
const Vector<S2> &,
- const S1) const;
+ const S1,
+ const std::vector<unsigned int>&) const;
template void SparseMatrix<S1>::
precondition_SOR<S2> (Vector<S2> &,
projects / dealii-svn.git / commitdiff