FEI Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fei_LinProbMgr_EpetraBasic.cpp
Go to the documentation of this file.
1 /*
2 // @HEADER
3 // ************************************************************************
4 // FEI: Finite Element Interface to Linear Solvers
5 // Copyright (2005) Sandia Corporation.
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
8 // U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Alan Williams (william@sandia.gov)
38 //
39 // ************************************************************************
40 // @HEADER
41 */
42 
43 
44 #include <fei_trilinos_macros.hpp>
45 #include <fei_SparseRowGraph.hpp>
46 
47 #ifdef HAVE_FEI_EPETRA
48 
50 
52  : comm_(comm),
53  ownedRows_(),
54  epetra_comm_(),
55  epetra_rowmap_(),
56  fei_srgraph_(),
57  crsgraph_(),
58  A_(),
59  numVectors_(1),
60  x_(),
61  b_()
62 {
63 #ifdef FEI_SER
65 #else
67 #endif
68 
69  epetra_comm_ = ecomm;
70 }
71 
73 {
74 }
75 
76 
78 ::setRowDistribution(const std::vector<int>& ownedGlobalRows)
79 {
80  if (ownedGlobalRows == ownedRows_) {
81  return;
82  }
83 
84  if (!ownedRows_.empty()) {
85  throw std::runtime_error("setRowDistribution called multiple times with different distributions. not allowed.");
86  }
87 
88  int* rows = const_cast<int*>(&ownedGlobalRows[0]);
89  epetra_rowmap_.reset(new Epetra_Map(-1, ownedGlobalRows.size(),
90  rows, 0, //indexBase
91  *epetra_comm_));
92 
93  x_.reset(new Epetra_MultiVector(*epetra_rowmap_, numVectors_));
94 
95  b_.reset(new Epetra_MultiVector(*epetra_rowmap_, numVectors_));
96 }
97 
100 {
101  if (fei_srgraph_.get() != NULL) {
102  if (*fei_srgraph_ != *matrixGraph) {
103  throw std::runtime_error("setMatrixGraph called multiple times with different graphs. not allowed.");
104  }
105  return;
106  }
107  else {
108  fei_srgraph_ = matrixGraph;
109  }
110 
111  if (epetra_rowmap_.get() == NULL) {
112  setRowDistribution(matrixGraph->rowNumbers);
113  }
114 
115  if ((int)fei_srgraph_->rowNumbers.size() != epetra_rowmap_->NumMyElements()) {
116  throw std::runtime_error("setMatrixGraph: num-rows not consistent with value from setRowDistribution");
117  }
118 
119  //We'll create and populate a Epetra_CrsGraph object.
120 
121  std::vector<int>& rowNumbers = fei_srgraph_->rowNumbers;
122  std::vector<int>& rowOffsets = fei_srgraph_->rowOffsets;
123 
124  //First create an array of num-indices-per-row.
125  std::vector<int> numIndicesPerRow; numIndicesPerRow.reserve(rowNumbers.size());
126 
127  int err;
128  unsigned i;
129  for(i=0; i<numIndicesPerRow.size(); ++i) {
130  numIndicesPerRow.push_back(rowOffsets[i+1] - rowOffsets[i]);
131  }
132 
133  bool static_profile = true;
134 
135  crsgraph_.reset(new Epetra_CrsGraph(Copy, *epetra_rowmap_,
136  &numIndicesPerRow[0], static_profile));
137 
138  //Now put in all the column-indices
139 
140  std::vector<int>& colIndices = fei_srgraph_->packedColumnIndices;
141  for(i=0; i<rowNumbers.size(); ++i) {
142  int offset = rowOffsets[i];
143  err = crsgraph_->InsertGlobalIndices(rowNumbers[i], numIndicesPerRow[i],
144  &colIndices[offset]);
145  if (err != 0) {
146  throw std::runtime_error("setMatrixGraph: err from Epetra_CrsGraph::InsertGlobalIndices.");
147  }
148  }
149 
150  err = crsgraph_->FillComplete();
151  if (err != 0) {
152  throw std::runtime_error("setMatrixGraph: err from Epetra_CrsGraph::FillComplete.");
153  }
154 
155  //and finally, create a matrix.
156  A_.reset(new Epetra_CrsMatrix(Copy, *crsgraph_));
157 }
158 
159 void LinProbMgr_EpetraBasic::setMatrixValues(double scalar)
160 {
161  int err = A_->PutScalar(scalar);
162  if (err != 0) {
163  throw std::runtime_error("error in Epetra_CrsMatrix->PutScalar");
164  }
165 }
166 
167 void LinProbMgr_EpetraBasic::setVectorValues(double scalar,
168  bool soln_vector)
169 {
170  int err = soln_vector ?
171  x_->PutScalar(scalar) : b_->PutScalar(scalar);
172  if (err != 0) {
173  throw std::runtime_error("error in Epetra_MultiVector->PutScalar");
174  }
175 }
176 
178 {
179  if (epetra_rowmap_.get() == NULL) return(-1);
180 
181  return(epetra_rowmap_->NumMyElements());
182 }
183 
185 {
186  if (A_.get() == NULL) return(-1);
187 
188  return( A_->NumGlobalEntries(row) );
189 }
190 
191 int LinProbMgr_EpetraBasic::copyOutMatrixRow(int row, int len,
192  double* coefs, int* indices)
193 {
194  int dummy;
195  return( A_->ExtractGlobalRowCopy(row, len, dummy, coefs, indices) );
196 }
197 
198 int LinProbMgr_EpetraBasic::insertMatrixValues(int numRows, const int* rows,
199  int numCols, const int* cols,
200  const double* const* values,
201  bool sum_into)
202 {
203  int* nc_cols = const_cast<int*>(cols);
204  double** nc_values = const_cast<double**>(values);
205  int err = 0;
206  if (sum_into) {
207  for(int i=0; i<numRows; ++i) {
208  err = A_->SumIntoGlobalValues(rows[i], numCols, nc_values[i], nc_cols);
209  if (err < 0) {
210  return(err);
211  }
212  }
213  }
214  else {
215  for(int i=0; i<numRows; ++i) {
216  err = A_->ReplaceGlobalValues(rows[i], numCols, nc_values[i], nc_cols);
217  if (err < 0) {
218  return(err);
219  }
220  }
221  }
222  return(err);
223 }
224 
226  const int* globalIndices,
227  const double* values,
228  bool sum_into,
229  bool soln_vector,
230  int vectorIndex)
231 {
232  double* localvaluesptr = soln_vector ?
233  x_->Pointers()[vectorIndex] : b_->Pointers()[vectorIndex];
234 
235  int min_my_gid = epetra_rowmap_->MinMyGID();
236  int returnValue = 0;
237 
238  if (sum_into) {
239  for(int i=0; i<numValues; ++i) {
240  int offset = globalIndices[i] - min_my_gid;
241  if (offset < 0) {
242  returnValue = 1;
243  continue;
244  }
245  localvaluesptr[offset] += values[i];
246  }
247  }
248  else {
249  for(int i=0; i<numValues; ++i) {
250  int offset = globalIndices[i] - min_my_gid;
251  if (offset < 0) {
252  returnValue = 1;
253  continue;
254  }
255  localvaluesptr[offset] = values[i];
256  }
257  }
258 
259  return(returnValue);
260 }
261 
263  const int* globalIndices,
264  double* values,
265  bool soln_vector,
266  int vectorIndex)
267 {
268  double* localvaluesptr = soln_vector ?
269  x_->Pointers()[vectorIndex] : b_->Pointers()[vectorIndex];
270 
271  int min_my_gid = epetra_rowmap_->MinMyGID();
272 
273  for(int i=0; i<numValues; ++i) {
274  int offset = globalIndices[i] - min_my_gid;
275  values[i] = localvaluesptr[offset];
276  }
277  return(0);
278 }
279 
280 double* LinProbMgr_EpetraBasic::getLocalVectorValuesPtr(bool soln_vector,
281  int vectorIndex)
282 {
283  double* localvaluesptr = soln_vector ?
284  x_->Pointers()[vectorIndex] : b_->Pointers()[vectorIndex];
285 
286  return(localvaluesptr);
287 }
288 
290 {
291  if (!A_->Filled()) {
292  //assumes the matrix is square!
293  int err = A_->FillComplete();
294  if (err != 0) {
295  return(err);
296  }
297  }
298 
299  if (!A_->StorageOptimized()) {
300  A_->OptimizeStorage();
301  }
302 
303  return(0);
304 }
305 
308 {
309  return( A_ );
310 }
311 
314 {
315  return( b_ );
316 }
317 
320 {
321  return( x_ );
322 }
323 
324 //HAVE_FEI_EPETRA
325 #endif
326 
void setMatrixGraph(fei::SharedPtr< fei::SparseRowGraph > matrixGraph)
int NumGlobalEntries(long long Row) const
void setRowDistribution(const std::vector< int > &ownedGlobalRows)
LinProbMgr_EpetraBasic(MPI_Comm comm)
bool StorageOptimized() const
int ExtractGlobalRowCopy(int_type Row, int Length, int &NumEntries, double *values, int_type *Indices) const
fei::SharedPtr< Epetra_MultiVector > x_
virtual ~LinProbMgr_EpetraBasic()
fei::SharedPtr< Epetra_CrsMatrix > get_A_matrix()
double * getLocalVectorValuesPtr(bool soln_vector, int vectorIndex=0)
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
std::vector< int > rowNumbers
int FillComplete(bool OptimizeDataStorage=true)
#define MPI_Comm
Definition: fei_mpi.h:56
int copyOutVectorValues(int numValues, const int *globalIndices, double *values, bool soln_vector, int vectorIndex=0)
fei::SharedPtr< Epetra_Map > epetra_rowmap_
int PutScalar(double ScalarConstant)
int NumMyElements() const
void setMatrixValues(double scalar)
int insertMatrixValues(int numRows, const int *rows, int numCols, const int *cols, const double *const *values, bool sum_into)
fei::SharedPtr< Epetra_MultiVector > get_rhs_vector()
T * get() const
int MinMyGID() const
int getRowLength(int row)
fei::SharedPtr< Epetra_CrsMatrix > A_
bool Filled() const
int copyOutMatrixRow(int row, int len, double *coefs, int *indices)
fei::SharedPtr< Epetra_MultiVector > b_
int insertVectorValues(int numValues, const int *globalIndices, const double *values, bool sum_into, bool soln_vector, int vectorIndex=0)
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
void setVectorValues(double scalar, bool soln_vector)
fei::SharedPtr< Epetra_MultiVector > get_solution_vector()