const std::vector<unsigned int> &n_entries_per_row);
/**
- * Destructor. Made virtual so that one
- * can use pointers to this class.
+ * This function is similar to the
+ * one above, but it now takes two
+ * different Epetra maps for rows
+ * and columns. This interface is
+ * meant to be used for generating
+ * rectangular matrices, where one
+ * map specifies the parallel
+ * distribution of rows and the other
+ * the one of the columns. This is
+ * in contrast to the first
+ * constructor, where the same map
+ * is used for both the number of
+ * rows and the number of columns.
+ * The number
+ * of columns per row is specified
+ * by the maximum number of entries.
*/
- virtual ~SparseMatrix ();
+ SparseMatrix (const Epetra_Map &InputRowMap,
+ const Epetra_Map &InputColMap,
+ const unsigned int n_max_entries_per_row);
+
/**
- * This operator assigns a scalar to a
- * matrix. Since this does usually not
- * make much sense (should we set all
- * matrix entries to this value? Only
- * the nonzero entries of the sparsity
- * pattern?), this operation is only
- * allowed if the actual value to be
- * assigned is zero. This operator only
- * exists to allow for the obvious
- * notation <tt>matrix=0</tt>, which
- * sets all elements of the matrix to
- * zero, but keeps the sparsity pattern
- * previously used.
+ * This function is similar to the
+ * one above, but it now takes two
+ * different Epetra maps for rows
+ * and columns. This interface is
+ * meant to be used for generating
+ * rectangular matrices, where one
+ * map specifies the parallel
+ * distribution of rows and the other
+ * the one of the columns. The
+ * vector n_entries_per_row specifies
+ * the number of entries in each
+ * row of the newly generated matrix.
*/
- SparseMatrix &
- operator = (const double d);
+ SparseMatrix (const Epetra_Map &InputRowMap,
+ const Epetra_Map &InputColMap,
+ const std::vector<unsigned int> &n_entries_per_row);
/**
+ * Destructor. Made virtual so that one
+ * can use pointers to this class.
+ */
+ virtual ~SparseMatrix ();
+
+ /**
* This function initializes the
* Trilinos matrix by attaching all
* the elements to the sparsity
* e.g. when the grid is refined.
*/
void reinit (const Epetra_Map &input_map,
+ const CompressedSparsityPattern &sparsity_pattern);
+
+ /**
+ * This function is similar to the
+ * other initialization function above,
+ * but now also reassigns the matrix
+ * rows and columns according to
+ * two user-supplied Epetra maps.
+ * To be used e.g. for rectangular
+ * matrices after remeshing.
+ */
+ void reinit (const Epetra_Map &input_row_map,
+ const Epetra_Map &input_col_map,
const CompressedSparsityPattern &sparsity_pattern);
/**
*/
void clear ();
+ /**
+ * This operator assigns a scalar to a
+ * matrix. Since this does usually not
+ * make much sense (should we set all
+ * matrix entries to this value? Only
+ * the nonzero entries of the sparsity
+ * pattern?), this operation is only
+ * allowed if the actual value to be
+ * assigned is zero. This operator only
+ * exists to allow for the obvious
+ * notation <tt>matrix=0</tt>, which
+ * sets all elements of the matrix to
+ * zero, but keeps the sparsity pattern
+ * previously used.
+ */
+ SparseMatrix &
+ operator = (const double d);
+
/**
* Set the element (<i>i,j</i>)
* to @p value.
/**
* Exception
*/
- DeclException4 (ExcAccessToNonlocalElement,
+ DeclException4 (ExcAccessToNonLocalElement,
int, int, int, int,
<< "You tried to access element (" << arg1
<< "/" << arg2 << ")"
<< " are stored locally and can be accessed.");
/**
- * The Epetra Trilinos mapping that
+ * Exception
+ */
+ DeclException2 (ExcAccessToNonPresentElement,
+ int, int,
+ << "You tried to access element (" << arg1
+ << "/" << arg2 << ")"
+ << " of a sparse matrix, but it appears to not "
+ << " exist in the Trilinos sparsity pattern.");
+
+ /**
+ * Epetra Trilinos mapping of the
+ * matrix rows that
+ * assigns parts of the matrix to
+ * the individual processes.
+ * TODO: is it possible to only use
+ * a pointer?
+ */
+ Epetra_Map row_map;
+
+ /**
+ * Pointer to the user-supplied
+ * Epetra Trilinos mapping of the
+ * matrix columns that
* assigns parts of the matrix to
* the individual processes.
*/
- Epetra_Map map;
+ Epetra_Map col_map;
/**
* A sparse matrix object in
inline
const_iterator::
- const_iterator(const SparseMatrix *matrix,
+ const_iterator(const SparseMatrix *matrix,
const unsigned int row,
const unsigned int index)
:
SparseMatrix::SparseMatrix ()
:
#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
- map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)),
+ row_map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD)),
#else
- map (0,0,Epetra_SerialComm()),
+ row_map (0, 0, Epetra_SerialComm()),
#endif
+ col_map (row_map),
matrix (std::auto_ptr<Epetra_FECrsMatrix>
- (new Epetra_FECrsMatrix(Copy, map, 0))),
+ (new Epetra_FECrsMatrix(Copy, row_map, 0))),
last_action (Insert)
{}
-
SparseMatrix::SparseMatrix (const Epetra_Map &InputMap,
const unsigned int n_max_entries_per_row)
:
- map (InputMap),
+ row_map (InputMap),
+ col_map (row_map),
matrix (std::auto_ptr<Epetra_FECrsMatrix>
- (new Epetra_FECrsMatrix(Copy, map,
+ (new Epetra_FECrsMatrix(Copy, row_map,
int(n_max_entries_per_row), false))),
last_action (Insert)
{}
SparseMatrix::SparseMatrix (const Epetra_Map &InputMap,
const std::vector<unsigned int> &n_entries_per_row)
:
- map (InputMap),
+ row_map (InputMap),
+ col_map (row_map),
+ matrix (std::auto_ptr<Epetra_FECrsMatrix>
+ (new Epetra_FECrsMatrix(Copy, row_map,
+ (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
+ false))),
+ last_action (Insert)
+ {}
+
+ SparseMatrix::SparseMatrix (const Epetra_Map &InputRowMap,
+ const Epetra_Map &InputColMap,
+ const unsigned int n_max_entries_per_row)
+ :
+ row_map (InputRowMap),
+ col_map (InputColMap),
+ matrix (std::auto_ptr<Epetra_FECrsMatrix>
+ (new Epetra_FECrsMatrix(Copy, row_map, col_map,
+ int(n_max_entries_per_row), false))),
+ last_action (Insert)
+ {}
+
+ SparseMatrix::SparseMatrix (const Epetra_Map &InputRowMap,
+ const Epetra_Map &InputColMap,
+ const std::vector<unsigned int> &n_entries_per_row)
+ :
+ row_map (InputRowMap),
+ col_map (InputColMap),
matrix (std::auto_ptr<Epetra_FECrsMatrix>
- (new Epetra_FECrsMatrix(Copy, map,
+ (new Epetra_FECrsMatrix(Copy, row_map, col_map,
(int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
- true))),
+ false))),
last_action (Insert)
{}
+
SparseMatrix::~SparseMatrix ()
{
}
{
unsigned int n_rows = sparsity_pattern.n_rows();
-
+
Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
ExcDimensionMismatch (matrix->NumGlobalRows(),
sparsity_pattern.n_rows()));
-
+
std::vector<double> values(n_max_entries_per_row, 0.);
std::vector<int> row_indices(n_max_entries_per_row);
-
+
for (unsigned int row=0; row<n_rows; ++row)
{
const int row_length = sparsity_pattern.row_length(row);
&values[0], &row_indices[0]);
}
- // In the end, the matrix is to
+ // In the end, the matrix needs to
// be compressed in order to be
// really ready. However, that is
// a collective operation, so it
// types are all serial.
}
-
-
+
+
void
SparseMatrix::reinit (const CompressedSparsityPattern &sparsity_pattern)
{
reinit (sparsity_pattern, n_max_entries_per_row);
}
-
-
-
+
+
+
void
SparseMatrix::reinit (const Epetra_Map &input_map,
const CompressedSparsityPattern &sparsity_pattern)
unsigned int n_rows = sparsity_pattern.n_rows();
matrix.reset();
- map = input_map;
+ row_map = input_map;
+
Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(),
ExcDimensionMismatch (input_map.NumGlobalElements(),
sparsity_pattern.n_rows()));
+ Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(),
+ ExcDimensionMismatch (input_map.NumGlobalElements(),
+ 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);
+
+ matrix = std::auto_ptr<Epetra_FECrsMatrix> (new Epetra_FECrsMatrix (
+ Copy, row_map, &n_entries_per_row[0], true));
+
+ const unsigned int n_max_entries_per_row = *std::max_element (
+ &n_entries_per_row[0], &n_entries_per_row[n_rows-1]);
+
+ reinit (sparsity_pattern, n_max_entries_per_row);
+ }
+
+
+
+ void
+ SparseMatrix::reinit (const Epetra_Map &input_row_map,
+ const Epetra_Map &input_col_map,
+ const CompressedSparsityPattern &sparsity_pattern)
+ {
+
+ unsigned int n_rows = sparsity_pattern.n_rows();
+ matrix.reset();
+ row_map = input_row_map;
+ col_map = input_col_map;
+
+ 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()));
std::vector<int> n_entries_per_row(n_rows);
n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
matrix = std::auto_ptr<Epetra_FECrsMatrix>
- (new Epetra_FECrsMatrix(Copy, map, &n_entries_per_row[0], true));
+ (new Epetra_FECrsMatrix(Copy, row_map, col_map,
+ &n_entries_per_row[0], true));
const unsigned int n_max_entries_per_row = *std::max_element (
&n_entries_per_row[0], &n_entries_per_row[n_rows-1]);
// generate an empty matrix.
matrix.reset();
#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
- map = Epetra_Map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)),
+ row_map = Epetra_Map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)),
#else
- map = Epetra_Map (0,0,Epetra_SerialComm()),
+ row_map = Epetra_Map (0,0,Epetra_SerialComm()),
#endif
matrix = std::auto_ptr<Epetra_FECrsMatrix>
- (new Epetra_FECrsMatrix(Copy, map, 0));
+ (new Epetra_FECrsMatrix(Copy, row_map, 0));
}
{
// flush buffers
int ierr;
- ierr = matrix->GlobalAssemble (true);
+ if (row_map.SameAs(col_map))
+ ierr = matrix->GlobalAssemble ();
+ else
+ ierr = matrix->GlobalAssemble (row_map, col_map);
+
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
ierr = matrix->OptimizeStorage ();
if (last_action == Add)
{
int ierr;
- ierr = matrix->GlobalAssemble(false);
+ if (row_map.SameAs(col_map))
+ ierr = matrix->GlobalAssemble (false);
+ else
+ ierr = matrix->GlobalAssemble(row_map, col_map, false);
+
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
last_action = Insert;
}
const_cast<double*>(&value),
&trilinos_j);
+ AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j));
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
}
if (last_action == Insert)
{
int ierr;
- ierr = matrix->GlobalAssemble(false);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ if (row_map.SameAs(col_map))
+ ierr = matrix->GlobalAssemble (false);
+ else
+ ierr = matrix->GlobalAssemble(row_map, col_map, false);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
last_action = Add;
}
const_cast<double*>(&value),
&trilinos_j);
+ AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j));
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
}
// continue.
if ((trilinos_i == -1 ) || (trilinos_j == -1))
{
- Assert (false, ExcAccessToNonlocalElement(i, j, local_range().first,
+ Assert (false, ExcAccessToNonLocalElement(i, j, local_range().first,
local_range().second));
}
else