* holds the sparsity_pattern
* structure because each
* processor sets its rows.
+ *
+ * This is a
+ * collective operation that
+ * needs to be called on all
+ * processors in order to avoid a
+ * dead lock.
*/
void reinit (const SparsityPattern &sparsity_pattern);
* necessary that each processor
* holds the sparsity_pattern
* structure because each
- * processor sets its rows.
+ * processor sets its
+ * rows.
+ *
+ * This is a
+ * collective operation that
+ * needs to be called on all
+ * processors in order to avoid a
+ * dead lock.
*/
void reinit (const Epetra_Map &input_map,
const SparsityPattern &sparsity_pattern);
* user-supplied Epetra maps.
* To be used for rectangular
* matrices.
+ *
+ * This is a
+ * collective operation that
+ * needs to be called on all
+ * processors in order to avoid a
+ * dead lock.
*/
void reinit (const Epetra_Map &input_row_map,
const Epetra_Map &input_col_map,
* content in
* <tt>sparse_matrix</tt> to
* the calling matrix.
+ *
+ * This is a
+ * collective operation that
+ * needs to be called on all
+ * processors in order to avoid a
+ * dead lock.
*/
void reinit (const SparseMatrix &sparse_matrix);
* threshold (so zeros in the
* deal.II matrix can be
* filtered away).
+ *
+ * This is a
+ * collective operation that
+ * needs to be called on all
+ * processors in order to avoid a
+ * dead lock.
*/
void reinit (const Epetra_Map &input_map,
const ::dealii::SparseMatrix<double> &dealii_sparse_matrix,
* the rows and the columns of
* the matrix. Chosen for
* rectangular matrices.
+ *
+ * This is a
+ * collective operation that
+ * needs to be called on all
+ * processors in order to avoid a
+ * dead lock.
*/
void reinit (const Epetra_Map &input_row_map,
const Epetra_Map &input_col_map,
* return to a state just like
* after having called the
* default constructor.
+ *
+ * This is a
+ * collective operation that
+ * needs to be called on all
+ * processors in order to avoid a
+ * dead lock.
*/
void clear ();
* Set the element (<i>i,j</i>)
* to @p value.
*
- * This function adds a new entry
- * to the matrix if it didn't
- * exist before, very much in
- * contrast to the SparseMatrix
- * class which throws an error if
- * the entry does not exist. If
+ * Just as the respective call in
+ * deal.II SparseMatrix<Number>
+ * class (but in contrast to the
+ * situation for PETSc based
+ * matrices), this function
+ * throws an exception if an
+ * entry does not exist in the
+ * sparsity pattern. Moreover, if
* <tt>value</tt> is not a finite
* number an exception is thrown.
*/
* Add @p value to the element
* (<i>i,j</i>).
*
- * This function adds a new
- * entry to the matrix if it
- * didn't exist before, very
- * much in contrast to the
- * SparseMatrix class which
- * throws an error if the entry
- * does not exist. If
- * <tt>value</tt> is not a
- * finite number an exception
- * is thrown.
+ * Just as the respective call in
+ * deal.II SparseMatrix<Number>
+ * class (but in contrast to the
+ * situation for PETSc based
+ * matrices), this function
+ * throws an exception if an
+ * entry does not exist in the
+ * sparsity pattern. Moreover, if
+ * <tt>value</tt> is not a finite
+ * number an exception is thrown.
*/
void add (const unsigned int i,
const unsigned int j,
* operation and you should
* always take care where to
* call this function. As in
- * the &p SparseMatrix<Number>
+ * the deal.II sparse matrix
* class, we throw an exception
* if the respective entry
* doesn't exist in the
* when the requested element
* is not saved on the calling
* process.
- *
- * This function is therefore
- * exactly equivalent to the
- * <tt>el()</tt> function.
*/
TrilinosScalar operator () (const unsigned int i,
const unsigned int j) const;
* <tt>n=local_size()</tt>.
*/
std::pair<unsigned int, unsigned int>
- local_range () const;
+ local_range () const;
/**
* Return whether @p index is
unsigned int row_length (const unsigned int row) const;
/**
- * Return the l1-norm of the
- * matrix, that is
- * $|M|_1=max_{all columns
- * j}\sum_{all rows i} |M_ij|$,
- * (max. sum of columns). This
- * is the natural matrix norm
- * that is compatible to the
- * l1-norm for vectors, i.e.
- * $|Mv|_1\leq |M|_1 |v|_1$.
+ * Return the
+ * <i>l</i><sub>1</sub>-norm of
+ * the matrix, that is
+ * $|M|_1=\max_{\mathrm{all
+ * columns} j}\sum_{\mathrm{all
+ * rows} i} |M_{ij}|$, (max. sum
+ * of columns). This is the
+ * natural matrix norm that is
+ * compatible to the l1-norm for
+ * vectors, i.e. $|Mv|_1\leq
+ * |M|_1 |v|_1$.
* (cf. Haemmerlin-Hoffmann:
* Numerische Mathematik)
*/
TrilinosScalar l1_norm () const;
/**
- * Return the linfty-norm of
- * the matrix, that is
- * $|M|_infty=max_{all rows
- * i}\sum_{all columns j}
- * |M_ij|$, (max. sum of rows).
- * This is the natural matrix
- * norm that is compatible to
- * the linfty-norm of vectors,
- * i.e. $|Mv|_infty \leq
- * |M|_infty |v|_infty$.
+ * Return the linfty-norm of the
+ * matrix, that is
+ * $|M|_\infty=\max_{\mathrm{all
+ * rows} i}\sum_{\mathrm{all
+ * columns} j} |M_{ij}|$,
+ * (max. sum of rows). This is
+ * the natural matrix norm that
+ * is compatible to the
+ * linfty-norm of vectors, i.e.
+ * $|Mv|_\infty \leq |M|_\infty
+ * |v|_\infty$.
* (cf. Haemmerlin-Hoffmann:
* Numerische Mathematik)
*/
SparseMatrix & operator /= (const TrilinosScalar factor);
/**
- * Matrix-vector
- * multiplication: let <i>dst =
- * M*src</i> with <i>M</i>
- * being this matrix.
+ * Matrix-vector multiplication:
+ * let <i>dst = M*src</i> with
+ * <i>M</i> being this matrix.
*
* Source and destination must
* not be the same vector.
*
- * Note that both vectors have
- * to be distributed vectors
+ * Note that both vectors have to
+ * be distributed vectors
* generated using the same Map
- * as was used for the matrix
- * in case you work on a
- * distributed memory
- * architecture, using the
- * interface of the
+ * as was used for the matrix in
+ * case you work on a distributed
+ * memory architecture, using the
+ * interface in the
* TrilinosWrappers::Vector
* class.
*/
const Vector &src) const;
/**
- * Matrix-vector
- * multiplication: let <i>dst =
+ * Matrix-vector multiplication:
+ * let <i>dst =
* M<sup>T</sup>*src</i> with
* <i>M</i> being this
- * matrix. This function does
- * the same as vmult() but
- * takes the transposed matrix.
+ * matrix. This function does the
+ * same as vmult() but takes the
+ * transposed matrix.
*
* Source and destination must
* not be the same vector.
*
- * Note that both vectors have
- * to be distributed vectors
+ * Note that both vectors have to
+ * be distributed vectors
* generated using the same Map
- * as was used for the matrix
- * in case you work on a
- * distributed memory
- * architecture, using the
+ * as was used for the matrix in
+ * case you work on a distributed
+ * memory architecture, using the
* interface of the
* TrilinosWrappers::Vector
* class.
* Source and destination must
* not be the same vector.
*
- * Note that both vectors have
- * to be distributed vectors
+ * Note that both vectors have to
+ * be distributed vectors
* generated using the same Map
- * as was used for the matrix
- * in case you work on a
- * distributed memory
- * architecture, using the
+ * as was used for the matrix in
+ * case you work on a distributed
+ * memory architecture, using the
* interface of the
* TrilinosWrappers::Vector
* class.
* Adding Matrix-vector
* multiplication. Add
* <i>M<sup>T</sup>*src</i> to
- * <i>dst</i> with <i>M</i>
- * being this matrix. This
- * function does the same as
- * vmult_add() but takes the
- * transposed matrix.
+ * <i>dst</i> with <i>M</i> being
+ * this matrix. This function
+ * does the same as vmult_add()
+ * but takes the transposed
+ * matrix.
*
* Source and destination must
* not be the same vector.
*
- * Note that both vectors have
- * to be distributed vectors
+ * Note that both vectors have to
+ * be distributed vectors
* generated using the same Map
- * as was used for the matrix
- * in case you work on a
- * distributed memory
- * architecture, using the
+ * as was used for the matrix in
+ * case you work on a distributed
+ * memory architecture, using the
* interface of the
* TrilinosWrappers::Vector
* class.
const Vector &src) const;
/**
- * Return the square of the
- * norm of the vector $v$ with
- * respect to the norm induced
- * by this matrix,
- * i.e. $\left(v,Mv\right)$. This
- * is useful, e.g. in the
- * finite element context,
- * where the $L_2$ norm of a
- * function equals the matrix
- * norm with respect to the
- * mass matrix of the vector
- * representing the nodal
- * values of the finite element
- * function.
+ * Return the square of the norm
+ * of the vector $v$ with respect
+ * to the norm induced by this
+ * matrix, i.e.,
+ * $\left(v,Mv\right)$. This is
+ * useful, e.g. in the finite
+ * element context, where the
+ * $L_2$ norm of a function
+ * equals the matrix norm with
+ * respect to the mass matrix of
+ * the vector representing the
+ * nodal values of the finite
+ * element function.
*
- * Obviously, the matrix needs
- * to be quadratic for this
+ * Obviously, the matrix needs to
+ * be quadratic for this
* operation.
*
* The implementation of this
* as the one in the @p
* SparseMatrix class used in
* deal.II (i.e. the original
- * one, not the Trilinos
- * wrapper class) since
- * Trilinos doesn't support
- * this operation and needs a
- * temporary vector.
+ * one, not the Trilinos wrapper
+ * class) since Trilinos doesn't
+ * support this operation and
+ * needs a temporary vector.
*
* Note that both vectors have
* to be distributed vectors
matrix (std::auto_ptr<Epetra_FECrsMatrix>
(new Epetra_FECrsMatrix(Copy, row_map,
(int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
- true)))
+ true)))
{}
SparseMatrix::SparseMatrix (const Epetra_Map &InputRowMap,
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)))
+ true)))
{}
unsigned int n_rows = sparsity_pattern.n_rows();
// TODO: As of now, assume that the
- // sparsity pattern only sits at the
- // zeroth processor (completely), let
- // this determine the Trilinos
- // sparsity pattern on that
- // processor, and only broadcast the
- // pattern afterwards.
+ // sparsity pattern sits at the all
+ // processors (completely), and let
+ // each processor set its rows. Since
+ // this is wasteful, a better solution
+ // should be found in the future.
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.
- //
+ // Trilinos has 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
// correct values.
matrix = std::auto_ptr<Epetra_FECrsMatrix>
(new Epetra_FECrsMatrix(Copy, row_map,
- &n_entries_per_row[0], false));
+ &n_entries_per_row[0], true));
reinit (sparsity_pattern);
}
// correct values.
matrix = std::auto_ptr<Epetra_FECrsMatrix>
(new Epetra_FECrsMatrix(Copy, row_map,
- &n_entries_per_row[0], false));
+ &n_entries_per_row[0], true));
std::vector<double> values;
std::vector<int> row_indices;
const_cast<double*>(&value),
&trilinos_j);
- if (ierr > 0)
- std::cout << ierr << " " << m() << " " << n() << std::endl;
AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j));
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
}
{
vmult (dst, x);
dst -= b;
- dst *= -1;
+ dst *= -1.;
return dst.l2_norm();
}
bool
SparseMatrix::is_symmetric (const double tolerance)
{
- //bool truth;
- if (tolerance == 0)
- Assert (false, ExcNotImplemented());
+ (void)tolerance;
+
+ Assert (false, ExcNotImplemented());
return false;
}
bool
SparseMatrix::is_hermitian ()
{
- //bool truth;
-
Assert (false, ExcNotImplemented());
return false;
}