// then eliminate these rows and set
// their diagonal entry to what we have
- // determined above. note that for trilinos
- // matrices interleaving read with write
- // operations is very expensive. thus, we
- // here always replace the diagonal
- // element, rather than first checking
- // whether it is nonzero and in that case
- // preserving it. this is different from
- // the case of deal.II sparse matrices
- // treated in the other functions.
-//TODO: clear_row is not currently implemented for Trilinos
- Assert (false, ExcInternalError());
-// matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
+ // determined above. if the value already
+ // is nonzero, it will be preserved,
+ // in accordance with the basic
+ // matrix classes in deal.II.
+ matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
// the next thing is to set right hand
// side to the wanted value. there's one
*/
void clear ();
+ /**
+ * Trilinos matrices store their own
+ * sparsity patterns. So, in analogy to
+ * our own SparsityPattern class,
+ * this function compresses the
+ * sparsity pattern and allows the
+ * resulting matrix to be used in all
+ * other operations where before only
+ * assembly functions were
+ * allowed. This function must
+ * therefore be called once you have
+ * assembled the matrix.
+ */
+ void compress ();
+
/**
* This operator assigns a scalar to a
* matrix. Since this does usually not
* pattern, though (but retains
* the allocated memory in case
* new entries are again added
- * later).
+ * later). Note that this
+ * is a global operation, so this
+ * needs to be done on all
+ * MPI processes.
*
* This operation is used in
* eliminating constraints (e.g. due to
* the matrix without having to read
* entries (such as the locations of
* non-zero elements) from it --
- * without this operation, removing
+ * without this o peration, removing
* constraints on parallel matrices is
* a rather complicated procedure.
*
* to a value different from zero. The
* default is to set it to zero.
*/
- // void clear_row (const unsigned int row,
- // const TrilinosScalar new_diag_value = 0);
+ void clear_row (const unsigned int row,
+ const TrilinosScalar new_diag_value = 0);
/**
* Same as clear_row(), except that it
* the diagonal entries, you have to
* set them by hand.
*/
- // void clear_rows (const std::vector<unsigned int> &rows,
- // const TrilinosScalar new_diag_value = 0);
-
- /**
- * Trilinos matrices store their own
- * sparsity patterns. So, in analogy to
- * our own SparsityPattern class,
- * this function compresses the
- * sparsity pattern and allows the
- * resulting matrix to be used in all
- * other operations where before only
- * assembly functions were
- * allowed. This function must
- * therefore be called once you have
- * assembled the matrix.
- */
- void compress ();
-
+ void clear_rows (const std::vector<unsigned int> &rows,
+ const TrilinosScalar new_diag_value = 0);
+
/**
* Return the value of the entry
* (<i>i,j</i>). This may be an
*/
void write_ascii ();
+ /**
+ * Print the matrix to the given
+ * stream, using the format
+ * <tt>(line,col) value</tt>,
+ * i.e. one nonzero entry of the
+ * matrix per line.
+ */
+ void print (std::ostream &out) const;
+
// TODO: Write an overloading
// of the operator << for output.
// Since the underlying Trilinos
void ratio (const Vector &a,
const Vector &b);
+ /**
+ * Output of vector in user-defined
+ * format in analogy to the
+ * dealii::Vector<number> class.
+ */
+ void print (const char* format = 0) const;
+
/**
* Print to a
* stream. @p precision denotes
*
* @ingroup TrilinosWrappers
* @relates TrilinosWrappers::Vector
- * @author Wolfgang Bangerth, 2004
+ * @author Martin Kronbichler, Wolfgang Bangerth, 2008
*/
inline
void swap (Vector &u, Vector &v)
Epetra_Vector LHS (View, multigrid_operator->OperatorDomainMap(),
dst.begin());
- Epetra_Vector RHS (View, multigrid_operator->OperatorDomainMap(),
+ Epetra_Vector RHS (View, multigrid_operator->OperatorRangeMap(),
const_cast<double*>(src.begin()));
const int res = multigrid_operator->ApplyInverse (RHS, LHS);
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.
- //
- // this bug is filed in the Sandia
- // bugzilla under #4123 and should be
- // fixed for version 9.0
+ // Trilinos seems to have a bug for
+ // rectangular matrices at this point,
+ // so do not check for consistent
+ // column numbers here.
// Assert (matrix->NumGlobalCols() == (int)sparsity_pattern.n_cols(),
// ExcDimensionMismatch (matrix->NumGlobalCols(),
// sparsity_pattern.n_cols()));
for (dealii::SparseMatrix<double>::const_iterator
p = deal_ii_sparse_matrix.begin(row);
p != deal_ii_sparse_matrix.end(row); ++p)
- if (std::abs(p->value()) > drop_tolerance)
+ if (std::fabs(p->value()) > drop_tolerance)
{
row_indices[index] = p->column();
values[index] = p->value();
// case to be consistent with the MPI
// communication model (see the
// comments in the documentation of
- // TrilinosWrappers::MPI::Vector), but we
+ // TrilinosWrappers::Vector), but we
// can save some work if the addend is
// zero
if (value == 0)
+ void
+ SparseMatrix::clear_row (const unsigned int row,
+ const TrilinosScalar new_diag_value)
+ {
+ Assert (matrix->Filled()==true,
+ ExcMessage("Matrix must be compressed before invoking clear_row."));
+
+ // Only do this on the rows
+ // owned locally on this processor.
+ int local_row = matrix->LRID(row);
+ if (local_row >= 0)
+ {
+ TrilinosScalar *values;
+ int *col_indices;
+ int num_entries;
+ const int ierr = matrix->ExtractMyRowView(local_row, num_entries,
+ values, col_indices);
+
+ Assert (ierr == 0,
+ ExcTrilinosError(ierr));
+
+ int* diag_find = std::find(col_indices,col_indices+num_entries,
+ local_row);
+ int diag_index = (int)(diag_find - col_indices);
+
+ for (int j=0; j<num_entries; ++j)
+ if (diag_index != col_indices[j])
+ values[j] = 0.;
+
+ if (diag_find && std::fabs(values[diag_index]) > 0.)
+ values[diag_index] = new_diag_value;
+ }
+ }
+
+
+
+ void
+ SparseMatrix::clear_rows (const std::vector<unsigned int> &rows,
+ const TrilinosScalar new_diag_value)
+ {
+ for (unsigned int row=0; row<rows.size(); ++row)
+ clear_row(rows[row], new_diag_value);
+ }
+
+
+
TrilinosScalar
SparseMatrix::el (const unsigned int i,
const unsigned int j) const
// 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 can't
// continue.
// Search the index where we
// look for the value, and then
// finally get it.
- int* index = std::find(&col_indices[0],&col_indices[0] + nnz_present,
- trilinos_j);
+
+ int* el_find = std::find(&col_indices[0],&col_indices[0] + nnz_present,
+ trilinos_j);
+
+ int el_index = (int)(el_find - col_indices);
- int position;
- if (!index)
+ if (!el_find)
value = 0;
else
{
- position = (int)(index - &(col_indices[0]));
- value = values[position];
+ value = values[el_index];
}
}
{
Assert (m() == n(), ExcNotQuadratic());
- // this doesn't seem to work any
- // different than any other element
return el(i,i);
}
{
int begin, end;
begin = matrix->RowMap().MinMyGID();
- end = matrix->RowMap().MaxMyGID();
+ end = matrix->RowMap().MaxMyGID()+1;
return std::make_pair (begin, end);
}
return false;
}
+
+
void
SparseMatrix::write_ascii ()
{
Assert (false, ExcNotImplemented());
}
+
+
+ // As of now, no particularly neat
+ // ouput is generated in case of
+ // multiple processors.
+ void SparseMatrix::print (std::ostream &out) const
+ {
+ double * values;
+ int * indices;
+ int num_entries;
+
+ for (int i=0; i<matrix->NumMyRows(); ++i)
+ {
+ matrix->ExtractMyRowView (i, num_entries, values, indices);
+ for (int j=0; j<num_entries; ++j)
+ out << "(" << i << "," << indices[matrix->GRID(j)] << ") "
+ << values[j] << std::endl;
+ }
+
+ AssertThrow (out, ExcIO());
+ }
+
}
DEAL_II_NAMESPACE_CLOSE
{
int begin, end;
begin = vector->Map().MinMyGID();
- end = vector->Map().MaxMyGID();
+ end = vector->Map().MaxMyGID()+1;
return std::make_pair (begin, end);
}
ExcMessage("The given value is not finite but "
"either infinite or Not A Number (NaN)"));
- std::vector<TrilinosScalar> list (size(), s);
-
- int* index = new int[size()];
- for (unsigned int i=0; i<size(); i++)
- index[i]=i;
-
- const int ierr = vector->SumIntoGlobalValues(size(), index, &list[0]);
+ unsigned int n_local = local_size();
+ int ierr;
- delete[] index;
-
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ for (unsigned int i=0; i<n_local; i++)
+ {
+ ierr = vector->SumIntoMyValue(i,0,s);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
}
void
Vector::add (const TrilinosScalar a,
- const Vector &v)
+ const Vector &v)
{
Assert (numbers::is_finite(a),
ExcMessage("The given value is not finite but "
"either infinite or Not A Number (NaN)"));
- const int ierr = vector->Update(a, *(v.vector), 1);
+ const int ierr = vector->Update(a, *(v.vector), 1.);
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
}
void
Vector::add (const TrilinosScalar a,
- const Vector &v,
- const TrilinosScalar b,
- const Vector &w)
+ const Vector &v,
+ const TrilinosScalar b,
+ const Vector &w)
{
Assert (numbers::is_finite(a),
ExcMessage("The given value is not finite but "
"either infinite or Not A Number (NaN)"));
- const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.0);
+ const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
}
void
Vector::sadd (const TrilinosScalar s,
- const Vector &v)
+ const Vector &v)
{
Assert (numbers::is_finite(s),
ExcMessage("The given value is not finite but "
"either infinite or Not A Number (NaN)"));
- const int ierr = vector->Update(1.0, *(v.vector), s);
+ const int ierr = vector->Update(1., *(v.vector), s);
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
}
void
Vector::sadd (const TrilinosScalar s,
- const TrilinosScalar a,
- const Vector &v)
+ const TrilinosScalar a,
+ const Vector &v)
{
Assert (numbers::is_finite(s),
void
Vector::sadd (const TrilinosScalar s,
- const TrilinosScalar a,
- const Vector &v,
- const TrilinosScalar b,
- const Vector &w)
+ const TrilinosScalar a,
+ const Vector &v,
+ const TrilinosScalar b,
+ const Vector &w)
{
Assert (numbers::is_finite(s),
void
Vector::sadd (const TrilinosScalar s,
- const TrilinosScalar a,
- const Vector &v,
- const TrilinosScalar b,
- const Vector &w,
- const TrilinosScalar c,
- const Vector &x)
+ const TrilinosScalar a,
+ const Vector &v,
+ const TrilinosScalar b,
+ const Vector &w,
+ const TrilinosScalar c,
+ const Vector &x)
{
Assert (numbers::is_finite(s),
ExcMessage("The given value is not finite but "
"either infinite or Not A Number (NaN)"));
- // Update member can only input two other vectors so
+ // Update member can only input
+ // two other vectors so
// do it in two steps
const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- const int jerr = vector->Update(c, *(x.vector), 1.0);
+
+ const int jerr = vector->Update(c, *(x.vector), 1.);
AssertThrow (jerr == 0, ExcTrilinosError(jerr));
}
ExcMessage("The given value is not finite but "
"either infinite or Not A Number (NaN)"));
- Assert (size() == v.size(),
- ExcDimensionMismatch (size(), v.size()));
-
- const int ierr = vector->Update(a, *(v.vector), 0.0);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ *vector = *v.vector;
+ map = v.map;
+
+ *this *= a;
}
ExcMessage("The given value is not finite but "
"either infinite or Not A Number (NaN)"));
- Assert (size() == v.size(),
- ExcDimensionMismatch (size(), v.size()));
-
- Assert (size() == w.size(),
- ExcDimensionMismatch (size(), w.size()));
-
- const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 0.0);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ Assert (v.size() == w.size(),
+ ExcDimensionMismatch (v.size(), w.size()));
+ *vector = *v.vector;
+ map = v.map;
+ sadd (a, b, w);
}
Vector::ratio (const Vector &v,
const Vector &w)
{
- Assert (size() == v.size(),
- ExcDimensionMismatch (size(), v.size()));
+ Assert (v.size() == w.size(),
+ ExcDimensionMismatch (v.size(), w.size()));
Assert (size() == w.size(),
ExcDimensionMismatch (size(), w.size()));
}
+
+ // TODO: up to now only local data
+ // printed out! Find a way to neatly
+ // output distributed data...
+ void
+ Vector::print (const char *format) const
+ {
+ Assert (vector->GlobalLength()!=0, ExcEmptyObject());
+
+ for (unsigned int j=0; j<size(); ++j)
+ {
+ double t = (*vector)[0][j];
+
+ if (format != 0)
+ std::printf (format, t);
+ else
+ std::printf (" %5.2f", double(t));
+ }
+ std::printf ("\n");
+ }
+
+
void
Vector::print (std::ostream &out,
// get a representation of the vector and
// loop over all the elements
// TODO: up to now only local data
- // printed out!
+ // printed out! Find a way to neatly
+ // output distributed data...
TrilinosScalar *val;
int leading_dimension;
int ierr = vector->ExtractView (&val, &leading_dimension);
// Just swap the pointers to the
// two Epetra vectors that hold all
// the data.
- std::auto_ptr<Epetra_FEVector> tmp;
- tmp = v.vector;
- v.vector = vector;
- vector = tmp;
+ Vector *p_v = &v, *p_this = this;
+ Vector* tmp = p_v;
+ p_v = p_this;
+ p_this = tmp;
}