#include <lac/block_sparse_matrix.h>
#include <lac/trilinos_sparse_matrix.h>
#include <lac/trilinos_block_vector.h>
+#include <lac/full_matrix.h>
#include <lac/exceptions.h>
#include <cmath>
*/
unsigned int n_nonzero_elements () const;
+ /**
+ * Set the element (<i>i,j</i>)
+ * to @p value.
+ */
+ void set (const unsigned int i,
+ const unsigned int j,
+ const TrilinosScalar value);
+
+ /**
+ * Set all elements given in a
+ * FullMatrix<double> into the sparse
+ * matrix, according to row_indices
+ * and col_indices.
+ */
+ void set (const std::vector<unsigned int> &row_indices,
+ const std::vector<unsigned int> &col_indices,
+ const FullMatrix<TrilinosScalar> &full_matrix);
+
+ /**
+ * Set several elements in the
+ * specified row of the matrix with
+ * column indices as given by
+ * <tt>col_indices</tt> to the
+ * respective value.
+ */
+ void set (const unsigned int row,
+ const std::vector<unsigned int> &col_indices,
+ const std::vector<TrilinosScalar> &values);
+
+ /**
+ * Set several elements to values
+ * given by <tt>values</tt> at
+ * locations specified by row_indices
+ * and col_indices in the sparse
+ * matrix. The array <tt>values</tt>
+ * is assumed to store data in C
+ * style, i.e., row-wise.
+ */
+ void set (const unsigned int n_rows,
+ const unsigned int *row_indices,
+ const unsigned int n_cols,
+ const unsigned int *col_indices,
+ const TrilinosScalar *values);
+
+ /**
+ * Add @p value to the element
+ * (<i>i,j</i>).
+ */
+ void add (const unsigned int i,
+ const unsigned int j,
+ const TrilinosScalar value);
+
+ /**
+ * Add all elements in a FullMatrix
+ * <tt>full_matrix</tt> (cell matrix)
+ * into the sparse matrix, at
+ * locations according to row_indices
+ * and col_indices.
+ */
+ void add (const std::vector<unsigned int> &row_indices,
+ const std::vector<unsigned int> &col_indices,
+ const FullMatrix<TrilinosScalar> &full_matrix);
+
+ /**
+ * Set several elements in the
+ * specified row of the matrix with
+ * column indices as given by
+ * <tt>col_indices</tt> to the
+ * respective value.
+ */
+ void add (const unsigned int row,
+ const std::vector<unsigned int> &col_indices,
+ const std::vector<TrilinosScalar> &values);
+
+ /**
+ * Set several elements to values
+ * given by <tt>values</tt> at
+ * locations specified by row_indices
+ * and col_indices in the sparse
+ * matrix. The array <tt>values</tt>
+ * is assumed to store data in C
+ * style, i.e., row-wise.
+ */
+ void add (const unsigned int n_rows,
+ const unsigned int *row_indices,
+ const unsigned int n_cols,
+ const unsigned int *col_indices,
+ const TrilinosScalar *values);
+
/**
* Matrix-vector multiplication:
* let $dst = M*src$ with $M$
<< "The blocks [" << arg1 << ',' << arg2 << "] and ["
<< arg3 << ',' << arg4 << "] have differing column numbers.");
///@}
+
+ private:
+ /**
+ * A few temporary vectors for
+ * writing data from local to global
+ * matrices.
+ */
+ std::vector<unsigned int> block_col_indices, local_row_length, local_col_indices;
};
#include <base/config.h>
#include <base/subscriptor.h>
#include <lac/sparse_matrix.h>
+#include <lac/full_matrix.h>
#include <lac/exceptions.h>
#include <lac/trilinos_vector_base.h>
const unsigned int j,
const TrilinosScalar value);
+ /**
+ * Set all elements given in a
+ * FullMatrix<double> into the sparse
+ * matrix, according to row_indices
+ * and col_indices.
+ *
+ * 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.
+ */
+ void set (const std::vector<unsigned int> &row_indices,
+ const std::vector<unsigned int> &col_indices,
+ const FullMatrix<TrilinosScalar> &full_matrix);
+
+ /**
+ * Set several elements in the
+ * specified row of the matrix with
+ * column indices as given by
+ * <tt>col_indices</tt> to the
+ * respective value.
+ *
+ * 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.
+ */
+ void set (const unsigned int row,
+ const std::vector<unsigned int> &col_indices,
+ const std::vector<TrilinosScalar> &values);
+
+ /**
+ * Set several elements to values
+ * given by <tt>values</tt> at
+ * locations specified by row_indices
+ * and col_indices in the sparse
+ * matrix. The array <tt>values</tt>
+ * is assumed to store data in C
+ * style, i.e., row-wise.
+ *
+ * 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.
+ */
+ void set (const unsigned int n_rows,
+ const unsigned int *row_indices,
+ const unsigned int n_cols,
+ const unsigned int *col_indices,
+ const TrilinosScalar *values);
+
/**
* Add @p value to the element
* (<i>i,j</i>).
void add (const unsigned int i,
const unsigned int j,
const TrilinosScalar value);
-
+
+ /**
+ * Add all elements in a FullMatrix
+ * <tt>full_matrix</tt> (cell matrix)
+ * into the sparse matrix, at
+ * locations according to row_indices
+ * and col_indices.
+ *
+ * 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.
+ */
+ void add (const std::vector<unsigned int> &row_indices,
+ const std::vector<unsigned int> &col_indices,
+ const FullMatrix<TrilinosScalar> &full_matrix);
+
+ /**
+ * Set several elements in the
+ * specified row of the matrix with
+ * column indices as given by
+ * <tt>col_indices</tt> to the
+ * respective 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.
+ */
+ void add (const unsigned int row,
+ const std::vector<unsigned int> &col_indices,
+ const std::vector<TrilinosScalar> &values);
+
+ /**
+ * Set several elements to values
+ * given by <tt>values</tt> at
+ * locations specified by row_indices
+ * and col_indices in the sparse
+ * matrix. The array <tt>values</tt>
+ * is assumed to store data in C
+ * style, i.e., row-wise.
+ *
+ * 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.
+ */
+ void add (const unsigned int n_rows,
+ const unsigned int *row_indices,
+ const unsigned int n_cols,
+ const unsigned int *col_indices,
+ const TrilinosScalar *values);
+
/**
* Multiply the entire matrix
* by a fixed factor.
+ void
+ BlockSparseMatrix::set (const unsigned int i,
+ const unsigned int j,
+ const TrilinosScalar value)
+ {
+ BaseClass::set (i, j, value);
+ }
+
+
+
+ void
+ BlockSparseMatrix::set (const std::vector<unsigned int> &row_indices,
+ const std::vector<unsigned int> &col_indices,
+ const FullMatrix<TrilinosScalar> &values)
+ {
+ Assert (row_indices.size() == values.m(),
+ ExcDimensionMismatch(row_indices.size(), values.m()));
+ Assert (col_indices.size() == values.n(),
+ ExcDimensionMismatch(col_indices.size(), values.n()));
+
+ set (row_indices.size(), &row_indices[0],
+ col_indices.size(), &col_indices[0], &values(0,0));
+ }
+
+
+
+ void
+ BlockSparseMatrix::set (const unsigned int row,
+ const std::vector<unsigned int> &col_indices,
+ const std::vector<TrilinosScalar> &values)
+ {
+ Assert (col_indices.size() == values.size(),
+ ExcDimensionMismatch(col_indices.size(), values.size()));
+
+ set (1, &row, col_indices.size(), &col_indices[0], &values[0]);
+ }
+
+
+
+ void
+ BlockSparseMatrix::set (const unsigned int n_rows,
+ const unsigned int *row_indices,
+ const unsigned int n_cols,
+ const unsigned int *col_indices,
+ const TrilinosScalar *values)
+ {
+ // Resize scratch arrays
+ block_col_indices.resize (this->n_block_cols());
+ local_row_length.resize (this->n_block_cols());
+ local_col_indices.resize (n_cols);
+
+ // Clear the content in local_row_length
+ for (unsigned int i=0; i<this->n_block_cols(); ++i)
+ local_row_length[i] = 0;
+
+ // Go through the column indices to find
+ // out which portions of the values
+ // should be written into which block
+ // matrix. This can be done before
+ // starting the loop over all the rows,
+ // since we assume a rectangular set of
+ // matrix data.
+ {
+ unsigned int current_block = 0, row_length = 0;
+ block_col_indices[0] = 0;
+ for (unsigned int j=0; j<n_cols; ++j, ++row_length)
+ {
+ const std::pair<unsigned int, unsigned int>
+ col_index = this->column_block_indices.global_to_local(col_indices[j]);
+ local_col_indices[j] = col_index.second;
+ if (col_index.first > current_block)
+ {
+ local_row_length[current_block] = row_length;
+ row_length = 0;
+ while (col_index.first > current_block)
+ current_block++;
+ block_col_indices[current_block] = j;
+ }
+
+ Assert (col_index.first == current_block,
+ ExcInternalError());
+ }
+ local_row_length[current_block] = row_length;
+ Assert (current_block < this->n_block_cols(),
+ ExcInternalError());
+
+#ifdef DEBUG
+ // If in debug mode, do a check whether
+ // the right length has been obtained.
+ unsigned int length = 0;
+ for (unsigned int i=0; i<this->n_block_cols(); ++i)
+ length += local_row_length[i];
+ Assert (length == n_cols,
+ ExcDimensionMismatch(length, n_cols));
+#endif
+ }
+
+ // Now we found out about where the
+ // individual columns should start and
+ // where we should start reading out
+ // data. Now let's write the data into
+ // the individual blocks!
+ for (unsigned int i=0; i<n_rows; ++i)
+ {
+ const std::pair<unsigned int,unsigned int>
+ row_index = this->row_block_indices.global_to_local (row_indices[i]);
+ for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
+ {
+ if (local_row_length[block_col] == 0)
+ continue;
+
+ block(row_index.first, block_col).set
+ (1, &row_index.second,
+ local_row_length[block_col],
+ &local_col_indices[block_col_indices[block_col]],
+ &values[n_cols*i + block_col_indices[block_col]]);
+ }
+ }
+ }
+
+
+
+ void
+ BlockSparseMatrix::add (const unsigned int i,
+ const unsigned int j,
+ const TrilinosScalar value)
+ {
+ // For adding one single element, it is
+ // faster to rely on the BaseClass add
+ // function, than doing all the strange
+ // operations below.
+ BaseClass::add (i, j, value);
+ }
+
+
+
+ void
+ BlockSparseMatrix::add (const std::vector<unsigned int> &row_indices,
+ const std::vector<unsigned int> &col_indices,
+ const FullMatrix<TrilinosScalar> &values)
+ {
+ Assert (row_indices.size() == values.m(),
+ ExcDimensionMismatch(row_indices.size(), values.m()));
+ Assert (col_indices.size() == values.n(),
+ ExcDimensionMismatch(col_indices.size(), values.n()));
+
+ add (row_indices.size(), &row_indices[0],
+ col_indices.size(), &col_indices[0], &values(0,0));
+ }
+
+
+
+ void
+ BlockSparseMatrix::add (const unsigned int row,
+ const std::vector<unsigned int> &col_indices,
+ const std::vector<TrilinosScalar> &values)
+ {
+ Assert (col_indices.size() == values.size(),
+ ExcDimensionMismatch(col_indices.size(), values.size()));
+
+ add (1, &row, col_indices.size(), &col_indices[0], &values[0]);
+ }
+
+
+
+ void
+ BlockSparseMatrix::add (const unsigned int n_rows,
+ const unsigned int *row_indices,
+ const unsigned int n_cols,
+ const unsigned int *col_indices,
+ const TrilinosScalar *values)
+ {
+ // TODO: Look over this to find out
+ // whether we can do that more
+ // efficiently.
+
+ // Resize scratch arrays
+ block_col_indices.resize (this->n_block_cols());
+ local_row_length.resize (this->n_block_cols());
+ local_col_indices.resize (n_cols);
+
+ // Clear the content in local_row_length
+ for (unsigned int i=0; i<this->n_block_cols(); ++i)
+ local_row_length[i] = 0;
+
+ // Go through the column indices to find
+ // out which portions of the values
+ // should be written into which block
+ // matrix. This can be done before
+ // starting the loop over all the rows,
+ // since we assume a rectangular set of
+ // matrix data.
+ {
+ unsigned int current_block = 0, row_length = 0;
+ block_col_indices[0] = 0;
+ for (unsigned int j=0; j<n_cols; ++j, ++row_length)
+ {
+ const std::pair<unsigned int, unsigned int>
+ col_index = this->column_block_indices.global_to_local(col_indices[j]);
+ local_col_indices[j] = col_index.second;
+ if (col_index.first > current_block)
+ {
+ local_row_length[current_block] = row_length;
+ row_length = 0;
+ while (col_index.first > current_block)
+ current_block++;
+ block_col_indices[current_block] = j;
+ }
+
+ Assert (col_index.first == current_block,
+ ExcInternalError());
+ }
+ local_row_length[current_block] = row_length;
+ Assert (current_block < this->n_block_cols(),
+ ExcInternalError());
+
+#ifdef DEBUG
+ // If in debug mode, do a check whether
+ // the right length has been obtained.
+ unsigned int length = 0;
+ for (unsigned int i=0; i<this->n_block_cols(); ++i)
+ length += local_row_length[i];
+ Assert (length == n_cols,
+ ExcDimensionMismatch(length, n_cols));
+#endif
+ }
+
+ // Now we found out about where the
+ // individual columns should start and
+ // where we should start reading out
+ // data. Now let's write the data into
+ // the individual blocks!
+ for (unsigned int i=0; i<n_rows; ++i)
+ {
+ const std::pair<unsigned int,unsigned int>
+ row_index = this->row_block_indices.global_to_local (row_indices[i]);
+ for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
+ {
+ if (local_row_length[block_col] == 0)
+ continue;
+
+ block(row_index.first, block_col).add
+ (1, &row_index.second,
+ local_row_length[block_col],
+ &local_col_indices[block_col_indices[block_col]],
+ &values[n_cols*i + block_col_indices[block_col]]);
+ }
+ }
+ }
+
+
+
TrilinosScalar
BlockSparseMatrix::residual (MPI::BlockVector &dst,
const MPI::BlockVector &x,
ExcMessage("The given value is not finite but either "
"infinite or Not A Number (NaN)"));
- // If the last thing we did was an
- // addition, we need to do all the
- // data exchange.
+ set (1, &i, 1, &j, &value);
+ }
+
+
+
+ void
+ SparseMatrix::set (const std::vector<unsigned int> &row_indices,
+ const std::vector<unsigned int> &col_indices,
+ const FullMatrix<TrilinosScalar> &values)
+ {
+ Assert (row_indices.size() == values.m(),
+ ExcDimensionMismatch(row_indices.size(), values.m()));
+ Assert (col_indices.size() == values.n(),
+ ExcDimensionMismatch(col_indices.size(), values.n()));
+
+ set (row_indices.size(), &row_indices[0],
+ col_indices.size(), &col_indices[0], &values(0,0));
+ }
+
+
+
+ void
+ SparseMatrix::set (const unsigned int row,
+ const std::vector<unsigned int> &col_indices,
+ const std::vector<TrilinosScalar> &values)
+ {
+ Assert (col_indices.size() == values.size(),
+ ExcDimensionMismatch(col_indices.size(), values.size()));
+
+ set (1, &row, col_indices.size(), &col_indices[0], &values[0]);
+ }
+
+
+
+ void
+ SparseMatrix::set (const unsigned int n_rows,
+ const unsigned int *row_indices,
+ const unsigned int n_cols,
+ const unsigned int *col_indices,
+ const TrilinosScalar *values)
+ {
+ int ierr;
if (last_action == Add)
{
- int ierr;
ierr = matrix->GlobalAssemble(col_map, row_map, false);
-
+
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
}
- // If our last action was not an
- // addition, we tell that at this
- // point here.
if (last_action != Insert)
last_action = Insert;
- compressed = false;
-
- int ierr;
-
- // In case this matrix owns the
- // actual row, we can save some
- // computing time by directly
- // accessing the Epetra_CrsMatrix
- // set functions.
- if (row_map.MyGID(i) == true)
+ // Now go through all rows that are
+ // present in the input data.
+ for (unsigned int i=0; i<n_rows; ++i)
{
- // 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)
+ // If the calling matrix owns the row to
+ // which we want to add values, we
+ // can directly call the Epetra_CrsMatrix
+ // input function, which is much faster
+ // than the Epetra_FECrsMatrix function.
+ if (row_map.MyGID(i) == true)
{
- ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues (
- (int)i, 1,
- const_cast<TrilinosScalar*>(&value),
- (int*)&j);
+ if (matrix->Filled() == false)
+ {
+ ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues(
+ (int)row_indices[i], (int)n_cols,
+ const_cast<TrilinosScalar*>(&values[i*n_cols]),
+ (int*)&col_indices[0]);
// When adding up elements, we do
// not want to create exceptions in
// the case when adding elements.
- if (ierr > 0)
- ierr = 0;
+ if (ierr > 0)
+ ierr = 0;
+ }
+ else
+ ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues(
+ (int)row_indices[i], (int)n_cols,
+ const_cast<TrilinosScalar*>(&values[i*n_cols]),
+ (int*)&col_indices[0]);
}
else
{
- ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues (
- (int)i, 1,
- const_cast<TrilinosScalar*>(&value),
- (int*)&j);
+ // When we're at off-processor data, we
+ // have to stick with the standard
+ // SumIntoGlobalValues
+ // function. Nevertheless, the way we
+ // call it is the fastest one (any other
+ // will lead to repeated allocation and
+ // deallocation of memory in order to
+ // call the function we already use,
+ // which is very unefficient if writing
+ // one element at a time).
+ compressed = false;
+
+ const TrilinosScalar* value_ptr = &values[i*n_cols];
+
+ if (matrix->Filled() == false)
+ {
+ ierr = matrix->InsertGlobalValues (1, (int*)&i,
+ (int)n_cols, (int*)&col_indices[0],
+ &value_ptr,
+ Epetra_FECrsMatrix::ROW_MAJOR);
+ if (ierr > 0)
+ ierr = 0;
+ }
+ else
+ ierr = matrix->ReplaceGlobalValues (1, (int*)&i,
+ (int)n_cols, (int*)&col_indices[0],
+ &value_ptr,
+ Epetra_FECrsMatrix::ROW_MAJOR);
}
+
+ Assert (ierr <= 0, ExcAccessToNonPresentElement(row_indices[i],
+ col_indices[0]));
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
}
- else
- {
- int trilinos_i = i;
- int trilinos_j = j;
- const TrilinosScalar* value_ptr = &value;
+ }
- // 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 (1, &trilinos_i,
- 1, &trilinos_j,
- &value_ptr,
- Epetra_FECrsMatrix::ROW_MAJOR);
- // When adding up elements, we do
- // not want to create exceptions in
- // the case when adding elements.
- if (ierr > 0)
- ierr = 0;
- }
- else
+ void
+ SparseMatrix::add (const unsigned int i,
+ const unsigned int j,
+ const TrilinosScalar value)
+ {
+
+ Assert (numbers::is_finite(value),
+ ExcMessage("The given value is not finite but either "
+ "infinite or Not A Number (NaN)"));
+
+ if (value == 0)
+ {
+ // we have to do above actions in any case
+ // to be consistent with the MPI
+ // communication model (see the comments
+ // in the documentation of
+ // TrilinosWrappers::Vector), but we can
+ // save some work if the addend is
+ // zero. However, these actions are done
+ // in case we pass on to the other
+ // function.
+ if (last_action == Insert)
{
- ierr = matrix->ReplaceGlobalValues (1, &trilinos_i,
- 1, &trilinos_j,
- &value_ptr,
- Epetra_FECrsMatrix::ROW_MAJOR);
+ int ierr;
+ ierr = matrix->GlobalAssemble(col_map, row_map, false);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
}
+
+ if (last_action != Add)
+ last_action = Add;
+
+ return;
}
+ else
+ add (1, &i, 1, &j, &value);
+ }
- Assert (ierr <= 0, ExcAccessToNonPresentElement(i,j));
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ void
+ SparseMatrix::add (const std::vector<unsigned int> &row_indices,
+ const std::vector<unsigned int> &col_indices,
+ const FullMatrix<TrilinosScalar> &values)
+ {
+ Assert (row_indices.size() == values.m(),
+ ExcDimensionMismatch(row_indices.size(), values.m()));
+ Assert (col_indices.size() == values.n(),
+ ExcDimensionMismatch(col_indices.size(), values.n()));
+
+ add (row_indices.size(), &row_indices[0],
+ col_indices.size(), &col_indices[0], &values(0,0));
}
void
- SparseMatrix::add (const unsigned int i,
- const unsigned int j,
- const TrilinosScalar value)
+ SparseMatrix::add (const unsigned int row,
+ const std::vector<unsigned int> &col_indices,
+ const std::vector<TrilinosScalar> &values)
{
+ Assert (col_indices.size() == values.size(),
+ ExcDimensionMismatch(col_indices.size(), values.size()));
+
+ add (1, &row, col_indices.size(), &col_indices[0], &values[0]);
+ }
+
- Assert (numbers::is_finite(value),
- ExcMessage("The given value is not finite but either "
- "infinite or Not A Number (NaN)"));
+ void
+ SparseMatrix::add (const unsigned int n_rows,
+ const unsigned int *row_indices,
+ const unsigned int n_cols,
+ const unsigned int *col_indices,
+ const TrilinosScalar *values)
+ {
+ int ierr;
if (last_action == Insert)
{
- int ierr;
ierr = matrix->GlobalAssemble(col_map, row_map, false);
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
+ }
if (last_action != Add)
last_action = Add;
- // we have to do above actions in any
- // case to be consistent with the MPI
- // communication model (see the
- // comments in the documentation of
- // TrilinosWrappers::Vector), but we
- // can save some work if the addend is
- // zero
- if (value == 0)
- return;
-
- int ierr;
+ // Now go through all rows that are
+ // present in the input data.
+ for (unsigned int i=0; i<n_rows; ++i)
+ {
// If the calling matrix owns the row to
- // which we want to add this value, we
+ // which we want to add values, we
// can directly call the Epetra_CrsMatrix
// input function, which is much faster
// than the Epetra_FECrsMatrix function.
- if (row_map.MyGID(i) == true)
- {
- ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues(
- (int)i, 1,
- const_cast<TrilinosScalar*>(&value),
- (int*)&j);
- }
- else
- {
+ if (row_map.MyGID(i) == true)
+ {
+ ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues(
+ (int)row_indices[i], (int)n_cols,
+ const_cast<TrilinosScalar*>(&values[i*n_cols]),
+ (int*)&col_indices[0]);
+ }
+ else
+ {
// When we're at off-processor data, we
// have to stick with the standard
// SumIntoGlobalValues
// call the function we already use,
// which is very unefficient if writing
// one element at a time).
- compressed = false;
+ compressed = false;
- int trilinos_i = i;
- int trilinos_j = j;
+ const TrilinosScalar* value_ptr = &values[i*n_cols];
- const TrilinosScalar* value_ptr = &value;
+ ierr = matrix->SumIntoGlobalValues (1, (int*)&i,
+ (int)n_cols, (int*)&col_indices[0],
+ &value_ptr,
+ Epetra_FECrsMatrix::ROW_MAJOR);
+ }
- ierr = matrix->SumIntoGlobalValues (1, &trilinos_i,
- 1, &trilinos_j,
- &value_ptr,
- Epetra_FECrsMatrix::ROW_MAJOR);
+ Assert (ierr <= 0, ExcAccessToNonPresentElement(row_indices[i],
+ col_indices[0]));
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
}
-
- AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j));
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
}
{
Assert (&src != &dst, ExcSourceEqualsDestination());
+ // Choose to reinit the vector with fast
+ // argument set, which does not overwrite
+ // the content -- this is what we need
+ // since we're going to overwrite that
+ // anyway in the vmult operation.
temp_vector.reinit(dst, true);
vmult (temp_vector, src);