*/
~BlockSparseMatrix ();
- /**
+ /**
* Pseudo copy operator only copying
* empty objects. The sizes of the block
* matrices need to be the same.
class BlockVector;
}
class BlockVector;
+ class BlockSparseMatrix;
namespace MPI
*/
void reinit (const unsigned int num_blocks);
- /**
+ /**
+ * This reinit function is
+ * meant to be used for
+ * parallel calculations where
+ * some non-local data has to
+ * be used. The typical
+ * situation where one needs
+ * this function is the call of
+ * the
+ * FEValues<dim>::get_function_values
+ * function (or of some
+ * derivatives) in
+ * parallel. Since it is
+ * usually faster to retrieve
+ * the data in advance, this
+ * function can be called
+ * before the assembly forks
+ * out to the different
+ * processors. What this
+ * function does is the
+ * following: It takes the
+ * information in the columns
+ * of the given matrix and
+ * looks which data couples
+ * between the different
+ * processors. That data is
+ * then queried from the input
+ * vector. Note that you should
+ * not write to the resulting
+ * vector any more, since the
+ * some data can be stored
+ * several times on different
+ * processors, leading to
+ * unpredictable results. In
+ * particular, such a vector
+ * cannot be used for
+ * matrix-vector products as
+ * for example done during the
+ * solution of linear systems.
+ */
+ void do_data_exchange (const TrilinosWrappers::BlockSparseMatrix &m,
+ const BlockVector &v);
+
+ /**
* Compresses all the components
* after assembling together all
* elements.
/**
* Returns the state of the
- * matrix, i.e., whether
+ * vector, i.e., whether
* compress() needs to be
* called after an operation
* requiring data
friend class SolverBase;
friend class SolverSaddlePoint;
- protected:
+ //protected:
/**
* This is a pointer to the
* preconditioner object.
* AdditionalData structure for details), and a parameter
* <tt>overlap</tt> that determines if and how much overlap there
* should be between the matrix partitions on the various MPI
- * processes. The default settings are 1 for the relaxation parameter,
- * 0 for the diagonal augmentation and 0 for the overlap.
+ * processes. The default settings are 0 for the additional fill-in, 0
+ * for the absolute augmentation tolerance, 1 for the relative
+ * augmentation tolerance, 0 for the overlap.
*
* Note that a parallel applicatoin of the IC preconditioner is
* actually a block-Jacobi preconditioner with block size equal to the
* AdditionalData structure for details), and a parameter
* <tt>overlap</tt> that determines if and how much overlap there
* should be between the matrix partitions on the various MPI
- * processes. The default settings are 1 for the relaxation parameter,
- * 0 for the diagonal augmentation and 0 for the overlap.
+ * processes. The default settings are 0 for the additional fill-in, 0
+ * for the absolute augmentation tolerance, 1 for the relative
+ * augmentation tolerance, 0 for the overlap.
*
* Note that a parallel applicatoin of the ILU preconditioner is
* actually a block-Jacobi preconditioner with block size equal to the
void initialize (const SparseMatrix &matrix,
const AdditionalData &additional_data = AdditionalData());
};
+
+
+
+/**
+ * A wrapper class for a sparse direct LU decomposition on parallel
+ * blocks for Trilinos matrices. When run in serial, this corresponds
+ * to a direct solve on the matrix.
+ *
+ * The AdditionalData data structure allows to set preconditioner
+ * options.
+ *
+ * Note that a parallel applicatoin of the block direct solve
+ * preconditioner is actually a block-Jacobi preconditioner with block
+ * size equal to the local matrix size. Spoken more technically, this
+ * parallel operation is an <a
+ * href="http://en.wikipedia.org/wiki/Additive_Schwarz_method">additive
+ * Schwarz method</a> with an <em>exact solve</em> as inner solver,
+ * based on the (outer) parallel partitioning.
+ *
+ * @ingroup TrilinosWrappers
+ * @ingroup Preconditioners
+ * @author Martin Kronbichler, 2008
+ */
+ class PreconditionBlockDirect : public PreconditionBase
+ {
+ public:
+ /**
+ * Standardized data struct to
+ * pipe additional parameters
+ * to the preconditioner.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor.
+ */
+ AdditionalData (const unsigned int overlap = 0);
+
+ /**
+ * Block direct parameters and overlap.
+ */
+ unsigned int overlap;
+ };
+
+ /**
+ * Initialize function. Takes
+ * the matrix which is used to
+ * form the preconditioner, and
+ * additional flags if there
+ * are any.
+ */
+ void initialize (const SparseMatrix &matrix,
+ const AdditionalData &additional_data = AdditionalData());
+ };
}
//---------------------------------------------------------------------------
#include <lac/trilinos_block_vector.h>
+#include <lac/trilinos_block_sparse_matrix.h>
#ifdef DEAL_II_USE_TRILINOS
+ void
+ BlockVector::do_data_exchange (const TrilinosWrappers::BlockSparseMatrix &m,
+ const BlockVector &v)
+ {
+ Assert (m.n_block_rows() == v.n_blocks(),
+ ExcDimensionMismatch(m.n_block_rows(),v.n_blocks()));
+ Assert (m.n_block_cols() == v.n_blocks(),
+ ExcDimensionMismatch(m.n_block_cols(),v.n_blocks()));
+
+ if (v.n_blocks() != n_blocks())
+ {
+ block_indices = v.get_block_indices();
+ components.resize(v.n_blocks());
+ }
+
+ for (unsigned int i=0; i<this->n_blocks(); ++i)
+ components[i].do_data_exchange(m.block(i,i), v.block(i));
+
+ collect_sizes();
+ }
+
+
+
void
BlockVector::compress ()
{
int ierr;
Teuchos::ParameterList parameter_list;
+ parameter_list.set ("relaxation: sweeps", 1);
parameter_list.set ("relaxation: type", "Jacobi");
parameter_list.set ("relaxation: damping factor", additional_data.omega);
parameter_list.set ("relaxation: min diagonal value",
int ierr;
Teuchos::ParameterList parameter_list;
+ parameter_list.set ("relaxation: sweeps", 1);
parameter_list.set ("relaxation: type", "symmetric Gauss-Seidel");
parameter_list.set ("relaxation: damping factor", additional_data.omega);
parameter_list.set ("relaxation: min diagonal value",
int ierr;
Teuchos::ParameterList parameter_list;
+ parameter_list.set ("relaxation: sweeps", 1);
parameter_list.set ("relaxation: type", "Gauss-Seidel");
parameter_list.set ("relaxation: damping factor", additional_data.omega);
parameter_list.set ("relaxation: min diagonal value",
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
}
+
+
+/* ---------------------- PreconditionBlockDirect --------------------- */
+
+ PreconditionBlockDirect::AdditionalData::
+ AdditionalData (const unsigned int overlap)
+ :
+ overlap (overlap)
+ {}
+
+
+
+ void
+ PreconditionBlockDirect::initialize (const SparseMatrix &matrix,
+ const AdditionalData &additional_data)
+ {
+ preconditioner.release();
+
+ preconditioner = Teuchos::rcp (Ifpack().Create ("Amesos", &*matrix.matrix,
+ additional_data.overlap));
+ Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
+ "preconditioner"));
+
+ int ierr;
+
+ Teuchos::ParameterList parameter_list;
+ parameter_list.set ("schwarz: combine mode", "Add");
+
+ ierr = preconditioner->SetParameters(parameter_list);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ ierr = preconditioner->Initialize();
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ ierr = preconditioner->Compute();
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
}
DEAL_II_NAMESPACE_CLOSE
return *this;
}
+
+
+ void
+ Vector::do_data_exchange (const TrilinosWrappers::SparseMatrix &m,
+ const Vector &v)
+ {
+ Assert (m.matrix->Filled() == true,
+ ExcMessage ("Matrix is not compressed. "
+ "Cannot find exchange information!"));
+ Assert (v.vector->Map().UniqueGIDs() == true,
+ ExcMessage ("The input vector has overlapping data, "
+ "which is not allowed."));
+
+ if (vector->Map().SameAs(m.matrix->ColMap()) == false)
+ {
+ map = m.matrix->ColMap();
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map));
+ }
+
+ Epetra_Import data_exchange (vector->Map(), v.vector->Map());
+ const int ierr = vector->Import(*v.vector, data_exchange, Insert);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ last_action = Insert;
+ }
+
} /* end of namespace MPI */