#include <base/config.h>
#include <base/table.h>
#include <lac/block_matrix_base.h>
+#include <lac/block_sparsity_pattern.h>
#include <lac/trilinos_sparse_matrix.h>
#include <lac/trilinos_block_vector.h>
#include <lac/exceptions.h>
* assumes that a quadratic block
* matrix is generated.
*/
- //void reinit (const std::vector<Epetra_Map> &input_maps);
+ void reinit (const std::vector<Epetra_Map> &input_maps,
+ const BlockSparsityPattern &block_sparsity_pattern);
/**
* This function calls the compress()
inline
- BlockVector::BlockVector (const unsigned int num_blocks = 0)
+ BlockVector::BlockVector (const unsigned int num_blocks)
{
reinit (num_blocks);
}
* TrilinosWrappers::SparseMatrix.
*/
void initialize (const SparseMatrix &matrix,
- const std::vector<double> &null_space,
- const unsigned int null_space_dimension,
- const bool higher_order_elements,
- const bool elliptic,
- const bool output_details);
+ const bool elliptic = true,
+ const bool higher_order_elements = false,
+ const std::vector<double> &null_space = std::vector<double>(),
+ const unsigned int null_space_dimension = 1,
+ const bool output_details = false);
/**
* Let Trilinos compute a
* inefficient.
*/
void initialize (const dealii::SparseMatrix<double> &matrix,
- const std::vector<double> &null_space,
- const unsigned int null_space_dimension,
- const bool higher_order_elements,
- const bool elliptic,
- const bool output_details);
+ const bool elliptic = true,
+ const bool higher_order_elements = false,
+ const std::vector<double> &null_space = std::vector<double>(),
+ const unsigned int null_space_dimension = 1,
+ const bool output_details = false);
/**
* This function can be used
*/
void reinit (const Epetra_Map &input_map,
const ::dealii::SparseMatrix<double> &deal_ii_sparse_matrix,
- const double drop_tolerance=1e-15);
+ const double drop_tolerance=1e-13);
/**
* This function is similar to the
/**
* Copy constructor. Sets the dimension
* to that of the given vector and uses
- * the map of that vector, and
- * copies all elements.
+ * the map of that vector, but
+ * does not copies any element. Instead,
+ * the memory will remain untouched
+ * in case <tt>fast</tt> is false and
+ * initialized with zero otherwise.
*/
- Vector (const Vector &v);
+ Vector (const Vector &v,
+ const bool fast = false);
/**
* Destructor
* copies the vector v to the current
* one.
*/
- void reinit (const Vector &v);
+ void reinit (const Vector &v,
+ const bool fast = false);
/**
* Release all memory and return
*/
void swap (Vector &v);
+ /**
+ * Estimate for the memory
+ * consumption (not implemented
+ * for this class).
+ */
+ unsigned int memory_consumption () const;
+
/**
* Exception
*/
{
// if the vectors have different sizes,
// then first resize the present one
- reinit (v);
+ if (!map.SameAs(v.map))
+ {
+ map = v.map;
+ vector = std::auto_ptr<Epetra_FEVector>
+ (new Epetra_FEVector(*v.vector));
+ }
+ else
+ *vector = *v.vector;
+
+ last_action = Insert;
return *this;
}
}
+
void
BlockSparseMatrix::
reinit (const unsigned int n_block_rows,
}
}
-
+
+
+ void
+ BlockSparseMatrix::
+ reinit (const std::vector<Epetra_Map> &input_maps,
+ const BlockSparsityPattern &block_sparsity_pattern)
+ {
+ const unsigned int n_block_rows = input_maps.size();
+
+ // Call the other reinit function ...
+ reinit (n_block_rows, n_block_rows);
+
+ // ... and now assign the correct
+ // blocks.
+
+ for (unsigned int r=0; r<this->n_block_rows(); ++r)
+ for (unsigned int c=0; c<this->n_block_cols(); ++c)
+ {
+ this->block(r,c).reinit(input_maps[r],input_maps[c],
+ block_sparsity_pattern.block(r,c));
+ }
+ }
+
+
+
void
BlockSparseMatrix::compress()
{
void
PreconditionAMG::
initialize (const SparseMatrix &matrix,
- const std::vector<double> &null_space,
- const unsigned int null_space_dimension,
const bool elliptic,
const bool higher_order_elements,
+ const std::vector<double> &null_space,
+ const unsigned int null_space_dimension,
const bool output_details)
{
const unsigned int n_rows = matrix.m();
-
- if (!Matrix)
- {
- Matrix = boost::shared_ptr<SparseMatrix>
- (const_cast<SparseMatrix*>(&matrix));
- Map = boost::shared_ptr<Epetra_Map>
- (const_cast<Epetra_Map*>(&matrix.row_map));
- }
// Build the AMG preconditioner.
Teuchos::ParameterList parameter_list;
parameter_list.set("aggregation: block scaling", true);
}
+ parameter_list.set("aggregation: threshold", 1e-12);
+
if (output_details)
parameter_list.set("ML output", 10);
else
if (higher_order_elements)
parameter_list.set("aggregation: type", "MIS");
- Assert (n_rows * null_space_dimension == null_space.size(),
- ExcDimensionMismatch(n_rows * null_space_dimension,
- null_space.size()));
-
if (null_space_dimension > 1)
{
+ Assert (n_rows * null_space_dimension == null_space.size(),
+ ExcDimensionMismatch(n_rows * null_space_dimension,
+ null_space.size()));
+
parameter_list.set("null space: type", "pre-computed");
parameter_list.set("null space: dimension", int(null_space_dimension));
parameter_list.set("null space: vectors", (double *)&null_space[0]);
void
PreconditionAMG::
initialize (const dealii::SparseMatrix<double> &deal_ii_sparse_matrix,
- const std::vector<double> &null_space,
- const unsigned int null_space_dimension,
const bool elliptic,
const bool higher_order_elements,
+ const std::vector<double> &null_space,
+ const unsigned int null_space_dimension,
const bool output_details)
{
const unsigned int n_rows = deal_ii_sparse_matrix.m();
Matrix->reinit (*Map, deal_ii_sparse_matrix);
Matrix->compress();
- initialize (*Matrix, null_space, null_space_dimension, elliptic,
- higher_order_elements, output_details);
+ initialize (*Matrix, elliptic, higher_order_elements, null_space,
+ null_space_dimension, output_details);
}
void PreconditionAMG::vmult (Vector &dst,
const Vector &src) const
{
- const int ierr = multigrid_operator->ApplyInverse (*src.vector, *dst.vector);
+ const int ierr = multigrid_operator->ApplyInverse (*src.vector,
+ *dst.vector);
Assert (ierr == 0, ExcTrilinosError(ierr));
}
"a map that is not compatible. Check ML preconditioner "
"setup."));
- Epetra_Vector LHS (View, *Map, dst.begin());
- Epetra_Vector RHS (View, *Map, const_cast<double*>(src.begin()));
+ Epetra_Vector LHS (View, multigrid_operator->OperatorDomainMap(),
+ dst.begin());
+ Epetra_Vector RHS (View, multigrid_operator->OperatorDomainMap(),
+ const_cast<double*>(src.begin()));
const int res = multigrid_operator->ApplyInverse (RHS, LHS);
{}
- Vector::Vector (const Vector &v)
+ Vector::Vector (const Vector &v,
+ const bool fast)
:
map (v.map),
vector(std::auto_ptr<Epetra_FEVector>
- (new Epetra_FEVector(*(v.vector)))),
+ (new Epetra_FEVector(v.map,!fast))),
last_action (Insert)
{}
void
- Vector::reinit (const Vector &v)
+ Vector::reinit (const Vector &v,
+ const bool fast)
{
vector.reset();
if (!map.SameAs(v.map))
map = v.map;
- vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(*v.vector));
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(v.map,!fast));
last_action = Insert;
}
vector = tmp;
}
+
+ unsigned int
+ Vector::memory_consumption () const
+ {
+ AssertThrow(false, ExcNotImplemented() );
+ return 0;
+ }
+
}
DEAL_II_NAMESPACE_CLOSE
#include <lac/petsc_block_vector.h>
#include <lac/petsc_parallel_vector.h>
#include <lac/petsc_parallel_block_vector.h>
+#include <lac/trilinos_vector.h>
+#include <lac/trilinos_block_vector.h>
DEAL_II_NAMESPACE_OPEN
template class GrowingVectorMemory<PETScWrappers::MPI::BlockVector>;
#endif
+#ifdef DEAL_II_USE_TRILINOS
+ template class VectorMemory<TrilinosWrappers::Vector>;
+ template class GrowingVectorMemory<TrilinosWrappers::Vector>;
+
+ template class VectorMemory<TrilinosWrappers::BlockVector>;
+ template class GrowingVectorMemory<TrilinosWrappers::BlockVector>;
+#endif
+
DEAL_II_NAMESPACE_CLOSE