const SmartPointer<const InverseMatrix<TrilinosWrappers::SparseMatrix,
PreconditionerMp > > m_inverse;
const PreconditionerA &a_preconditioner;
-
- mutable TrilinosWrappers::Vector tmp;
-
};
:
stokes_matrix (&S),
m_inverse (&Mpinv),
- a_preconditioner (Apreconditioner),
- tmp (stokes_matrix->block(1,1).row_map)
+ a_preconditioner (Apreconditioner)
{}
template <class PreconditionerA, class PreconditionerMp>
const TrilinosWrappers::BlockVector &src) const
{
a_preconditioner.vmult (dst.block(0), src.block(0));
+ TrilinosWrappers::Vector tmp (src.block(1), true);
stokes_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1));
tmp *= -1;
m_inverse->vmult (dst.block(1), tmp);
{
stokes_matrix.clear ();
+ BlockSparsityPattern stokes_sparsity_pattern (2,2);
+
if (trilinos_communicator.MyPID() == 0)
{
BlockCompressedSetSparsityPattern csp (2,2);
DoFTools::make_sparsity_pattern (stokes_dof_handler, coupling, csp);
stokes_constraints.condense (csp);
- BlockSparsityPattern stokes_sparsity_pattern;
stokes_sparsity_pattern.copy_from (csp);
-
- stokes_matrix.reinit (stokes_partitioner, stokes_sparsity_pattern);
}
- stokes_matrix.collect_sizes();
+ stokes_matrix.reinit (stokes_partitioner, stokes_sparsity_pattern);
}
{
Mp_preconditioner.reset ();
stokes_preconditioner_matrix.clear ();
+ BlockSparsityPattern stokes_preconditioner_sparsity_pattern (2,2);
+
if (trilinos_communicator.MyPID() == 0)
{
BlockCompressedSetSparsityPattern csp (2,2);
DoFTools::make_sparsity_pattern (stokes_dof_handler, coupling, csp);
stokes_constraints.condense (csp);
- BlockSparsityPattern stokes_preconditioner_sparsity_pattern;
stokes_preconditioner_sparsity_pattern.copy_from (csp);
-
- stokes_preconditioner_matrix.reinit (stokes_partitioner,
- stokes_preconditioner_sparsity_pattern);
}
- stokes_preconditioner_matrix.collect_sizes();
+ stokes_preconditioner_matrix.reinit (stokes_partitioner,
+ stokes_preconditioner_sparsity_pattern);
+
}
+ std::cout << " bla " << std::flush;
temperature_partitioner = Epetra_Map (n_T, 0, trilinos_communicator);
{
temperature_stiffness_matrix.clear ();
temperature_matrix.clear ();
+ SparsityPattern temperature_sparsity_pattern;
+
if (trilinos_communicator.MyPID() == 0)
{
CompressedSetSparsityPattern csp (n_T, n_T);
DoFTools::make_sparsity_pattern (temperature_dof_handler, csp);
temperature_constraints.condense (csp);
- SparsityPattern temperature_sparsity_pattern;
temperature_sparsity_pattern.copy_from (csp);
-
- temperature_matrix.reinit (temperature_partitioner,
- temperature_sparsity_pattern);
- temperature_mass_matrix.reinit (temperature_partitioner,
- temperature_sparsity_pattern);
- temperature_stiffness_matrix.reinit (temperature_partitioner,
- temperature_sparsity_pattern);
}
- temperature_matrix.compress();
- temperature_mass_matrix.compress();
- temperature_stiffness_matrix.compress();
+
+ temperature_matrix.reinit (temperature_partitioner,
+ temperature_sparsity_pattern);
+ temperature_mass_matrix.reinit (temperature_partitioner,
+ temperature_sparsity_pattern);
+ temperature_stiffness_matrix.reinit (temperature_partitioner,
+ temperature_sparsity_pattern);
}
stokes_solution.reinit (stokes_partitioner);
cell->get_dof_indices (local_dof_indices);
stokes_constraints.distribute_local_to_global (local_matrix,
- local_dof_indices,
- stokes_preconditioner_matrix);
+ local_dof_indices,
+ stokes_preconditioner_matrix);
}
stokes_preconditioner_matrix.compress();
old_old_temperature_solution = old_temperature_solution;
old_temperature_solution = temperature_solution;
}
- while (time <= 100);
+ while (time <= 1);
}
(void)argv;
#endif
+ sleep (20);
try
{
deallog.depth_console (0);
int, int,
<< "You tried to access element (" << arg1
<< "/" << arg2 << ")"
- << " of a sparse matrix, but it appears to not "
+ << " of a sparse matrix, but it appears to not"
<< " exist in the Trilinos sparsity pattern.");
+ private:
/**
* Epetra Trilinos mapping of the
* matrix rows that
*/
Epetra_Map col_map;
- /**
- * A sparse matrix object in
- * Trilinos to be used for
- * finite element based problems
- * which allows for assembling into
- * non-local elements.
- * The actual type, a sparse
- * matrix, is set in the constructor.
- */
- std::auto_ptr<Epetra_FECrsMatrix> matrix;
-
- protected:
/**
* Trilinos doesn't allow to mix additions
* to matrix entries and overwriting
* or an insertion.
*/
Epetra_CombineMode last_action;
+
+ public:
+ /**
+ * A sparse matrix object in
+ * Trilinos to be used for
+ * finite element based problems
+ * which allows for assembling into
+ * non-local elements.
+ * The actual type, a sparse
+ * matrix, is set in the constructor.
+ */
+ std::auto_ptr<Epetra_FECrsMatrix> matrix;
+
};
* as to only allow the actual vector
* class to create it.
*/
- VectorReference (Vector &vector,
+ VectorReference (Vector &vector,
const unsigned int index);
public:
<< " are stored locally and can be accessed.");
+ private:
/**
* An Epetra map used to map vector data
* accross multiple processes. This is
*/
Epetra_Map map;
- /**
- * An Epetra distibuted vector type.
- * Requires an existing Epetra_Map for
- * storing data.
- */
- std::auto_ptr<Epetra_FEVector> vector;
-
-
- private:
- /**
+ /**
* Trilinos doesn't allow to mix additions
* to matrix entries and overwriting
* them (to make synchronisation of
*/
Epetra_CombineMode last_action;
+ public:
+ /**
+ * An Epetra distibuted vector type.
+ * Requires an existing Epetra_Map for
+ * storing data.
+ * TODO: Should become private at some
+ * point.
+ */
+ std::auto_ptr<Epetra_FEVector> vector;
+
+
/**
* Make the reference class a friend.
*/
{
const unsigned int n_block_rows = input_maps.size();
- Assert (n_block_rows == block_sparsity_pattern.n_block_rows(),
- ExcDimensionMismatch (n_block_rows,
- block_sparsity_pattern.n_block_rows()));
- Assert (n_block_rows == block_sparsity_pattern.n_block_cols(),
- ExcDimensionMismatch (n_block_rows,
- block_sparsity_pattern.n_block_cols()));
+ if (input_maps[0].Comm().MyPID()==0)
+ {
+ Assert (n_block_rows == block_sparsity_pattern.n_block_rows(),
+ ExcDimensionMismatch (n_block_rows,
+ block_sparsity_pattern.n_block_rows()));
+ Assert (n_block_rows == block_sparsity_pattern.n_block_cols(),
+ ExcDimensionMismatch (n_block_rows,
+ block_sparsity_pattern.n_block_cols()));
+ }
// Call the other basic reinit function ...
reinit (n_block_rows, n_block_rows);
this->block(r,c).reinit(input_maps[r],input_maps[c],
block_sparsity_pattern.block(r,c));
}
+
+ collect_sizes();
}
dealii_block_sparse_matrix.block(r,c),
drop_tolerance);
}
+
+ collect_sizes();
}
void PreconditionAMG::vmult (dealii::Vector<double> &dst,
const dealii::Vector<double> &src) const
{
- Assert (Map->SameAs(Matrix->row_map),
+ Assert (Map->SameAs(Matrix->matrix->RowMap()),
ExcMessage("The sparse matrix given to the preconditioner uses "
"a map that is not compatible. Check ML preconditioner "
"setup."));
- Assert (Map->SameAs(Matrix->col_map),
+ Assert (Map->SameAs(Matrix->matrix->RowMap()),
ExcMessage("The sparse matrix given to the preconditioner uses "
"a map that is not compatible. Check ML preconditioner "
"setup."));
row_map (0, 0, Epetra_SerialComm()),
#endif
col_map (row_map),
+ last_action (Insert),
matrix (std::auto_ptr<Epetra_FECrsMatrix>
- (new Epetra_FECrsMatrix(Copy, row_map, 0))),
- last_action (Insert)
+ (new Epetra_FECrsMatrix(Copy, row_map, 0)))
{}
SparseMatrix::SparseMatrix (const Epetra_Map &InputMap,
:
row_map (InputMap),
col_map (row_map),
+ last_action (Insert),
matrix (std::auto_ptr<Epetra_FECrsMatrix>
(new Epetra_FECrsMatrix(Copy, row_map,
- int(n_max_entries_per_row), false))),
- last_action (Insert)
+ int(n_max_entries_per_row), false)))
{}
SparseMatrix::SparseMatrix (const Epetra_Map &InputMap,
:
row_map (InputMap),
col_map (row_map),
+ last_action (Insert),
matrix (std::auto_ptr<Epetra_FECrsMatrix>
(new Epetra_FECrsMatrix(Copy, row_map,
(int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
- true))),
- last_action (Insert)
+ true)))
{}
SparseMatrix::SparseMatrix (const Epetra_Map &InputRowMap,
:
row_map (InputRowMap),
col_map (InputColMap),
+ last_action (Insert),
matrix (std::auto_ptr<Epetra_FECrsMatrix>
(new Epetra_FECrsMatrix(Copy, row_map, col_map,
- int(n_max_entries_per_row), false))),
- last_action (Insert)
+ int(n_max_entries_per_row), false)))
{}
SparseMatrix::SparseMatrix (const Epetra_Map &InputRowMap,
:
row_map (InputRowMap),
col_map (InputColMap),
+ last_action (Insert),
matrix (std::auto_ptr<Epetra_FECrsMatrix>
(new Epetra_FECrsMatrix(Copy, row_map, col_map,
(int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
- true))),
- last_action (Insert)
+ true)))
{}
{
unsigned int n_rows = sparsity_pattern.n_rows();
- Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
- ExcDimensionMismatch (matrix->NumGlobalRows(),
- sparsity_pattern.n_rows()));
-
- // Trilinos seems to have a bug for
- // rectangular matrices at this point,
- // so do not check for consistent
- // column numbers here.
- //
- //
- // this bug is filed in the Sandia
- // bugzilla under #4123 and should be
- // fixed for version 9.0
-// Assert (matrix->NumGlobalCols() == (int)sparsity_pattern.n_cols(),
-// ExcDimensionMismatch (matrix->NumGlobalCols(),
-// sparsity_pattern.n_cols()));
-
- std::vector<int> n_entries_per_row(n_rows);
-
- for (unsigned int row=0; row<n_rows; ++row)
- n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
-
- std::vector<double> values;
- std::vector<int> row_indices;
-
- for (unsigned int row=0; row<n_rows; ++row)
+ if (row_map.Comm().MyPID() == 0)
{
- const int row_length = sparsity_pattern.row_length(row);
- row_indices.resize (row_length, 0);
- values.resize (row_length, 0.);
-
- for (int col=0; col< row_length; ++col)
- row_indices[col] = sparsity_pattern.column_number (row, col);
-
- matrix->InsertGlobalValues(row, row_length,
- &values[0], &row_indices[0]);
+ Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
+ ExcDimensionMismatch (matrix->NumGlobalRows(),
+ sparsity_pattern.n_rows()));
+
+ // Trilinos seems to have a bug for
+ // rectangular matrices at this point,
+ // so do not check for consistent
+ // column numbers here.
+ //
+ //
+ // this bug is filed in the Sandia
+ // bugzilla under #4123 and should be
+ // fixed for version 9.0
+// Assert (matrix->NumGlobalCols() == (int)sparsity_pattern.n_cols(),
+// ExcDimensionMismatch (matrix->NumGlobalCols(),
+// sparsity_pattern.n_cols()));
+
+ std::vector<int> n_entries_per_row(n_rows);
+
+ for (unsigned int row=0; row<n_rows; ++row)
+ n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
+
+ std::vector<double> values;
+ std::vector<int> row_indices;
+
+ for (unsigned int row=0; row<n_rows; ++row)
+ {
+ const int row_length = sparsity_pattern.row_length(row);
+ row_indices.resize (row_length, 0);
+ values.resize (row_length, 0.);
+
+ for (int col=0; col< row_length; ++col)
+ row_indices[col] = sparsity_pattern.column_number (row, col);
+
+ matrix->InsertGlobalValues(row, row_length,
+ &values[0], &row_indices[0]);
+ }
}
+ compress();
// In the end, the matrix needs to
// be compressed in order to be
// really ready. However, that is
unsigned int n_rows = sparsity_pattern.n_rows();
+ if (row_map.Comm().MyPID() == 0)
+ {
Assert (input_row_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(),
ExcDimensionMismatch (input_row_map.NumGlobalElements(),
sparsity_pattern.n_rows()));
Assert (input_col_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(),
ExcDimensionMismatch (input_col_map.NumGlobalElements(),
sparsity_pattern.n_cols()));
+ }
row_map = input_row_map;
col_map = input_col_map;
const int n_row_entries = index;
matrix->InsertGlobalValues(row, n_row_entries,
- &values[0], &row_indices[0]);
+ &values[0], &row_indices[0]);
}
}
{
Assert (&src != &dst, ExcSourceEqualsDestination());
- Assert (col_map.SameAs(src.map),
+ Assert (col_map.SameAs(src.vector->Map()),
ExcMessage ("Column map of matrix does not fit with vector map!"));
- Assert (row_map.SameAs(dst.map),
+ Assert (row_map.SameAs(dst.vector->Map()),
ExcMessage ("Row map of matrix does not fit with vector map!"));
if (!matrix->Filled())
{
Assert (&src != &dst, ExcSourceEqualsDestination());
- Assert (row_map.SameAs(src.map),
+ Assert (row_map.SameAs(src.vector->Map()),
ExcMessage ("Row map of matrix does not fit with vector map!"));
- Assert (col_map.SameAs(dst.map),
+ Assert (col_map.SameAs(dst.vector->Map()),
ExcMessage ("Column map of matrix does not fit with vector map!"));
if (!matrix->Filled())
TrilinosScalar
SparseMatrix::matrix_norm_square (const Vector &v) const
{
- Vector tmp(v.map);
+ Vector tmp(v,true);
vmult (tmp, v);
return tmp*v;
}
TrilinosScalar
SparseMatrix::matrix_scalar_product (const Vector &u,
- const Vector &v) const
+ const Vector &v) const
{
- Vector tmp(v.map);
+ Vector tmp(v,true);
vmult (tmp, v);
return u*tmp;
}
#else
map (0,0,Epetra_SerialComm()),
#endif
+ last_action (Insert),
vector(std::auto_ptr<Epetra_FEVector>
- (new Epetra_FEVector(map))),
- last_action (Insert)
+ (new Epetra_FEVector(map)))
{}
Vector::Vector (unsigned int GlobalSize, Epetra_Comm &Comm)
:
map (GlobalSize, 0, Comm),
+ last_action (Insert),
vector (std::auto_ptr<Epetra_FEVector>
- (new Epetra_FEVector(map))),
- last_action (Insert)
+ (new Epetra_FEVector(map)))
{}
Vector::Vector (const Epetra_Map &InputMap)
:
map (InputMap),
+ last_action (Insert),
vector (std::auto_ptr<Epetra_FEVector>
- (new Epetra_FEVector(map))),
- last_action (Insert)
+ (new Epetra_FEVector(map)))
{}
const bool fast)
:
map (v.map),
+ last_action (Insert),
vector(std::auto_ptr<Epetra_FEVector>
- (new Epetra_FEVector(v.map,!fast))),
- last_action (Insert)
+ (new Epetra_FEVector(v.map,!fast)))
{}