#include <lac/vector.h>
#include <lac/block_indices.h>
#include <lac/block_vector_base.h>
+#include <lac/trilinos_block_vector.h>
#include <cstdio>
#include <vector>
DEAL_II_NAMESPACE_OPEN
+#ifdef DEAL_II_USE_TRILINOS
+namespace TrilinosWrappers
+{
+ class Vector;
+}
+#endif
+
+
+
/*! @addtogroup Vectors
*@{
*/
BlockVector (const BlockVector<OtherNumber> &v);
#endif
+
+#ifdef DEAL_II_USE_TRILINOS
+ /**
+ * A copy constructor taking a
+ * (parallel) Trilinos block
+ * vector and copying it into
+ * the deal.II own format.
+ */
+ BlockVector (const TrilinosWrappers::BlockVector &v);
+
+#endif
/**
* Constructor. Set the number of
* blocks to
BlockVector &
operator= (const Vector<Number> &V);
+#ifdef DEAL_II_USE_TRILINOS
+ /**
+ * A copy constructor from a
+ * Trilinos block vector to a
+ * deal.II block vector.
+ */
+ BlockVector &
+ operator= (const TrilinosWrappers::BlockVector &V);
+#endif
/**
* Reinitialize the BlockVector to
* contain <tt>num_blocks</tt> blocks of
/*----------------------- Inline functions ----------------------------------*/
+
template <typename Number>
template <typename InputIterator>
BlockVector<Number>::BlockVector (const std::vector<unsigned int> &n,
}
+#ifdef DEAL_II_USE_TRILINOS
+template <typename Number>
+inline
+BlockVector<Number> &
+BlockVector<Number>::operator = (const TrilinosWrappers::BlockVector &v)
+{
+ BaseClass::operator = (v);
+ return *this;
+}
+#endif
+
template <typename Number>
void BlockVector<Number>::scale (const value_type factor)
{
row_map (0, 0, Epetra_SerialComm()),
#endif
col_map (row_map),
- last_action (Insert),
+ last_action (Zero),
matrix (std::auto_ptr<Epetra_FECrsMatrix>
(new Epetra_FECrsMatrix(Copy, row_map, 0)))
{}
:
row_map (InputMap),
col_map (row_map),
- last_action (Insert),
+ last_action (Zero),
matrix (std::auto_ptr<Epetra_FECrsMatrix>
(new Epetra_FECrsMatrix(Copy, row_map,
int(n_max_entries_per_row), false)))
:
row_map (InputMap),
col_map (row_map),
- last_action (Insert),
+ last_action (Zero),
matrix (std::auto_ptr<Epetra_FECrsMatrix>
(new Epetra_FECrsMatrix(Copy, row_map,
(int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
:
row_map (InputRowMap),
col_map (InputColMap),
- last_action (Insert),
+ last_action (Zero),
matrix (std::auto_ptr<Epetra_FECrsMatrix>
(new Epetra_FECrsMatrix(Copy, row_map, col_map,
int(n_max_entries_per_row), false)))
:
row_map (InputRowMap),
col_map (InputColMap),
- last_action (Insert),
+ last_action (Zero),
matrix (std::auto_ptr<Epetra_FECrsMatrix>
(new Epetra_FECrsMatrix(Copy, row_map, col_map,
(int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
// sparsity pattern on that
// processor, and only broadcast the
// pattern afterwards.
- if (row_map.Comm().MyPID() == 0)
- {
- Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
- ExcDimensionMismatch (matrix->NumGlobalRows(),
- sparsity_pattern.n_rows()));
+ Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
+ ExcDimensionMismatch (matrix->NumGlobalRows(),
+ sparsity_pattern.n_rows()));
// Trilinos seems to have a bug for
// rectangular matrices at this point,
// ExcDimensionMismatch (matrix->NumGlobalCols(),
// sparsity_pattern.n_cols()));
- std::vector<double> values;
- std::vector<int> row_indices;
+ std::vector<double> values;
+ std::vector<int> row_indices;
- for (unsigned int row=0; row<n_rows; ++row)
- {
- const int row_length = sparsity_pattern.row_length(row);
- row_indices.resize (row_length, -1);
- values.resize (row_length, 0.);
-
- for (int col=0; col < row_length; ++col)
- row_indices[col] = sparsity_pattern.column_number (row, col);
-
- matrix->InsertGlobalValues (row, row_length,
- &values[0], &row_indices[0]);
- }
- }
+ for (unsigned int row=0; row<n_rows; ++row)
+ if (row_map.MyGID(row))
+ {
+ const int row_length = sparsity_pattern.row_length(row);
+ row_indices.resize (row_length, -1);
+ values.resize (row_length, 0.);
+
+ for (int col=0; col < row_length; ++col)
+ row_indices[col] = sparsity_pattern.column_number (row, col);
+
+ matrix->InsertGlobalValues (row, row_length,
+ &values[0], &row_indices[0]);
+ }
- last_action = Insert;
+ last_action = Zero;
// In the end, the matrix needs to
// be compressed in order to be
n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
// TODO: There seems to be problem
- // in Epetra when a quadratic matrix
- // is initialized with both row and
+ // in Epetra when a matrix is
+ // initialized with both row and
// column map. Maybe find something
- // more out about this...
- if (row_map.SameAs(col_map) == true)
- matrix = std::auto_ptr<Epetra_FECrsMatrix>
+ // more out about this... It could
+ // be related to the bug #4123. For
+ // the moment, just ignore the
+ // column information and generate
+ // the matrix as if it were
+ // square. The call to
+ // GlobalAssemble later will set the
+ // correct values.
+ matrix = std::auto_ptr<Epetra_FECrsMatrix>
(new Epetra_FECrsMatrix(Copy, row_map,
&n_entries_per_row[0], false));
- else
- matrix = std::auto_ptr<Epetra_FECrsMatrix>
- (new Epetra_FECrsMatrix(Copy, row_map, col_map,
- &n_entries_per_row[0], false));
reinit (sparsity_pattern);
}
n_entries_per_row[(int)row] =
dealii_sparse_matrix.get_sparsity_pattern().row_length(row);
+ // TODO: There seems to be problem
+ // in Epetra when a matrix is
+ // initialized with both row and
+ // column map. Maybe find something
+ // more out about this... It could
+ // be related to the bug #4123. For
+ // the moment, just ignore the
+ // column information and generate
+ // the matrix as if it were
+ // square. The call to
+ // GlobalAssemble later will set the
+ // correct values.
matrix = std::auto_ptr<Epetra_FECrsMatrix>
- (new Epetra_FECrsMatrix(Copy, row_map, col_map,
- &n_entries_per_row[0], true));
+ (new Epetra_FECrsMatrix(Copy, row_map,
+ &n_entries_per_row[0], false));
std::vector<double> values;
std::vector<int> row_indices;
Assert (numbers::is_finite(value),
ExcMessage("The given value is not finite but either "
- "infinite or Not A Number (NaN)"));
+ "infinite or Not A Number (NaN)"));
if (last_action == Insert)
{
ierr = matrix->GlobalAssemble (false);
else
ierr = matrix->GlobalAssemble(col_map, row_map, false);
-
+
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
last_action = Add;
- }
+ }
// we have to do above actions in any
// case to be consistent with the MPI
const_cast<double*>(&value),
&trilinos_j);
+ if (ierr > 0)
+ std::cout << ierr << " " << m() << " " << n() << std::endl;
AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j));
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
}