FEI Version of the Day
fei_MatrixTraits_Epetra.hpp
00001 /*
00002 // @HEADER
00003 // ************************************************************************
00004 //             FEI: Finite Element Interface to Linear Solvers
00005 //                  Copyright (2005) Sandia Corporation.
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
00008 // U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Alan Williams (william@sandia.gov) 
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 */
00042 
00043 
00044 #ifndef _fei_MatrixTraits_Epetra_h_
00045 #define _fei_MatrixTraits_Epetra_h_
00046 
00047 #include <fei_trilinos_macros.hpp>
00048 
00049 #ifdef HAVE_FEI_EPETRA
00050 
00051 //
00052 //IMPORTANT NOTE: Make sure that wherever this file is included from, it
00053 //appears BEFORE any include of fei_Vector_Impl.hpp or fei_Matrix_Impl.hpp !!!
00054 //
00055 #include <fei_MatrixTraits.hpp>
00056 #include <snl_fei_BlockMatrixTraits.hpp>
00057 #include <fei_VectorTraits_Epetra.hpp>
00058 #include <fei_Include_Trilinos.hpp>
00059 #include <fei_Vector_Impl.hpp>
00060 
00061 namespace fei {
00068   template<>
00069   struct MatrixTraits<Epetra_CrsMatrix> {
00070     static const char* typeName()
00071       { return("Epetra_CrsMatrix"); }
00072 
00073     static double* getBeginPointer(Epetra_CrsMatrix* mat)
00074       {
00075         return (*mat)[0];
00076       }
00077 
00078     static int getOffset(Epetra_CrsMatrix* A, int row, int col)
00079       {
00080         const Epetra_Map& erowmap = A->RowMap();
00081         const Epetra_Map& ecolmap = A->ColMap();
00082         int local_row = erowmap.LID(row);
00083         int local_col = ecolmap.LID(col);
00084     
00085         int* rowOffsets;
00086         int* colIndices;
00087         double* coefs;
00088         A->ExtractCrsDataPointers(rowOffsets, colIndices, coefs);
00089     
00090         int* row_ptr = &colIndices[rowOffsets[local_row]];
00091         int* end_row = &colIndices[rowOffsets[local_row+1]];
00092     
00093         int col_offset = 0;
00094         for(; row_ptr != end_row; ++row_ptr) {
00095           if (*row_ptr == local_col) break;
00096           ++col_offset;
00097         }
00098     
00099         return rowOffsets[local_row] + col_offset;
00100       }
00101 
00102     static int setValues(Epetra_CrsMatrix* mat, double scalar)
00103       {
00104         return( mat->PutScalar(scalar) );
00105       }
00106 
00107     static int getNumLocalRows(Epetra_CrsMatrix* mat, int& numRows)
00108     {
00109       numRows = mat->NumMyRows();
00110       return(0);
00111     }
00112 
00113     static int getRowLength(Epetra_CrsMatrix* mat, int row, int& length)
00114       {
00115         length = mat->NumGlobalEntries(row);
00116         if (length < 0) return(-1);
00117         return( 0 );
00118       }
00119 
00120     static int copyOutRow(Epetra_CrsMatrix* mat,
00121                       int row, int len, double* coefs, int* indices)
00122       {
00123         int dummy;
00124         return(mat->ExtractGlobalRowCopy(row, len, dummy, coefs, indices));
00125       }
00126 
00127     static int putValuesIn(Epetra_CrsMatrix* mat,
00128                      int numRows, const int* rows,
00129                      int numCols, const int* cols,
00130                      const double* const* values,
00131                            bool sum_into)
00132       {
00134         static std::vector<int> idx;
00135         idx.resize(numRows+numCols);
00136         int* idx_row = &idx[0];
00137         int* idx_col = idx_row+numRows;
00138         for(int i=0; i<numRows; ++i) {
00139           idx_row[i] = mat->RowMap().LID(rows[i]);
00140         }
00141         for(int i=0; i<numCols; ++i) {
00142           idx_col[i] = mat->ColMap().LID(cols[i]);
00143         }
00144         if (sum_into) {
00145           for(int i=0; i<numRows; ++i) {
00146             int err = mat->SumIntoMyValues(idx_row[i], numCols,
00147                                                (double*)values[i],
00148                                                idx_col);
00149             if (err != 0) {
00150               return(err);
00151             }
00152           }
00153         }
00154         else {
00155           for(int i=0; i<numRows; ++i) {
00156             int err = mat->ReplaceMyValues(idx_row[i], numCols,
00157                                                (double*)values[i],
00158                                                idx_col);
00159             if (err != 0) {
00160               return(err);
00161             }
00162           }
00163         }
00164         return(0);
00165       }
00166 
00167     static int globalAssemble(Epetra_CrsMatrix* mat)
00168     {
00169       if (!mat->Filled()) {
00170         int err = mat->FillComplete();
00171         if (err != 0) {
00172           fei::console_out() << "MatrixTraits<Epetra_CrsMatrix>::globalAssemble"
00173                    << " ERROR in mat->FillComplete" << FEI_ENDL;
00174           return(-1);
00175         }
00176       }
00177 
00178       if (!mat->StorageOptimized()) {
00179         mat->OptimizeStorage();
00180       }
00181 
00182       return( 0 );
00183     }
00184 
00185     static int matvec(Epetra_CrsMatrix* mat,
00186                       fei::Vector* x,
00187                       fei::Vector* y)
00188     {
00189       fei::Vector_Impl<Epetra_MultiVector>* evx =
00190         dynamic_cast<fei::Vector_Impl<Epetra_MultiVector>* >(x);
00191       fei::Vector_Impl<Epetra_MultiVector>* evy =
00192         dynamic_cast<fei::Vector_Impl<Epetra_MultiVector>* >(y);
00193 
00194       if (evx == NULL || evy == NULL) {
00195         return(-1);
00196       }
00197 
00198       Epetra_MultiVector* ex = evx->getUnderlyingVector();
00199       Epetra_MultiVector* ey = evy->getUnderlyingVector();
00200 
00201       return( mat->Multiply(false, *ex, *ey) );
00202     }
00203 
00204   };//struct MatrixTraits<Epetra_CrsMatrix>
00205 }//namespace fei
00206 
00207 namespace snl_fei {
00214   template<>
00215   struct BlockMatrixTraits<Epetra_VbrMatrix> {
00216     static const char* typeName()
00217       { return("Epetra_VbrMatrix"); }
00218 
00219     static int putScalar(Epetra_VbrMatrix* mat, double scalar)
00220       {
00221         return( mat->PutScalar(scalar) );
00222       }
00223 
00224     static int getRowLength(Epetra_VbrMatrix* mat, int row, int& length)
00225       {
00226         length = mat->NumGlobalBlockEntries(row);
00227         return(0);
00228       }
00229 
00230     static int getPointRowLength(Epetra_VbrMatrix* mat, int row, int& length)
00231     {
00232       const Epetra_Map& map = mat->RowMatrixRowMap();
00233       int minLocalRow = map.MinMyGID();
00234       int localRow = row - minLocalRow;
00235       int error = mat->NumMyRowEntries(localRow, length);
00236       return(error);
00237     }
00238 
00239     static int copyOutRow(Epetra_VbrMatrix* mat,
00240                           int row, int numBlkCols,
00241                           int rowDim,
00242                           int* blkCols,
00243                           int* colDims,
00244                           double* coefs,
00245                           int coefsLen,
00246                           int& blkRowLength)
00247       {
00248         int checkRowDim;
00249         int error = mat->BeginExtractGlobalBlockRowCopy(row, numBlkCols,
00250                                                         checkRowDim,
00251                                                         blkRowLength,
00252                                                         blkCols, colDims);
00253         if (error != 0 || checkRowDim != rowDim || blkRowLength != numBlkCols) {
00254           return(error);
00255         }
00256 
00257         int offset = 0;
00258         for(int i=0; i<numBlkCols; ++i) {
00259           if (offset >= coefsLen) {
00260             cerr << "BlockMatrixTraits::copyOutRow ran off end of coefs array."
00261                  << endl;
00262             return(-2);
00263           }
00264           int numValues = rowDim*colDims[i];
00265           error = mat->ExtractEntryCopy(numValues, &(coefs[offset]),
00266                                         rowDim, false);
00267           if (error != 0) {
00268             return(error);
00269           }
00270           offset += numValues;
00271         }
00272 
00273         return(0);
00274       }
00275 
00276     static int copyOutPointRow(Epetra_VbrMatrix* mat,
00277                                int firstLocalOffset,
00278                                int row,
00279                                int len,
00280                                double* coefs,
00281                                int* indices,
00282                                int& rowLength)
00283       {
00284         int error = mat->ExtractMyRowCopy(row-firstLocalOffset,
00285                                           len, rowLength,
00286                                           coefs, indices);
00287 
00288         const Epetra_Map& colmap = mat->RowMatrixColMap();
00289         for(int i=0; i<len; ++i) {
00290           indices[i] = colmap.GID(indices[i]);
00291         }
00292 
00293         return(error);
00294       }
00295 
00296     static int sumIn(Epetra_VbrMatrix* mat,
00297                      int blockRow,
00298                      int rowDim,
00299                      int numBlockCols,
00300                      const int* blockCols,
00301                      const int* colDims,
00302                      int LDA,
00303                      const double* values)
00304     {
00305       int err, voffset = 0;
00306       for(int j=0; j<numBlockCols; ++j) {
00307         err = mat->DirectSubmitBlockEntry(blockRow, blockCols[j],
00308                                           &(values[voffset]), LDA,
00309                                          rowDim, colDims[j], true/*sum_into*/);
00310         if (err != 0) return(err);
00311 
00312         voffset += colDims[j]*LDA;
00313       }
00314 
00315       return(0);
00316     }
00317 
00318     static int copyIn(Epetra_VbrMatrix* mat,
00319                       int blockRow,
00320                       int rowDim,
00321                       int numBlockCols,
00322                       const int* blockCols,
00323                       const int* colDims,
00324                       int LDA,
00325                       const double* values)
00326     {
00327       int err, voffset = 0;
00328       for(int j=0; j<numBlockCols; ++j) {
00329         err = mat->DirectSubmitBlockEntry(blockRow, blockCols[j],
00330                                          &(values[voffset]), LDA,
00331                                     rowDim, colDims[j], false/*replace*/);
00332         if (err != 0) return(err);
00333 
00334         voffset += colDims[j]*LDA;
00335       }
00336 
00337       return(0);
00338     }
00339 
00340     static int sumIn(Epetra_VbrMatrix* mat,
00341                      int row,
00342                      int rowDim,
00343                      int numCols,
00344                      const int* cols,
00345                      const int* LDAs,
00346                      const int* colDims,
00347                      const double* const* values)
00348       {
00349         int err = 0;
00350         for(int i=0; i<numCols; ++i) {
00351           err = mat->DirectSubmitBlockEntry(row, cols[i], values[i],
00352                                             LDAs[i], rowDim, colDims[i], true);
00353           if (err != 0) return(err);
00354         }
00355 
00356         return(err);
00357       }
00358 
00359     static int copyIn(Epetra_VbrMatrix* mat,
00360                       int row,
00361                       int rowDim,
00362                       int numCols,
00363                       const int* cols,
00364                       const int* LDAs,
00365                       const int* colDims,
00366                       const double* const* values)
00367       {
00368         int err = 0;
00369         for(int i=0; i<numCols; ++i) {
00370           err = mat->DirectSubmitBlockEntry(row, cols[i], values[i],
00371                                           LDAs[i], rowDim, colDims[i], false);
00372           if (err != 0) return(err);
00373         }
00374 
00375         return(err);
00376       }
00377 
00378     static int globalAssemble(Epetra_VbrMatrix* mat)
00379     {
00380       return( mat->FillComplete() );
00381     }
00382   };//struct BlockMatrixTraits<Epetra_VbrMatrix>
00383 }//namespace snl_fei
00384 #endif //HAVE_FEI_EPETRA
00385 #endif // _fei_MatrixTraits_Epetra_hpp_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends