* FullMatrix<double> into the sparse
* matrix, according to row_indices
* and col_indices.
+ *
+ * The optional parameter
+ * <tt>elide_zero_values</tt> can be
+ * used to specify whether zero
+ * values should be inserted anyway
+ * or they should be filtered
+ * away. The default value is
+ * <tt>false</tt>, i.e., even zero
+ * values are inserted/replaced.
*/
void set (const std::vector<unsigned int> &row_indices,
const std::vector<unsigned int> &col_indices,
- const FullMatrix<TrilinosScalar> &full_matrix);
+ const FullMatrix<TrilinosScalar> &full_matrix,
+ const bool elide_zero_values = false);
/**
* Set several elements in the
* column indices as given by
* <tt>col_indices</tt> to the
* respective value.
+ *
+ * The optional parameter
+ * <tt>elide_zero_values</tt> can be
+ * used to specify whether zero
+ * values should be inserted anyway
+ * or they should be filtered
+ * away. The default value is
+ * <tt>false</tt>, i.e., even zero
+ * values are inserted/replaced.
*/
- void set (const unsigned int row,
+ void set (const unsigned int row,
const std::vector<unsigned int> &col_indices,
- const std::vector<TrilinosScalar> &values);
+ const std::vector<TrilinosScalar> &values,
+ const bool elide_zero_values = false);
/**
* 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.
+ * given by <tt>values</tt> in the
+ * given row at column locations
+ * specified by col_indices in the
+ * sparse matrix.
+ *
+ * The optional parameter
+ * <tt>elide_zero_values</tt> can be
+ * used to specify whether zero
+ * values should be inserted anyway
+ * or they should be filtered
+ * away. The default value is
+ * <tt>false</tt>, i.e., even zero
+ * values are inserted/replaced.
*/
- void set (const unsigned int n_rows,
- const unsigned int *row_indices,
+ void set (const unsigned int row,
const unsigned int n_cols,
const unsigned int *col_indices,
- const TrilinosScalar *values);
+ const TrilinosScalar *values,
+ const bool elide_zero_values = false);
/**
* Add @p value to the element
* into the sparse matrix, at
* locations according to row_indices
* and col_indices.
+ *
+ * The optional parameter
+ * <tt>elide_zero_values</tt> can be
+ * used to specify whether zero
+ * values should be added anyway or
+ * these should be filtered away and
+ * only non-zero data is added. The
+ * default value is <tt>true</tt>,
+ * i.e., zero values won't be added
+ * into the matrix.
*/
void add (const std::vector<unsigned int> &row_indices,
const std::vector<unsigned int> &col_indices,
- const FullMatrix<TrilinosScalar> &full_matrix);
+ const FullMatrix<TrilinosScalar> &full_matrix,
+ const bool elide_zero_values = true);
/**
* Set several elements in the
* column indices as given by
* <tt>col_indices</tt> to the
* respective value.
+ *
+ * The optional parameter
+ * <tt>elide_zero_values</tt> can be
+ * used to specify whether zero
+ * values should be added anyway or
+ * these should be filtered away and
+ * only non-zero data is added. The
+ * default value is <tt>true</tt>,
+ * i.e., zero values won't be added
+ * into the matrix.
*/
- void add (const unsigned int row,
+ void add (const unsigned int row,
const std::vector<unsigned int> &col_indices,
- const std::vector<TrilinosScalar> &values);
+ const std::vector<TrilinosScalar> &values,
+ const bool elide_zero_values = true);
/**
* 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.
+ * given by <tt>values</tt> in the
+ * given row at column locations
+ * specified by col_indices in the
+ * sparse matrix.
+ *
+ * The optional parameter
+ * <tt>elide_zero_values</tt> can be
+ * used to specify whether zero
+ * values should be added anyway or
+ * these should be filtered away and
+ * only non-zero data is added. The
+ * default value is <tt>true</tt>,
+ * i.e., zero values won't be added
+ * into the matrix.
*/
- void add (const unsigned int n_rows,
- const unsigned int *row_indices,
+ void add (const unsigned int row,
const unsigned int n_cols,
const unsigned int *col_indices,
- const TrilinosScalar *values);
+ const TrilinosScalar *values,
+ const bool elide_zero_values = true);
/**
* Matrix-vector multiplication:
* writing data from local to global
* matrices.
*/
- std::vector<unsigned int> block_col_indices, local_row_length, local_col_indices;
+ std::vector<unsigned int> block_col_indices, local_col_indices, local_row_length;
};
void
BlockSparseMatrix::set (const std::vector<unsigned int> &row_indices,
const std::vector<unsigned int> &col_indices,
- const FullMatrix<TrilinosScalar> &values)
+ const FullMatrix<TrilinosScalar> &values,
+ const bool elide_zero_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));
+ for (unsigned int i=0; i<row_indices.size(); ++i)
+ set (row_indices[i], col_indices.size(), &col_indices[0], &values(0,0),
+ elide_zero_values);
}
void
BlockSparseMatrix::set (const unsigned int row,
const std::vector<unsigned int> &col_indices,
- const std::vector<TrilinosScalar> &values)
+ const std::vector<TrilinosScalar> &values,
+ const bool elide_zero_values)
{
Assert (col_indices.size() == values.size(),
ExcDimensionMismatch(col_indices.size(), values.size()));
- set (1, &row, col_indices.size(), &col_indices[0], &values[0]);
+ set (row, col_indices.size(), &col_indices[0], &values[0],
+ elide_zero_values);
}
+ // This is a very messy function, since
+ // we need to calculate to each position
+ // the location in the global array.
inline
void
- BlockSparseMatrix::set (const unsigned int n_rows,
- const unsigned int *row_indices,
+ BlockSparseMatrix::set (const unsigned int row,
const unsigned int n_cols,
const unsigned int *col_indices,
- const TrilinosScalar *values)
+ const TrilinosScalar *values,
+ const bool elide_zero_values)
{
// Resize scratch arrays
+ local_row_length.resize (this->n_block_cols());
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
// 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.
+ // matrix.
{
unsigned int current_block = 0, row_length = 0;
block_col_indices[0] = 0;
// 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);
+ for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
{
- 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;
+ if (local_row_length[block_col] == 0)
+ continue;
- block(row_index.first, block_col).set
- (1, &row_index.second,
+ block(row_index.first, block_col).set
+ (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]]);
- }
+ &values[block_col_indices[block_col]],
+ elide_zero_values);
}
}
void
BlockSparseMatrix::add (const std::vector<unsigned int> &row_indices,
const std::vector<unsigned int> &col_indices,
- const FullMatrix<TrilinosScalar> &values)
+ const FullMatrix<TrilinosScalar> &values,
+ const bool elide_zero_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));
+ for (unsigned int i=0; i<row_indices.size(); ++i)
+ add(row_indices[i], col_indices.size(), &col_indices[0],
+ &values(0,0), elide_zero_values);
}
void
BlockSparseMatrix::add (const unsigned int row,
const std::vector<unsigned int> &col_indices,
- const std::vector<TrilinosScalar> &values)
+ const std::vector<TrilinosScalar> &values,
+ const bool elide_zero_values)
{
Assert (col_indices.size() == values.size(),
ExcDimensionMismatch(col_indices.size(), values.size()));
- add (1, &row, col_indices.size(), &col_indices[0], &values[0]);
+ add (row, col_indices.size(), &col_indices[0], &values[0],
+ elide_zero_values);
}
inline
void
- BlockSparseMatrix::add (const unsigned int n_rows,
- const unsigned int *row_indices,
+ BlockSparseMatrix::add (const unsigned int row,
const unsigned int n_cols,
const unsigned int *col_indices,
- const TrilinosScalar *values)
+ const TrilinosScalar *values,
+ const bool elide_zero_values)
{
// TODO: Look over this to find out
// whether we can do that more
// 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.
+ // matrix.
{
unsigned int current_block = 0, row_length = 0;
block_col_indices[0] = 0;
// 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);
+ for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
{
- 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;
+ if (local_row_length[block_col] == 0)
+ continue;
- block(row_index.first, block_col).add
- (1, &row_index.second,
+ block(row_index.first, block_col).add
+ (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]]);
- }
+ &values[block_col_indices[block_col]],
+ elide_zero_values);
}
}
* insertion of elements at positions
* which have not been initialized
* will throw an exception.
+ *
+ * The optional parameter
+ * <tt>elide_zero_values</tt> can be
+ * used to specify whether zero
+ * values should be inserted anyway
+ * or they should be filtered
+ * away. The default value is
+ * <tt>false</tt>, i.e., even zero
+ * values are inserted/replaced.
*/
void set (const std::vector<unsigned int> &row_indices,
const std::vector<unsigned int> &col_indices,
- const FullMatrix<TrilinosScalar> &full_matrix);
+ const FullMatrix<TrilinosScalar> &full_matrix,
+ const bool elide_zero_values = false);
/**
* Set several elements in the
* insertion of elements at positions
* which have not been initialized
* will throw an exception.
+ *
+ * The optional parameter
+ * <tt>elide_zero_values</tt> can be
+ * used to specify whether zero
+ * values should be inserted anyway
+ * or they should be filtered
+ * away. The default value is
+ * <tt>false</tt>, i.e., even zero
+ * values are inserted/replaced.
*/
void set (const unsigned int row,
const std::vector<unsigned int> &col_indices,
- const std::vector<TrilinosScalar> &values);
+ const std::vector<TrilinosScalar> &values,
+ const bool elide_zero_values = false);
/**
* 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.
+ * given by <tt>values</tt> in a
+ * given row in columns given by
+ * col_indices into the sparse
+ * matrix.
*
* This function is able to insert
* new elements into the matrix as
* insertion of elements at positions
* which have not been initialized
* will throw an exception.
+ *
+ * The optional parameter
+ * <tt>elide_zero_values</tt> can be
+ * used to specify whether zero
+ * values should be inserted anyway
+ * or they should be filtered
+ * away. The default value is
+ * <tt>false</tt>, i.e., even zero
+ * values are inserted/replaced.
*/
- void set (const unsigned int n_rows,
- const unsigned int *row_indices,
+ void set (const unsigned int row,
const unsigned int n_cols,
const unsigned int *col_indices,
- const TrilinosScalar *values);
+ const TrilinosScalar *values,
+ const bool elide_zero_values = false);
/**
* Add @p value to the element
* throws an exception if an
* entry does not exist in the
* sparsity pattern.
+ *
+ * The optional parameter
+ * <tt>elide_zero_values</tt> can be
+ * used to specify whether zero
+ * values should be added anyway or
+ * these should be filtered away and
+ * only non-zero data is added. The
+ * default value is <tt>true</tt>,
+ * i.e., zero values won't be added
+ * into the matrix.
*/
void add (const std::vector<unsigned int> &row_indices,
const std::vector<unsigned int> &col_indices,
- const FullMatrix<TrilinosScalar> &full_matrix);
+ const FullMatrix<TrilinosScalar> &full_matrix,
+ const bool elide_zero_values = true);
/**
* Set several elements in the
* throws an exception if an
* entry does not exist in the
* sparsity pattern.
+ *
+ * The optional parameter
+ * <tt>elide_zero_values</tt> can be
+ * used to specify whether zero
+ * values should be added anyway or
+ * these should be filtered away and
+ * only non-zero data is added. The
+ * default value is <tt>true</tt>,
+ * i.e., zero values won't be added
+ * into the matrix.
*/
void add (const unsigned int row,
const std::vector<unsigned int> &col_indices,
- const std::vector<TrilinosScalar> &values);
+ const std::vector<TrilinosScalar> &values,
+ const bool elide_zero_values = true);
/**
- * 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.
+ * Add an array of values given by
+ * <tt>values</tt> in the given
+ * global matrix row at columns
+ * specified by col_indices in the
+ * sparse matrix.
*
* 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
+ * 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.
+ *
+ * The optional parameter
+ * <tt>elide_zero_values</tt> can be
+ * used to specify whether zero
+ * values should be added anyway or
+ * these should be filtered away and
+ * only non-zero data is added. The
+ * default value is <tt>true</tt>,
+ * i.e., zero values won't be added
+ * into the matrix.
*/
- void add (const unsigned int n_rows,
- const unsigned int *row_indices,
+ void add (const unsigned int row,
const unsigned int n_cols,
const unsigned int *col_indices,
- const TrilinosScalar *values);
+ const TrilinosScalar *values,
+ const bool elide_zero_values = false);
/**
* Multiply the entire matrix
*/
mutable VectorBase temp_vector;
+ /**
+ * An internal array of integer
+ * values that is used to store the
+ * column indices when
+ * adding/inserting local data into
+ * the (large) sparse matrix.
+ */
+ std::vector<unsigned int> column_indices;
+
+ /**
+ * An internal array of double values
+ * that is used to store the column
+ * indices when adding/inserting
+ * local data into the (large) sparse
+ * matrix.
+ */
+ std::vector<TrilinosScalar> column_values;
+
public:
/**
* A sparse matrix object in
ExcMessage("The given value is not finite but either "
"infinite or Not A Number (NaN)"));
- set (1, &i, 1, &j, &value);
+ set (i, 1, &j, &value, false);
}
void
SparseMatrix::set (const std::vector<unsigned int> &row_indices,
const std::vector<unsigned int> &col_indices,
- const FullMatrix<TrilinosScalar> &values)
+ const FullMatrix<TrilinosScalar> &values,
+ const bool elide_zero_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));
+ for (unsigned int i=0; i<row_indices.size(); ++i)
+ set (row_indices[i], col_indices.size(), &col_indices[0], &values(0,0),
+ elide_zero_values);
}
void
SparseMatrix::set (const unsigned int row,
const std::vector<unsigned int> &col_indices,
- const std::vector<TrilinosScalar> &values)
+ const std::vector<TrilinosScalar> &values,
+ const bool elide_zero_values)
{
Assert (col_indices.size() == values.size(),
ExcDimensionMismatch(col_indices.size(), values.size()));
- set (1, &row, col_indices.size(), &col_indices[0], &values[0]);
+ set (row, col_indices.size(), &col_indices[0], &values[0],
+ elide_zero_values);
}
inline
void
- SparseMatrix::set (const unsigned int n_rows,
- const unsigned int *row_indices,
+ SparseMatrix::set (const unsigned int row,
const unsigned int n_cols,
const unsigned int *col_indices,
- const TrilinosScalar *values)
+ const TrilinosScalar *values,
+ const bool elide_zero_values)
{
int ierr;
if (last_action == Add)
if (last_action != Insert)
last_action = Insert;
- // Now go through all rows that are
- // present in the input data.
- for (unsigned int i=0; i<n_rows; ++i)
+ int * col_index_ptr;
+ TrilinosScalar const* col_value_ptr;
+
+ // If we don't elide zeros, the pointers
+ // are already available...
+ if (elide_zero_values == false)
+ {
+ col_index_ptr = (int*)col_indices;
+ col_value_ptr = values;
+ }
+ else
{
- const int row = row_indices[i];
+ // Otherwise, extract nonzero values in
+ // each row and pass on to the other
+ // function.
+ column_indices.resize(n_cols);
+ column_values.resize(n_cols);
+
+ unsigned int n_columns = 0;
+ for (unsigned int j=0; j<n_cols; ++j)
+ if (values[j] != 0)
+ {
+ column_indices[n_columns] = col_indices[j];
+ column_values[n_columns] = values[j];
+ n_columns++;
+ }
+ Assert(n_columns <= n_cols, ExcInternalError());
+
+ col_index_ptr = (int*)&column_indices[0];
+ col_value_ptr = &column_values[0];
+ }
+
// If the calling matrix owns the row to
// which we want to insert values, we
// add the possibility to insert new values,
// and in the second we just replace
// data.
- if (row_map.MyGID(row) == true)
+ if (row_map.MyGID(row) == true)
+ {
+ if (matrix->Filled() == false)
{
- if (matrix->Filled() == false)
- {
- ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues(
- row, (int)n_cols,
- const_cast<TrilinosScalar*>(&values[i*n_cols]),
- (int*)&col_indices[0]);
+ ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues(row, n_cols,
+ const_cast<double*>(col_value_ptr),
+ col_index_ptr);
// When inserting elements, we do
// not want to create exceptions in
// the case when inserting non-local
// data (since that's what we want
// to do right now).
- if (ierr > 0)
- ierr = 0;
- }
- else
- ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues(
- row, (int)n_cols,
- const_cast<TrilinosScalar*>(&values[i*n_cols]),
- (int*)&col_indices[0]);
+ if (ierr > 0)
+ ierr = 0;
}
else
- {
+ ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues(row, n_cols,
+ const_cast<double*>(col_value_ptr),
+ col_index_ptr);
+ }
+ else
+ {
// When we're at off-processor data, we
// have to stick with the standard
// Insert/ReplaceGlobalValues
// call the function we already use,
// which is very unefficient if writing
// one element at a time).
- compressed = false;
+ compressed = false;
- const TrilinosScalar* value_ptr = &values[i*n_cols];
-
- if (matrix->Filled() == false)
- {
- ierr = matrix->InsertGlobalValues (1, &row,
- (int)n_cols, (int*)&col_indices[0],
- &value_ptr,
- Epetra_FECrsMatrix::ROW_MAJOR);
- if (ierr > 0)
- ierr = 0;
- }
- else
- ierr = matrix->ReplaceGlobalValues (1, &row,
- (int)n_cols, (int*)&col_indices[0],
- &value_ptr,
- Epetra_FECrsMatrix::ROW_MAJOR);
+ if (matrix->Filled() == false)
+ {
+ ierr = matrix->InsertGlobalValues (1, (int*)&row,
+ (int)n_cols, col_index_ptr,
+ &col_value_ptr,
+ Epetra_FECrsMatrix::ROW_MAJOR);
+ if (ierr > 0)
+ ierr = 0;
}
-
- Assert (ierr <= 0, ExcAccessToNonPresentElement(row_indices[i],
- col_indices[0]));
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ else
+ ierr = matrix->ReplaceGlobalValues (1, (int*)&row,
+ (int)n_cols, col_index_ptr,
+ &col_value_ptr,
+ Epetra_FECrsMatrix::ROW_MAJOR);
}
+
+ Assert (ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0]));
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
}
if (value == 0)
{
- // we have to do above actions in any case
+ // we have to do checkings on Insert/Add
+ // in any case
// to be consistent with the MPI
// communication model (see the comments
// in the documentation of
return;
}
else
- add (1, &i, 1, &j, &value);
+ add (i, 1, &j, &value, false);
}
void
SparseMatrix::add (const std::vector<unsigned int> &row_indices,
const std::vector<unsigned int> &col_indices,
- const FullMatrix<TrilinosScalar> &values)
+ const FullMatrix<TrilinosScalar> &values,
+ const bool elide_zero_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));
+ for (unsigned int i=0; i<row_indices.size(); ++i)
+ add (row_indices[i], col_indices.size(), &col_indices[0], &values(0,0),
+ elide_zero_values);
}
void
SparseMatrix::add (const unsigned int row,
const std::vector<unsigned int> &col_indices,
- const std::vector<TrilinosScalar> &values)
+ const std::vector<TrilinosScalar> &values,
+ const bool elide_zero_values)
{
Assert (col_indices.size() == values.size(),
ExcDimensionMismatch(col_indices.size(), values.size()));
- add (1, &row, col_indices.size(), &col_indices[0], &values[0]);
+ add (row, col_indices.size(), &col_indices[0], &values[0],
+ elide_zero_values);
}
inline
void
- SparseMatrix::add (const unsigned int n_rows,
- const unsigned int *row_indices,
+ SparseMatrix::add (const unsigned int row,
const unsigned int n_cols,
const unsigned int *col_indices,
- const TrilinosScalar *values)
+ const TrilinosScalar *values,
+ const bool elide_zero_values)
{
int ierr;
if (last_action == Insert)
if (last_action != Add)
last_action = Add;
- // Now go through all rows that are
- // present in the input data.
- for (unsigned int i=0; i<n_rows; ++i)
+ int * col_index_ptr;
+ TrilinosScalar const* col_value_ptr;
+
+ // If we don't elide zeros, the pointers
+ // are already available...
+ if (elide_zero_values == false)
{
- const int row = row_indices[i];
+ col_index_ptr = (int*)&col_indices[0];
+ col_value_ptr = values;
+ }
+ else
+ {
+ // Otherwise, extract nonzero values in
+ // each row and pass on to the other
+ // function.
+ column_indices.resize(n_cols);
+ column_values.resize(n_cols);
+
+ unsigned int n_columns = 0;
+ for (unsigned int j=0; j<n_cols; ++j)
+ if (values[j] != 0)
+ {
+ column_indices[n_columns] = col_indices[j];
+ column_values[n_columns] = values[j];
+ n_columns++;
+ }
+ Assert(n_columns <= n_cols, ExcInternalError());
+
+ col_index_ptr = (int*)&column_indices[0];
+ col_value_ptr = &column_values[0];
+ }
+
// 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(row_indices[i]) == true)
- {
- ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues(
- row, (int)n_cols,
- const_cast<TrilinosScalar*>(&values[i*n_cols]),
- (int*)&col_indices[0]);
- }
- else
- {
+ if (row_map.MyGID(row) == true)
+ {
+ ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues(row, n_cols,
+ const_cast<double*>(col_value_ptr),
+ col_index_ptr);
+ }
+ 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;
-
- const TrilinosScalar* value_ptr = &values[i*n_cols];
-
- ierr = matrix->SumIntoGlobalValues (1, &row,
- (int)n_cols, (int*)&col_indices[0],
- &value_ptr,
- Epetra_FECrsMatrix::ROW_MAJOR);
- }
+ compressed = false;
- Assert (ierr <= 0, ExcAccessToNonPresentElement(row_indices[i],
- col_indices[0]));
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ ierr = matrix->SumIntoGlobalValues (1, (int*)&row, (int)n_cols,
+ col_index_ptr,
+ &col_value_ptr,
+ Epetra_FECrsMatrix::ROW_MAJOR);
}
+
+ Assert (ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0]));
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
}