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 int setValues(Epetra_CrsMatrix* mat, double scalar)
00074       {
00075         return( mat->PutScalar(scalar) );
00076       }
00077 
00078     static int getNumLocalRows(Epetra_CrsMatrix* mat, int& numRows)
00079     {
00080       numRows = mat->NumMyRows();
00081       return(0);
00082     }
00083 
00084     static int getRowLength(Epetra_CrsMatrix* mat, int row, int& length)
00085       {
00086         length = mat->NumGlobalEntries(row);
00087         if (length < 0) return(-1);
00088         return( 0 );
00089       }
00090 
00091     static int copyOutRow(Epetra_CrsMatrix* mat,
00092                       int row, int len, double* coefs, int* indices)
00093       {
00094         int dummy;
00095         return(mat->ExtractGlobalRowCopy(row, len, dummy, coefs, indices));
00096       }
00097 
00098     static int putValuesIn(Epetra_CrsMatrix* mat,
00099                      int numRows, const int* rows,
00100                      int numCols, const int* cols,
00101                      const double* const* values,
00102                            bool sum_into)
00103       {
00105         static std::vector<int> idx;
00106         idx.resize(numRows+numCols);
00107         int* idx_row = &idx[0];
00108         int* idx_col = idx_row+numRows;
00109         for(int i=0; i<numRows; ++i) {
00110           idx_row[i] = mat->RowMap().LID(rows[i]);
00111         }
00112         for(int i=0; i<numCols; ++i) {
00113           idx_col[i] = mat->ColMap().LID(cols[i]);
00114         }
00115         if (sum_into) {
00116           for(int i=0; i<numRows; ++i) {
00117             int err = mat->SumIntoMyValues(idx_row[i], numCols,
00118                                                (double*)values[i],
00119                                                idx_col);
00120             if (err != 0) {
00121               return(err);
00122             }
00123           }
00124         }
00125         else {
00126           for(int i=0; i<numRows; ++i) {
00127             int err = mat->ReplaceMyValues(idx_row[i], numCols,
00128                                                (double*)values[i],
00129                                                idx_col);
00130             if (err != 0) {
00131               return(err);
00132             }
00133           }
00134         }
00135         return(0);
00136       }
00137 
00138     static int globalAssemble(Epetra_CrsMatrix* mat)
00139     {
00140       if (!mat->Filled()) {
00141         int err = mat->FillComplete();
00142         if (err != 0) {
00143           fei::console_out() << "MatrixTraits<Epetra_CrsMatrix>::globalAssemble"
00144                    << " ERROR in mat->FillComplete" << FEI_ENDL;
00145           return(-1);
00146         }
00147       }
00148 
00149       if (!mat->StorageOptimized()) {
00150         mat->OptimizeStorage();
00151       }
00152 
00153       return( 0 );
00154     }
00155 
00156     static int matvec(Epetra_CrsMatrix* mat,
00157                       fei::Vector* x,
00158                       fei::Vector* y)
00159     {
00160       fei::Vector_Impl<Epetra_MultiVector>* evx =
00161         dynamic_cast<fei::Vector_Impl<Epetra_MultiVector>* >(x);
00162       fei::Vector_Impl<Epetra_MultiVector>* evy =
00163         dynamic_cast<fei::Vector_Impl<Epetra_MultiVector>* >(y);
00164 
00165       if (evx == NULL || evy == NULL) {
00166         return(-1);
00167       }
00168 
00169       Epetra_MultiVector* ex = evx->getUnderlyingVector();
00170       Epetra_MultiVector* ey = evy->getUnderlyingVector();
00171 
00172       return( mat->Multiply(false, *ex, *ey) );
00173     }
00174 
00175   };//struct MatrixTraits<Epetra_CrsMatrix>
00176 }//namespace fei
00177 
00178 namespace snl_fei {
00185   template<>
00186   struct BlockMatrixTraits<Epetra_VbrMatrix> {
00187     static const char* typeName()
00188       { return("Epetra_VbrMatrix"); }
00189 
00190     static int putScalar(Epetra_VbrMatrix* mat, double scalar)
00191       {
00192         return( mat->PutScalar(scalar) );
00193       }
00194 
00195     static int getRowLength(Epetra_VbrMatrix* mat, int row, int& length)
00196       {
00197         length = mat->NumGlobalBlockEntries(row);
00198         return(0);
00199       }
00200 
00201     static int getPointRowLength(Epetra_VbrMatrix* mat, int row, int& length)
00202     {
00203       const Epetra_Map& map = mat->RowMatrixRowMap();
00204       int minLocalRow = map.MinMyGID();
00205       int localRow = row - minLocalRow;
00206       int error = mat->NumMyRowEntries(localRow, length);
00207       return(error);
00208     }
00209 
00210     static int copyOutRow(Epetra_VbrMatrix* mat,
00211                           int row, int numBlkCols,
00212                           int rowDim,
00213                           int* blkCols,
00214                           int* colDims,
00215                           double* coefs,
00216                           int coefsLen,
00217                           int& blkRowLength)
00218       {
00219         int checkRowDim;
00220         int error = mat->BeginExtractGlobalBlockRowCopy(row, numBlkCols,
00221                                                         checkRowDim,
00222                                                         blkRowLength,
00223                                                         blkCols, colDims);
00224         if (error != 0 || checkRowDim != rowDim || blkRowLength != numBlkCols) {
00225           return(error);
00226         }
00227 
00228         int offset = 0;
00229         for(int i=0; i<numBlkCols; ++i) {
00230           if (offset >= coefsLen) {
00231             cerr << "BlockMatrixTraits::copyOutRow ran off end of coefs array."
00232                  << endl;
00233             return(-2);
00234           }
00235           int numValues = rowDim*colDims[i];
00236           error = mat->ExtractEntryCopy(numValues, &(coefs[offset]),
00237                                         rowDim, false);
00238           if (error != 0) {
00239             return(error);
00240           }
00241           offset += numValues;
00242         }
00243 
00244         return(0);
00245       }
00246 
00247     static int copyOutPointRow(Epetra_VbrMatrix* mat,
00248                                int firstLocalOffset,
00249                                int row,
00250                                int len,
00251                                double* coefs,
00252                                int* indices,
00253                                int& rowLength)
00254       {
00255         int error = mat->ExtractMyRowCopy(row-firstLocalOffset,
00256                                           len, rowLength,
00257                                           coefs, indices);
00258 
00259         const Epetra_Map& colmap = mat->RowMatrixColMap();
00260         for(int i=0; i<len; ++i) {
00261           indices[i] = colmap.GID(indices[i]);
00262         }
00263 
00264         return(error);
00265       }
00266 
00267     static int sumIn(Epetra_VbrMatrix* mat,
00268                      int blockRow,
00269                      int rowDim,
00270                      int numBlockCols,
00271                      const int* blockCols,
00272                      const int* colDims,
00273                      int LDA,
00274                      const double* values)
00275     {
00276       int err, voffset = 0;
00277       for(int j=0; j<numBlockCols; ++j) {
00278         err = mat->DirectSubmitBlockEntry(blockRow, blockCols[j],
00279                                           &(values[voffset]), LDA,
00280                                          rowDim, colDims[j], true/*sum_into*/);
00281         if (err != 0) return(err);
00282 
00283         voffset += colDims[j]*LDA;
00284       }
00285 
00286       return(0);
00287     }
00288 
00289     static int copyIn(Epetra_VbrMatrix* mat,
00290                       int blockRow,
00291                       int rowDim,
00292                       int numBlockCols,
00293                       const int* blockCols,
00294                       const int* colDims,
00295                       int LDA,
00296                       const double* values)
00297     {
00298       int err, voffset = 0;
00299       for(int j=0; j<numBlockCols; ++j) {
00300         err = mat->DirectSubmitBlockEntry(blockRow, blockCols[j],
00301                                          &(values[voffset]), LDA,
00302                                     rowDim, colDims[j], false/*replace*/);
00303         if (err != 0) return(err);
00304 
00305         voffset += colDims[j]*LDA;
00306       }
00307 
00308       return(0);
00309     }
00310 
00311     static int sumIn(Epetra_VbrMatrix* mat,
00312                      int row,
00313                      int rowDim,
00314                      int numCols,
00315                      const int* cols,
00316                      const int* LDAs,
00317                      const int* colDims,
00318                      const double* const* values)
00319       {
00320         int err = 0;
00321         for(int i=0; i<numCols; ++i) {
00322           err = mat->DirectSubmitBlockEntry(row, cols[i], values[i],
00323                                             LDAs[i], rowDim, colDims[i], true);
00324           if (err != 0) return(err);
00325         }
00326 
00327         return(err);
00328       }
00329 
00330     static int copyIn(Epetra_VbrMatrix* mat,
00331                       int row,
00332                       int rowDim,
00333                       int numCols,
00334                       const int* cols,
00335                       const int* LDAs,
00336                       const int* colDims,
00337                       const double* const* values)
00338       {
00339         int err = 0;
00340         for(int i=0; i<numCols; ++i) {
00341           err = mat->DirectSubmitBlockEntry(row, cols[i], values[i],
00342                                           LDAs[i], rowDim, colDims[i], false);
00343           if (err != 0) return(err);
00344         }
00345 
00346         return(err);
00347       }
00348 
00349     static int globalAssemble(Epetra_VbrMatrix* mat)
00350     {
00351       return( mat->FillComplete() );
00352     }
00353   };//struct BlockMatrixTraits<Epetra_VbrMatrix>
00354 }//namespace snl_fei
00355 #endif //HAVE_FEI_EPETRA
00356 #endif // _fei_MatrixTraits_Epetra_hpp_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends