typedef const internal::VectorReference const_reference;
/**
- * Default constructor. It
- * doesn't do anything, derived
- * classes will have to
- * initialize the data.
+ * Default constructor that
+ * generates an empty (zero
+ * size) vector. The function
+ * <tt>reinit()</tt> will have
+ * to give the vector the
+ * correct size and
+ * distribution among processes
+ * in case of an MPI run.
*/
Vector ();
* Epetra_Map that already
* knows how to distribute the
* individual components among
- * the MPI processors,
- * including the size of the
- * vector.
+ * the MPI processors. It also
+ * includes information about
+ * the size of the vector.
*/
Vector (const Epetra_Map &InputMap);
*/
virtual ~Vector ();
- /**
+ /**
* Reinit functionality. This function
* destroys the old vector content
* and generates a new one based on
*/
void reinit (const Epetra_Map &input_map);
- /**
+ /**
* Reinit functionality. This function
* copies the vector v to the current
* one.
// ExcDimensionMismatch (matrix->NumGlobalCols(),
// sparsity_pattern.n_cols()));
- std::vector<int> n_entries_per_row(n_rows);
-
- for (unsigned int row=0; row<n_rows; ++row)
- n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
-
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, 0);
+ row_indices.resize (row_length, -1);
values.resize (row_length, 0.);
- for (int col=0; col< row_length; ++col)
+ 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]);
+ matrix->InsertGlobalValues (row, row_length,
+ &values[0], &row_indices[0]);
}
}
SparseMatrix::reinit (const Epetra_Map &input_map,
const SparsityPattern &sparsity_pattern)
{
- // TODO: There seems to be problem
- // in Epetra when a quadratic matrix
- // is initialized with both row and
- // column map. Maybe find something
- // more out about this...
reinit (input_map, input_map, sparsity_pattern);
-
- /* matrix.reset();
-
- unsigned int n_rows = sparsity_pattern.n_rows();
-
- if (row_map.Comm().MyPID() == 0)
- {
- Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(),
- ExcDimensionMismatch (input_map.NumGlobalElements(),
- sparsity_pattern.n_rows()));
- Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(),
- ExcDimensionMismatch (input_map.NumGlobalElements(),
- sparsity_pattern.n_cols()));
- }
-
- row_map = input_map;
- col_map = row_map;
-
- std::vector<int> n_entries_per_row(n_rows);
-
- for (unsigned int row=0; row<n_rows; ++row)
- n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
-
- matrix = std::auto_ptr<Epetra_FECrsMatrix>
- (new Epetra_FECrsMatrix(Copy, row_map, &n_entries_per_row[0],
- false));
-
- reinit (sparsity_pattern);*/
}
for (unsigned int row=0; row<n_rows; ++row)
n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
- matrix = std::auto_ptr<Epetra_FECrsMatrix>
- (new Epetra_FECrsMatrix(Copy, row_map, col_map,
- &n_entries_per_row[0], false));
+ // TODO: There seems to be problem
+ // in Epetra when a quadratic 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>
+ (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);
}
{
// flush buffers
int ierr;
- //if (row_map.SameAs(col_map))
- //ierr = matrix->GlobalAssemble ();
- //else
+ if (row_map.SameAs(col_map))
+ ierr = matrix->GlobalAssemble ();
+ else
ierr = matrix->GlobalAssemble (col_map, row_map);
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
// If the data is not on the
// present processor, we can't
- // continue.
+ // continue. Just print out
+ // zero.
+
+ // TODO: Is this reasonable?
if ((trilinos_i == -1 ) || (trilinos_j == -1))
{
- Assert (false, ExcAccessToNonLocalElement(i, j, local_range().first,
- local_range().second));
+ return 0.;
+ //Assert (false, ExcAccessToNonLocalElement(i, j, local_range().first,
+ // local_range().second));
}
else
{
unsigned int
SparseMatrix::m () const
{
- int n_rows = matrix->NumGlobalRows();
+ int n_rows = matrix -> NumGlobalRows();
return n_rows;
}
SparseMatrix::local_range () const
{
int begin, end;
- begin = matrix->RowMap().MinMyGID();
- end = matrix->RowMap().MaxMyGID()+1;
+ begin = matrix -> RowMap().MinMyGID();
+ end = matrix -> RowMap().MaxMyGID()+1;
return std::make_pair (begin, end);
}