typedef TrilinosScalar value_type;
/**
- * Default constructor.
+ * Default
+ * constructor. Generates an
+ * empty (zero-size) matrix.
*/
SparseMatrix ();
const Epetra_Map &input_col_map,
const SparsityPattern &sparsity_pattern);
- /**
+ /**
* This function copies the
* content in
* <tt>sparse_matrix</tt> to
* the matrix without having to
* read entries (such as the
* locations of non-zero
- * elements) from it -- without
- * this o peration, removing
- * constraints on parallel
- * matrices is a rather
- * complicated procedure.
+ * elements) from it —
+ * without this operation,
+ * removing constraints on
+ * parallel matrices is a
+ * rather complicated
+ * procedure.
*
* The second parameter can be
* used to set the diagonal
* may be an expensive
* operation and you should
* always take care where to
- * call this function. In
- * contrast to the respective
- * function in the @p
- * SparseMatrix class, we don't
- * throw an exception if the
- * respective entry doesn't
- * exist in the sparsity
- * pattern of this class, since
- * Trilinos does not transmit
- * this information. On the
- * other hand, an exception
- * will be thrown when the
- * requested element is not
- * saved on the calling
+ * call this function. As in
+ * the &p SparseMatrix<Number>
+ * class, we throw an exception
+ * if the respective entry
+ * doesn't exist in the
+ * sparsity pattern of this
+ * class, which is requested
+ * from Trilinos. Moreover, an
+ * exception will be thrown
+ * when the requested element
+ * is not saved on the calling
* process.
*
* This function is therefore
<< "An error with error number " << arg1
<< " occured while calling a Trilinos function");
+ /**
+ * Exception
+ */
+ DeclException2 (ExcInvalidIndex,
+ int, int,
+ << "The entry with index <" << arg1 << ',' << arg2
+ << "> does not exist.");
+
/**
* Exception
*/
* elements. The actual type,
* a sparse matrix, is set in
* the constructor.
+ *
+ * TODO: This object should
+ * finally become private.
*/
std::auto_ptr<Epetra_FECrsMatrix> matrix;
}
+
inline
unsigned int
const_iterator::Accessor::column() const
}
+
inline
unsigned int
const_iterator::Accessor::index() const
}
+
inline
TrilinosScalar
const_iterator::Accessor::value() const
}
+
inline
const_iterator::
const_iterator(const SparseMatrix *matrix,
}
+
inline
const_iterator
const_iterator::operator++ (int)
}
+
inline
const const_iterator::Accessor &
const_iterator::operator* () const
}
+
inline
const const_iterator::Accessor *
const_iterator::operator-> () const
}
+
inline
bool
const_iterator::
}
+
inline
bool
const_iterator::
}
+
inline
bool
const_iterator::
}
}
-
-
- inline
- TrilinosScalar
- SparseMatrix::operator() (const unsigned int i,
- const unsigned int j) const
- {
- return el(i,j);
- }
+
-
inline
SparseMatrix::const_iterator
}
+
inline
SparseMatrix::const_iterator
SparseMatrix::end() const
}
+
inline
SparseMatrix::const_iterator
SparseMatrix::begin(const unsigned int r) const
}
+
inline
SparseMatrix::const_iterator
SparseMatrix::end(const unsigned int r) const
// flush buffers
int ierr;
if (row_map.SameAs(col_map))
- ierr = matrix->GlobalAssemble ();
+ ierr = matrix->GlobalAssemble (true);
else
- ierr = matrix->GlobalAssemble (col_map, row_map);
+ ierr = matrix->GlobalAssemble (col_map, row_map, true);
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ TrilinosScalar
+ SparseMatrix::operator() (const unsigned int i,
+ const unsigned int j) const
+ {
+ // Extract local indices in
+ // the matrix.
+ int trilinos_i = matrix->LRID(i), trilinos_j = matrix->LRID(j);
+ TrilinosScalar value = 0.;
+
+ // If the data is not on the
+ // present processor, we throw
+ // an exception. This is on of
+ // the two tiny differences to
+ // the el(i,j) call, which does
+ // not throw any assertions.
+ if ((trilinos_i == -1 ) || (trilinos_j == -1))
+ {
+ Assert (false, ExcAccessToNonLocalElement(i, j, local_range().first,
+ local_range().second));
+ }
+ else
+ {
+ // Check whether the matrix
+ // already is transformed to
+ // local indices.
+ if (matrix->Filled() == false)
+ matrix->FillComplete(true);
+
+ // Prepare pointers for extraction
+ // of a view of the row.
+ int nnz_present = matrix->NumMyEntries(trilinos_i);
+ int nnz_extracted;
+ int *col_indices;
+ TrilinosScalar *values;
+
+ // Generate the view and make
+ // sure that we have not generated
+ // an error.
+ int ierr = matrix->ExtractMyRowView(trilinos_i, nnz_extracted,
+ values, col_indices);
+ Assert (ierr==0, ExcTrilinosError(ierr));
+
+ Assert (nnz_present == nnz_extracted,
+ ExcDimensionMismatch(nnz_present, nnz_extracted));
+
+ // Search the index where we
+ // look for the value, and then
+ // finally get it.
+
+ int* el_find = std::find(&col_indices[0],&col_indices[0] + nnz_present,
+ trilinos_j);
+
+ int el_index = (int)(el_find - col_indices);
+
+ // This is actually the only
+ // difference to the el(i,j)
+ // function, which means that
+ // we throw an exception in
+ // this case instead of just
+ // returning zero for an
+ // element that is not present
+ // in the sparsity pattern.
+ if (!el_find)
+ {
+ Assert (false, ExcInvalidIndex (i,j));
+ }
+ else
+ value = values[el_index];
+
+ }
+
+ return value;
+ }
+
+
+
TrilinosScalar
SparseMatrix::el (const unsigned int i,
const unsigned int j) const
// continue. Just print out
// zero.
- // TODO: Is this reasonable?
+ // TODO: Is this reasonable? Or
+ // should we retain the assert
+ // call?
if ((trilinos_i == -1 ) || (trilinos_j == -1))
{
return 0.;
TrilinosScalar
SparseMatrix::l1_norm () const
{
- if (!matrix->Filled())
+ if (matrix->Filled() == false)
matrix->FillComplete();
TrilinosScalar result = matrix->NormOne();
TrilinosScalar
SparseMatrix::linfty_norm () const
{
- if (!matrix->Filled())
+ if (matrix->Filled() == false)
matrix->FillComplete();
TrilinosScalar result = matrix->NormInf();
TrilinosScalar
SparseMatrix::frobenius_norm () const
{
- if (!matrix->Filled())
+ if (matrix->Filled() == false)
matrix->FillComplete();
TrilinosScalar result = matrix->NormFrobenius();