preset_nonzero_locations);
}
+ template <typename SparsityType>
+ void
+ SparseMatrix::
+ reinit (const IndexSet & local_rows,
+ const IndexSet & local_columns,
+ const SparsityType &sparsity_pattern,
+ const MPI_Comm &communicator)
+ {
+ this->communicator = communicator;
+ // get rid of old matrix and generate a
+ // new one
+#if DEAL_II_PETSC_VERSION_LT(3,2,0)
+ const int ierr = MatDestroy (matrix);
+#else
+ const int ierr = MatDestroy (&matrix);
+#endif
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ do_reinit (local_rows, local_columns, sparsity_pattern);
+ }
void
SparseMatrix::do_reinit (const unsigned int m,
}
+ template <typename SparsityType>
+ void
+ SparseMatrix::
+ do_reinit (const IndexSet & local_rows,
+ const IndexSet & local_columns,
+ const SparsityType &sparsity_pattern)
+ {
+ Assert(sparsity_pattern.n_rows()==local_rows.size(),
+ ExcMessage("SparsityPattern and IndexSet have different number of rows"));
+ Assert(sparsity_pattern.n_cols()==local_columns.size(),
+ ExcMessage("SparsityPattern and IndexSet have different number of columns"));
+ Assert(local_rows.is_contiguous() && local_columns.is_contiguous(),
+ ExcMessage("PETSc only supports contiguous row/column ranges"));
+
+
+
+
+ // create the matrix. We do not set row length but set the
+ // correct SparsityPattern later.
+ int ierr;
+
+ ierr = MatCreate(communicator,&matrix);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ ierr = MatSetSizes(matrix,
+ local_rows.n_elements(),
+ local_columns.n_elements(),
+ sparsity_pattern.n_rows(),
+ sparsity_pattern.n_cols());
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ ierr = MatSetType(matrix,MATMPIAIJ);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+
+ // next preset the exact given matrix
+ // entries with zeros. this doesn't avoid any
+ // memory allocations, but it at least
+ // avoids some searches later on. the
+ // key here is that we can use the
+ // matrix set routines that set an
+ // entire row at once, not a single
+ // entry at a time
+ //
+ // for the usefulness of this option
+ // read the documentation of this
+ // class.
+ //if (preset_nonzero_locations == true)
+ {
+ // MatMPIAIJSetPreallocationCSR
+ // can be used to allocate the sparsity
+ // pattern of a matrix
+
+ const PetscInt local_row_start = local_rows.nth_index_in_set(0);
+ const PetscInt local_col_start = local_columns.nth_index_in_set(0);
+ const PetscInt
+ local_row_end = local_row_start + local_rows.n_elements();
+
+
+ // first set up the column number
+ // array for the rows to be stored
+ // on the local processor. have one
+ // dummy entry at the end to make
+ // sure petsc doesn't read past the
+ // end
+ std::vector<PetscInt>
+
+ rowstart_in_window (local_row_end - local_row_start + 1, 0),
+ colnums_in_window;
+ {
+ unsigned int n_cols = 0;
+ for (PetscInt i=local_row_start; i<local_row_end; ++i)
+ {
+ const PetscInt row_length = sparsity_pattern.row_length(i);
+ rowstart_in_window[i+1-local_row_start]
+ = rowstart_in_window[i-local_row_start] + row_length;
+ n_cols += row_length;
+ }
+ colnums_in_window.resize (n_cols+1, -1);
+ }
+
+ // now copy over the information
+ // from the sparsity pattern.
+ {
+ PetscInt* ptr = & colnums_in_window[0];
+
+ for (PetscInt i=local_row_start; i<local_row_end; ++i)
+ {
+ typename SparsityType::row_iterator
+ row_start = sparsity_pattern.row_begin(i),
+ row_end = sparsity_pattern.row_end(i);
+
+ std::copy(row_start, row_end, ptr);
+ ptr += row_end - row_start;
+ }
+ }
+
+
+ // then call the petsc function
+ // that summarily allocates these
+ // entries:
+ MatMPIAIJSetPreallocationCSR (matrix,
+ &rowstart_in_window[0],
+ &colnums_in_window[0],
+ 0);
+ compress ();
+
+
+ // Tell PETSc that we are not
+ // planning on adding new entries
+ // to the matrix. Generate errors
+ // in debug mode.
+ int ierr;
+ #if DEAL_II_PETSC_VERSION_LT(3,0,0)
+ #ifdef DEBUG
+ ierr = MatSetOption (matrix, MAT_NEW_NONZERO_LOCATION_ERR);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+ #else
+ ierr = MatSetOption (matrix, MAT_NO_NEW_NONZERO_LOCATIONS);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+ #endif
+ #else
+ #ifdef DEBUG
+ ierr = MatSetOption (matrix, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+ #else
+ ierr = MatSetOption (matrix, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+ #endif
+ #endif
+
+ // Tell PETSc to keep the
+ // SparsityPattern entries even if
+ // we delete a row with
+ // clear_rows() which calls
+ // MatZeroRows(). Otherwise one can
+ // not write into that row
+ // afterwards.
+ #if DEAL_II_PETSC_VERSION_LT(3,0,0)
+ ierr = MatSetOption (matrix, MAT_KEEP_ZEROED_ROWS);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+ #elif DEAL_II_PETSC_VERSION_LT(3,1,0)
+ ierr = MatSetOption (matrix, MAT_KEEP_ZEROED_ROWS, PETSC_TRUE);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+ #else
+ ierr = MatSetOption (matrix, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+ #endif
+
+ }
+
+ }
+
template <typename SparsityType>
void
const unsigned int,
const bool);
+ template void
+ SparseMatrix::
+ reinit (const IndexSet &,
+ const IndexSet &,
+ const CompressedSimpleSparsityPattern &,
+ const MPI_Comm &);
+
template void
SparseMatrix::do_reinit (const SparsityPattern &,
const std::vector<unsigned int> &,
const unsigned int ,
const bool);
+ template void
+ SparseMatrix::
+ do_reinit (const IndexSet &,
+ const IndexSet &,
+ const CompressedSimpleSparsityPattern &);
+
PetscScalar
SparseMatrix::matrix_norm_square (const Vector &v) const