* elements.
*/
BlockVector (const BlockVector &V);
-
+
/**
* Creates a block vector
* consisting of
* reinit of the blocks.
*/
BlockVector (const unsigned int num_blocks);
-
+
+ /**
+ * Constructor. Set the number of
+ * blocks to <tt>n.size()</tt> and
+ * initialize each block with
+ * <tt>n[i]</tt> zero elements.
+ *
+ * References BlockVector.reinit().
+ */
+ BlockVector (const std::vector<unsigned int> &N);
+
+ /**
+ * Constructor. Set the number of
+ * blocks to
+ * <tt>n.size()</tt>. Initialize the
+ * vector with the elements
+ * pointed to by the range of
+ * iterators given as second and
+ * third argument. Apart from the
+ * first argument, this
+ * constructor is in complete
+ * analogy to the respective
+ * constructor of the
+ * <tt>std::vector</tt> class, but the
+ * first argument is needed in
+ * order to know how to subdivide
+ * the block vector into
+ * different blocks.
+ */
+ template <typename InputIterator>
+ BlockVector (const std::vector<unsigned int> &n,
+ const InputIterator first,
+ const InputIterator end);
+
/**
* Destructor. Clears memory
*/
+ inline
+ BlockVector::BlockVector (const std::vector<unsigned int> &N)
+ {
+ reinit (N);
+ }
+
+
+
+ template <typename InputIterator>
+ BlockVector::BlockVector (const std::vector<unsigned int> &n,
+ const InputIterator first,
+ const InputIterator end)
+ {
+ // first set sizes of blocks, but
+ // don't initialize them as we will
+ // copy elements soon
+ reinit (n, true);
+ InputIterator start = first;
+ for (unsigned int b=0; b<n.size(); ++b)
+ {
+ InputIterator end = start;
+ std::advance (end, static_cast<signed int>(n[b]));
+
+ for (unsigned int i=0; i<n[b]; ++i, ++start)
+ this->block(b)(i) = *start;
+ }
+ Assert (start == end, ExcIteratorRangeDoesNotMatchVectorSize());
+ }
+
+
+
inline
BlockVector::BlockVector (const unsigned int num_blocks)
{
const unsigned int n_max_entries_per_row);
/**
- * This constructor is similar to
- * the one above, but it now
- * takes two different Epetra
- * maps for rows and
- * columns. This interface is
- * meant to be used for
- * generating rectangular
- * matrices, where one map
- * specifies the parallel
+ * This constructor is similar to the
+ * one above, but it now takes two
+ * different Epetra maps for rows and
+ * columns. This interface is meant
+ * to be used for generating
+ * rectangular matrices, where one
+ * map specifies the parallel
* distribution of rows and the
- * second one specifies the
- * number of columns in the
- * total matrix. It also provides
- * information for the internal
- * arrangement in matrix vector
- * products, but is not used for
- * the distribution of the
- * columns – rather, all
- * column elements of a row are
- * stored on the same
- * processor. The vector
- * n_entries_per_row specifies
- * the number of entries in each
- * row of the newly generated
+ * second one specifies the number of
+ * columns in the total matrix. It
+ * also provides information for the
+ * internal arrangement in matrix
+ * vector products, but is not used
+ * for the distribution of the
+ * columns – rather, all column
+ * elements of a row are stored on
+ * the same processor. The vector
+ * <tt>n_entries_per_row</tt>
+ * specifies the number of entries in
+ * each row of the newly generated
* matrix.
*/
SparseMatrix (const Epetra_Map &InputRowMap,
* #n columns. The resulting matrix
* will be completely stored locally.
*
- * The number of columns entries
- * per row is specified as the
- * maximum number of entries
- * argument.
+ * The number of columns entries per
+ * row is specified as the maximum
+ * number of entries argument.
*/
SparseMatrix (const unsigned int m,
const unsigned int n,
const unsigned int n_max_entries_per_row);
+ /**
+ * Generate a matrix that is completely
+ * stored locally, having #m rows and
+ * #n columns. The resulting matrix
+ * will be completely stored locally.
+ *
+ * The vector
+ * <tt>n_entries_per_row</tt>
+ * specifies the number of entries in
+ * each row.
+ */
+ SparseMatrix (const unsigned int m,
+ const unsigned int n,
+ const std::vector<unsigned int> &n_entries_per_row);
+
/**
* Copy constructor. Sets the
* calling matrix to be the same
:
map (InputMap)
{
- vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map, false));
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map));
const int min_my_id = map.MinMyGID();
- const unsigned int size = map.NumMyElements();
+ const int size = map.NumMyElements();
Assert (map.MaxLID() == size-1,
ExcDimensionMismatch(map.MaxLID(), size-1));
+ // Need to copy out values, since the
+ // deal.II might not use doubles, so
+ // that a direct access is not possible.
std::vector<int> indices (size);
+ std::vector<double> values (size);
for (unsigned int i=0; i<size; ++i)
- indices[i] = map.GID(i);
+ {
+ indices[i] = map.GID(i);
+ values[i] = v(i);
+ }
+ const int ierr = vector->ReplaceGlobalValues (size, &indices[0],
+ &values[0]);
- const int ierr = vector->ReplaceGlobalValues(size, v.begin()+min_my_id,
- &indices[0]);
AssertThrow (ierr == 0, ExcTrilinosError());
}
map (v.size(), 0, Epetra_SerialComm())
#endif
{
- vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map, false));
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map));
std::vector<int> indices (v.size());
+ std::vector<double> values (v.size());
for (unsigned int i=0; i<v.size(); ++i)
- indices[i] = map.GID(i);
+ {
+ indices[i] = map.GID(i);
+ values[i] = v(i);
+ }
+
+ const int ierr = vector->ReplaceGlobalValues((int)v.size(),
+ &indices[0], &values[0]);
- const int ierr = vector->ReplaceGlobalValues(v.size(), v.begin(),
- &indices[0]);
- AssertThrow (ierr == 0, ExcTrilinosError());
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
}
int(n_max_entries_per_row), false)))
{}
+ SparseMatrix::SparseMatrix (const unsigned int m,
+ const unsigned int n,
+ const std::vector<unsigned int> &n_entries_per_row)
+ :
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ row_map (m, 0, Epetra_MpiComm(MPI_COMM_WORLD)),
+ col_map (n, 0, Epetra_MpiComm(MPI_COMM_WORLD)),
+#else
+ row_map (m, 0, Epetra_SerialComm()),
+ col_map (n, 0, Epetra_SerialComm()),
+#endif
+ last_action (Zero),
+ compressed (true),
+ matrix (std::auto_ptr<Epetra_FECrsMatrix>
+ (new Epetra_FECrsMatrix(Copy, row_map,
+ (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
+ false)))
+ {}
+
SparseMatrix::SparseMatrix (const SparseMatrix &InputMatrix)
:
Subscriptor(),
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 the last thing we did was an
+ // addition, we need to do all the
+ // data exchange.
if (last_action == Add)
{
int ierr;
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;
// 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);
+ {
+ ierr = matrix->InsertGlobalValues (trilinos_i, 1,
+ const_cast<double*>(&value),
+ &trilinos_j);
+ if (ierr > 0)
+ ierr = 0;
+ }
else
- ierr = matrix->ReplaceGlobalValues (trilinos_i, 1,
- const_cast<double*>(&value),
- &trilinos_j);
+ {
+ ierr = matrix->ReplaceGlobalValues (trilinos_i, 1,
+ const_cast<double*>(&value),
+ &trilinos_j);
+ }
AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j));
AssertThrow (ierr == 0, ExcTrilinosError(ierr));