#include <deal.II/base/vectorization.h>
#include <deal.II/base/subscriptor.h>
+#include <deal.II/lac/diagonal_matrix.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/multigrid/mg_constrained_dofs.h>
#include <deal.II/matrix_free/matrix_free.h>
* the finest mesh or at a certain level in geometric multigrid.
*
* A derived class has to implement apply_add() method as well as
- * compute_diagonal() to fill the protected member inverse_diagonal_entries.
+ * compute_diagonal() to initialize the protected member inverse_diagonal_entries.
* In case of a non-symmetric operator, Tapply_add() should be additionally
* implemented.
*
/**
* Get read access to the inverse diagonal of this operator.
*/
- const LinearAlgebra::distributed::Vector<Number> &get_matrix_diagonal_inverse() const;
+ const DiagonalMatrix<LinearAlgebra::distributed::Vector<Number> > &get_matrix_diagonal_inverse() const;
/**
* Apply the Jacobi preconditioner, which multiplies every element of the
/**
* A vector to store inverse of diagonal elements.
*/
- LinearAlgebra::distributed::Vector<Number> inverse_diagonal_entries;
+ DiagonalMatrix<LinearAlgebra::distributed::Vector<Number> > inverse_diagonal_entries;
private:
Base<dim,Number>::clear ()
{
data = NULL;
- inverse_diagonal_entries.reinit(0);
+ inverse_diagonal_entries.clear();
}
{
(void) col;
Assert (row == col, ExcNotImplemented());
- Assert (inverse_diagonal_entries.size() > 0, ExcNotInitialized());
- return 1.0/inverse_diagonal_entries(row);
+ Assert (inverse_diagonal_entries.m() > 0, ExcNotInitialized());
+ return 1.0/inverse_diagonal_entries(row,row);
}
template <int dim, typename Number>
- const LinearAlgebra::distributed::Vector<Number> &
+ const DiagonalMatrix<LinearAlgebra::distributed::Vector<Number> > &
Base<dim,Number>::get_matrix_diagonal_inverse() const
{
- Assert(inverse_diagonal_entries.size() > 0, ExcNotInitialized());
+ Assert(inverse_diagonal_entries.m() > 0, ExcNotInitialized());
return inverse_diagonal_entries;
}
const LinearAlgebra::distributed::Vector<Number> &src,
const Number omega) const
{
- Assert(inverse_diagonal_entries.size() > 0, ExcNotInitialized());
-
- dst = src;
- dst.scale(inverse_diagonal_entries);
+ Assert(inverse_diagonal_entries.m() > 0, ExcNotInitialized());
+ inverse_diagonal_entries.vmult(dst,src);
dst*= omega;
}
Assert((Base<dim, Number>::data != NULL), ExcNotInitialized());
LinearAlgebra::distributed::Vector<Number> ones;
- Base<dim, Number>::initialize_dof_vector(Base<dim, Number>::inverse_diagonal_entries);
- Base<dim, Number>::initialize_dof_vector(ones);
+ LinearAlgebra::distributed::Vector<Number> &inverse_diagonal_entries = Base<dim, Number>::inverse_diagonal_entries.get_vector();
+ this->initialize_dof_vector(ones);
+ this->initialize_dof_vector(inverse_diagonal_entries);
ones = 1.;
ones.update_ghost_values();
- apply_add(Base<dim, Number>::inverse_diagonal_entries, ones);
+ apply_add(inverse_diagonal_entries, ones);
- const std::vector<unsigned int> &constrained_dofs
- = Base<dim, Number>::data->get_constrained_dofs();
- for (unsigned int i=0; i< constrained_dofs.size(); ++i)
- Base<dim, Number>::inverse_diagonal_entries.local_element(constrained_dofs[i]) = 1.;
+ this->set_constrained_entries_to_one(inverse_diagonal_entries);
- const unsigned int local_size = Base<dim, Number>::inverse_diagonal_entries.local_size();
+ const unsigned int local_size = inverse_diagonal_entries.local_size();
for (unsigned int i=0; i<local_size; ++i)
- Base<dim, Number>::inverse_diagonal_entries.local_element(i)
- =1./Base<dim, Number>::inverse_diagonal_entries.local_element(i);
+ inverse_diagonal_entries.local_element(i)
+ =1./inverse_diagonal_entries.local_element(i);
- Base<dim, Number>::inverse_diagonal_entries.compress(VectorOperation::insert);
- Base<dim, Number>::inverse_diagonal_entries.update_ghost_values();
+ inverse_diagonal_entries.compress(VectorOperation::insert);
+ inverse_diagonal_entries.update_ghost_values();
}
Assert((Base<dim, Number>::data != NULL), ExcNotInitialized());
unsigned int dummy = 0;
- LinearAlgebra::distributed::Vector<Number> &inverse_diagonal_entries = Base<dim,Number>::inverse_diagonal_entries;
+ LinearAlgebra::distributed::Vector<Number> &inverse_diagonal_entries = Base<dim,Number>::inverse_diagonal_entries.get_vector();
this->initialize_dof_vector(inverse_diagonal_entries);
Base<dim,Number>::
data->cell_loop (&LaplaceOperator::local_diagonal_cell,
else
inverse_diagonal_entries.local_element(i) = 1.;
- Base<dim, Number>::inverse_diagonal_entries.compress(VectorOperation::insert);
- Base<dim, Number>::inverse_diagonal_entries.update_ghost_values();
+ inverse_diagonal_entries.compress(VectorOperation::insert);
+ inverse_diagonal_entries.update_ghost_values();
}
//now invert the matrix
ReductionControl control(rhs.size(), 0., 1e-12, false, false);
SolverCG<LinearAlgebra::distributed::Vector<Number> > cg(control);
- typename PreconditionChebyshev<MatrixType, LocalVectorType>::AdditionalData data;
- data.matrix_diagonal_inverse = mass_matrix.get_matrix_diagonal_inverse();
- PreconditionChebyshev<MatrixType, LocalVectorType> preconditioner;
+ typename PreconditionJacobi<MatrixType>::AdditionalData data(0.8);
+ PreconditionJacobi<MatrixType> preconditioner;
preconditioner.initialize(mass_matrix, data);
cg.solve (mass_matrix, vec, rhs, preconditioner);
vec+=inhomogeneities;
//now invert the matrix
ReductionControl control(rhs.size(), 0., 1e-12, false, false);
SolverCG<LinearAlgebra::distributed::Vector<Number> > cg(control);
- typename PreconditionChebyshev<MatrixType, LocalVectorType>::AdditionalData data;
- data.matrix_diagonal_inverse = mass_matrix.get_matrix_diagonal_inverse();
- PreconditionChebyshev<MatrixType, LocalVectorType> preconditioner;
+ typename PreconditionJacobi<MatrixType>::AdditionalData data(0.8);
+ PreconditionJacobi<MatrixType> preconditioner;
preconditioner.initialize(mass_matrix, data);
cg.solve (mass_matrix, vec, rhs, preconditioner);
vec+=inhomogeneities;