FEI Version of the Day
Loading...
Searching...
No Matches
fei_LinProbMgr_EpetraBasic.cpp
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
49#include <fei_LinProbMgr_EpetraBasic.hpp>
50
51LinProbMgr_EpetraBasic::LinProbMgr_EpetraBasic(MPI_Comm comm)
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
72LinProbMgr_EpetraBasic::~LinProbMgr_EpetraBasic()
73{
74}
75
76
77void LinProbMgr_EpetraBasic
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
98void LinProbMgr_EpetraBasic
99::setMatrixGraph(fei::SharedPtr<fei::SparseRowGraph> matrixGraph)
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
159void 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
167void 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
177int LinProbMgr_EpetraBasic::getLocalNumRows()
178{
179 if (epetra_rowmap_.get() == NULL) return(-1);
180
181 return(epetra_rowmap_->NumMyElements());
182}
183
184int LinProbMgr_EpetraBasic::getRowLength(int row)
185{
186 if (A_.get() == NULL) return(-1);
187
188 return( A_->NumGlobalEntries(row) );
189}
190
191int 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
198int 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
225int LinProbMgr_EpetraBasic::insertVectorValues(int numValues,
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
262int LinProbMgr_EpetraBasic::copyOutVectorValues(int numValues,
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
280double* 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
289int LinProbMgr_EpetraBasic::globalAssemble()
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
307LinProbMgr_EpetraBasic::get_A_matrix()
308{
309 return( A_ );
310}
311
313LinProbMgr_EpetraBasic::get_rhs_vector()
314{
315 return( b_ );
316}
317
319LinProbMgr_EpetraBasic::get_solution_vector()
320{
321 return( x_ );
322}
323
324//HAVE_FEI_EPETRA
325#endif
326