* Set the element (<i>i,j</i>)
* to @p value.
*
- * 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
+ * This function is able to insert
+ * new elements into the matrix as
+ * long as compress() has not been
+ * called, so the sparsity pattern
+ * will be extended. When compress()
+ * is called for the first time, then
+ * this is no longer possible and an
+ * insertion of elements at positions
+ * which have not been initialized
+ * will throw an exception. Moreover,
+ * if <tt>value</tt> is not a finite
* number an exception is thrown.
*/
void set (const unsigned int i,
int trilinos_i = i;
int trilinos_j = j;
- const int ierr = matrix->ReplaceGlobalValues (trilinos_i, 1,
- const_cast<double*>(&value),
- &trilinos_j);
+ int ierr;
+
+ // If the matrix is not yet filled,
+ // we can insert new entries into
+ // the matrix using this
+ // command. Otherwise, we're just
+ // able to replace existing entries.
+ if (matrix->Filled() == false)
+ ierr = matrix->InsertGlobalValues (trilinos_i, 1,
+ const_cast<double*>(&value),
+ &trilinos_j);
+ else
+ ierr = matrix->ReplaceGlobalValues (trilinos_i, 1,
+ const_cast<double*>(&value),
+ &trilinos_j);
AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j));
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
// already is transformed to
// local indices.
if (matrix->Filled() == false)
- matrix->FillComplete(col_map, row_map, true);
+ matrix->GlobalAssemble(col_map, row_map, true);
// Prepare pointers for extraction
// of a view of the row.
// look for the value, and then
// finally get it.
- int* el_find = std::find(&col_indices[0],&col_indices[0] + nnz_present,
+ int* el_find = std::find(col_indices, col_indices + nnz_present,
trilinos_j);
-
- int el_index = (int)(el_find - col_indices);
+
+ int local_col_index = (int)(el_find - col_indices);
// This is actually the only
// difference to the el(i,j)
// returning zero for an
// element that is not present
// in the sparsity pattern.
- if (!el_find)
+ if (local_col_index == nnz_present)
{
Assert (false, ExcInvalidIndex (i,j));
}
else
- value = values[el_index];
-
+ value = values[local_col_index];
}
return value;
// already is transformed to
// local indices.
if (!matrix->Filled())
- matrix->FillComplete(col_map, row_map, true);
+ matrix->GlobalAssemble(col_map, row_map, true);
// Prepare pointers for extraction
// of a view of the row.
// 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);
+ int* el_find = std::find(col_indices, col_indices + nnz_present,
+ trilinos_j);
+
+ int local_col_index = (int)(el_find - col_indices);
- if (!el_find)
+
+ // This is actually the only
+ // difference to the () function
+ // querying (i,j), where we throw an
+ // exception instead of just
+ // returning zero for an element
+ // that is not present in the
+ // sparsity pattern.
+ if (local_col_index == nnz_present)
value = 0;
else
- {
- value = values[el_index];
- }
+ value = values[local_col_index];
}
return value;