FEI Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fei_MatrixTraits_Epetra.hpp
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 #ifndef _fei_MatrixTraits_Epetra_h_
45 #define _fei_MatrixTraits_Epetra_h_
46 
47 #include <fei_trilinos_macros.hpp>
48 
49 #ifdef HAVE_FEI_EPETRA
50 
51 //
52 //IMPORTANT NOTE: Make sure that wherever this file is included from, it
53 //appears BEFORE any include of fei_Vector_Impl.hpp or fei_Matrix_Impl.hpp !!!
54 //
55 #include <fei_MatrixTraits.hpp>
58 #include <fei_Include_Trilinos.hpp>
59 #include <fei_Vector_Impl.hpp>
60 
61 namespace fei {
68  template<>
69  struct MatrixTraits<Epetra_CrsMatrix> {
70  static const char* typeName()
71  { return("Epetra_CrsMatrix"); }
72 
73  static double* getBeginPointer(Epetra_CrsMatrix* mat)
74  {
75  return (*mat)[0];
76  }
77 
78  static int getOffset(Epetra_CrsMatrix* A, int row, int col)
79  {
80  const Epetra_Map& erowmap = A->RowMap();
81  const Epetra_Map& ecolmap = A->ColMap();
82  int local_row = erowmap.LID(row);
83  int local_col = ecolmap.LID(col);
84 
85  int* rowOffsets;
86  int* colIndices;
87  double* coefs;
88  A->ExtractCrsDataPointers(rowOffsets, colIndices, coefs);
89 
90  int* row_ptr = &colIndices[rowOffsets[local_row]];
91  int* end_row = &colIndices[rowOffsets[local_row+1]];
92 
93  int col_offset = 0;
94  for(; row_ptr != end_row; ++row_ptr) {
95  if (*row_ptr == local_col) break;
96  ++col_offset;
97  }
98 
99  return rowOffsets[local_row] + col_offset;
100  }
101 
102  static int setValues(Epetra_CrsMatrix* mat, double scalar)
103  {
104  return( mat->PutScalar(scalar) );
105  }
106 
107  static int getNumLocalRows(Epetra_CrsMatrix* mat, int& numRows)
108  {
109  numRows = mat->NumMyRows();
110  return(0);
111  }
112 
113  static int getRowLength(Epetra_CrsMatrix* mat, int row, int& length)
114  {
115  length = mat->NumGlobalEntries(row);
116  if (length < 0) return(-1);
117  return( 0 );
118  }
119 
120  static int copyOutRow(Epetra_CrsMatrix* mat,
121  int row, int len, double* coefs, int* indices)
122  {
123  int dummy;
124  return(mat->ExtractGlobalRowCopy(row, len, dummy, coefs, indices));
125  }
126 
127  static int putValuesIn(Epetra_CrsMatrix* mat,
128  int numRows, const int* rows,
129  int numCols, const int* cols,
130  const double* const* values,
131  bool sum_into)
132  {
134  static std::vector<int> idx;
135  idx.resize(numRows+numCols);
136  int* idx_row = &idx[0];
137  int* idx_col = idx_row+numRows;
138  for(int i=0; i<numRows; ++i) {
139  idx_row[i] = mat->RowMap().LID(rows[i]);
140  }
141  for(int i=0; i<numCols; ++i) {
142  idx_col[i] = mat->ColMap().LID(cols[i]);
143  }
144  if (sum_into) {
145  for(int i=0; i<numRows; ++i) {
146  int err = mat->SumIntoMyValues(idx_row[i], numCols,
147  (double*)values[i],
148  idx_col);
149  if (err != 0) {
150  return(err);
151  }
152  }
153  }
154  else {
155  for(int i=0; i<numRows; ++i) {
156  int err = mat->ReplaceMyValues(idx_row[i], numCols,
157  (double*)values[i],
158  idx_col);
159  if (err != 0) {
160  return(err);
161  }
162  }
163  }
164  return(0);
165  }
166 
167  static int globalAssemble(Epetra_CrsMatrix* mat)
168  {
169  if (!mat->Filled()) {
170  int err = mat->FillComplete();
171  if (err != 0) {
172  fei::console_out() << "MatrixTraits<Epetra_CrsMatrix>::globalAssemble"
173  << " ERROR in mat->FillComplete" << FEI_ENDL;
174  return(-1);
175  }
176  }
177 
178  if (!mat->StorageOptimized()) {
179  mat->OptimizeStorage();
180  }
181 
182  return( 0 );
183  }
184 
185  static int matvec(Epetra_CrsMatrix* mat,
186  fei::Vector* x,
187  fei::Vector* y)
188  {
190  dynamic_cast<fei::Vector_Impl<Epetra_MultiVector>* >(x);
192  dynamic_cast<fei::Vector_Impl<Epetra_MultiVector>* >(y);
193 
194  if (evx == NULL || evy == NULL) {
195  return(-1);
196  }
197 
200 
201  return( mat->Multiply(false, *ex, *ey) );
202  }
203 
204  };//struct MatrixTraits<Epetra_CrsMatrix>
205 }//namespace fei
206 
207 namespace snl_fei {
214  template<>
215  struct BlockMatrixTraits<Epetra_VbrMatrix> {
216  static const char* typeName()
217  { return("Epetra_VbrMatrix"); }
218 
219  static int putScalar(Epetra_VbrMatrix* mat, double scalar)
220  {
221  return( mat->PutScalar(scalar) );
222  }
223 
224  static int getRowLength(Epetra_VbrMatrix* mat, int row, int& length)
225  {
226  length = mat->NumGlobalBlockEntries(row);
227  return(0);
228  }
229 
230  static int getPointRowLength(Epetra_VbrMatrix* mat, int row, int& length)
231  {
232  const Epetra_Map& map = mat->RowMatrixRowMap();
233  int minLocalRow = map.MinMyGID();
234  int localRow = row - minLocalRow;
235  int error = mat->NumMyRowEntries(localRow, length);
236  return(error);
237  }
238 
239  static int copyOutRow(Epetra_VbrMatrix* mat,
240  int row, int numBlkCols,
241  int rowDim,
242  int* blkCols,
243  int* colDims,
244  double* coefs,
245  int coefsLen,
246  int& blkRowLength)
247  {
248  int checkRowDim;
249  int error = mat->BeginExtractGlobalBlockRowCopy(row, numBlkCols,
250  checkRowDim,
251  blkRowLength,
252  blkCols, colDims);
253  if (error != 0 || checkRowDim != rowDim || blkRowLength != numBlkCols) {
254  return(error);
255  }
256 
257  int offset = 0;
258  for(int i=0; i<numBlkCols; ++i) {
259  if (offset >= coefsLen) {
260  std::cerr << "BlockMatrixTraits::copyOutRow ran off end of coefs array."
261  << std::endl;
262  return(-2);
263  }
264  int numValues = rowDim*colDims[i];
265  error = mat->ExtractEntryCopy(numValues, &(coefs[offset]),
266  rowDim, false);
267  if (error != 0) {
268  return(error);
269  }
270  offset += numValues;
271  }
272 
273  return(0);
274  }
275 
276  static int copyOutPointRow(Epetra_VbrMatrix* mat,
277  int firstLocalOffset,
278  int row,
279  int len,
280  double* coefs,
281  int* indices,
282  int& rowLength)
283  {
284  int error = mat->ExtractMyRowCopy(row-firstLocalOffset,
285  len, rowLength,
286  coefs, indices);
287 
288  const Epetra_Map& colmap = mat->RowMatrixColMap();
289  for(int i=0; i<len; ++i) {
290  indices[i] = colmap.GID(indices[i]);
291  }
292 
293  return(error);
294  }
295 
296  static int sumIn(Epetra_VbrMatrix* mat,
297  int blockRow,
298  int rowDim,
299  int numBlockCols,
300  const int* blockCols,
301  const int* colDims,
302  int LDA,
303  const double* values)
304  {
305  int err, voffset = 0;
306  for(int j=0; j<numBlockCols; ++j) {
307  err = mat->DirectSubmitBlockEntry(blockRow, blockCols[j],
308  &(values[voffset]), LDA,
309  rowDim, colDims[j], true/*sum_into*/);
310  if (err != 0) return(err);
311 
312  voffset += colDims[j]*LDA;
313  }
314 
315  return(0);
316  }
317 
318  static int copyIn(Epetra_VbrMatrix* mat,
319  int blockRow,
320  int rowDim,
321  int numBlockCols,
322  const int* blockCols,
323  const int* colDims,
324  int LDA,
325  const double* values)
326  {
327  int err, voffset = 0;
328  for(int j=0; j<numBlockCols; ++j) {
329  err = mat->DirectSubmitBlockEntry(blockRow, blockCols[j],
330  &(values[voffset]), LDA,
331  rowDim, colDims[j], false/*replace*/);
332  if (err != 0) return(err);
333 
334  voffset += colDims[j]*LDA;
335  }
336 
337  return(0);
338  }
339 
340  static int sumIn(Epetra_VbrMatrix* mat,
341  int row,
342  int rowDim,
343  int numCols,
344  const int* cols,
345  const int* LDAs,
346  const int* colDims,
347  const double* const* values)
348  {
349  int err = 0;
350  for(int i=0; i<numCols; ++i) {
351  err = mat->DirectSubmitBlockEntry(row, cols[i], values[i],
352  LDAs[i], rowDim, colDims[i], true);
353  if (err != 0) return(err);
354  }
355 
356  return(err);
357  }
358 
359  static int copyIn(Epetra_VbrMatrix* mat,
360  int row,
361  int rowDim,
362  int numCols,
363  const int* cols,
364  const int* LDAs,
365  const int* colDims,
366  const double* const* values)
367  {
368  int err = 0;
369  for(int i=0; i<numCols; ++i) {
370  err = mat->DirectSubmitBlockEntry(row, cols[i], values[i],
371  LDAs[i], rowDim, colDims[i], false);
372  if (err != 0) return(err);
373  }
374 
375  return(err);
376  }
377 
378  static int globalAssemble(Epetra_VbrMatrix* mat)
379  {
380  return( mat->FillComplete() );
381  }
382  };//struct BlockMatrixTraits<Epetra_VbrMatrix>
383 }//namespace snl_fei
384 #endif //HAVE_FEI_EPETRA
385 #endif // _fei_MatrixTraits_Epetra_hpp_
int NumGlobalEntries(long long Row) const
bool StorageOptimized() const
virtual const Epetra_Map & RowMatrixRowMap() const =0
static int getPointRowLength(T *, int, int &)
int ExtractGlobalRowCopy(int_type Row, int Length, int &NumEntries, double *values, int_type *Indices) const
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
static int copyIn(T *, int, int, int, const int *, const int *, int, const double *)
static const char * typeName()
static double * getBeginPointer(T *)
static int matvec(T *A, fei::Vector *x, fei::Vector *y)
const Epetra_Map & ColMap() const
int FillComplete(bool OptimizeDataStorage=true)
int PutScalar(double ScalarConstant)
const Epetra_Map & RowMap() const
virtual int NumMyRowEntries(int MyRow, int &NumEntries) const =0
static int getRowLength(T *mat, int row, int &length)
static int copyOutPointRow(T *, int, int, int, double *, int *, int &)
static int getNumLocalRows(T *mat, int &numRows)
int GID(int LID) const
int NumMyRows() const
static int getRowLength(T *, int, int &)
int LID(int GID) const
std::string error
int MinMyGID() const
#define FEI_ENDL
int SumIntoMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
std::ostream & console_out()
static int copyOutRow(T *, int, int, int, int *, int *, double *, int, int &)
bool Filled() const
static int setValues(T *mat, double scalar)
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
virtual const Epetra_Map & RowMatrixColMap() const =0
static int sumIn(T *, int, int, int, const int *, const int *, int, const double *)
static int copyOutRow(T *mat, int row, int len, double *coefs, int *indices)
int ReplaceMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
static int putValuesIn(T *mat, int numRows, const int *rows, int numCols, const int *cols, const double *const *values, bool sum_into)
static int getOffset(T *, int, int)
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const
static int globalAssemble(T *A)