FEI Version of the Day
fei_MatrixTraits_Epetra.hpp
00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2005 Sandia Corporation.                              */
00003 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00004 /*    non-exclusive license for use of this work by or on behalf      */
00005 /*    of the U.S. Government.  Export of this program may require     */
00006 /*    a license from the United States Government.                    */
00007 /*--------------------------------------------------------------------*/
00008 
00009 #ifndef _fei_MatrixTraits_Epetra_h_
00010 #define _fei_MatrixTraits_Epetra_h_
00011 
00012 #include <fei_trilinos_macros.hpp>
00013 
00014 #ifdef HAVE_FEI_EPETRA
00015 
00016 //
00017 //IMPORTANT NOTE: Make sure that wherever this file is included from, it
00018 //appears BEFORE any include of fei_Vector_Impl.hpp or fei_Matrix_Impl.hpp !!!
00019 //
00020 #include <fei_MatrixTraits.hpp>
00021 #include <snl_fei_BlockMatrixTraits.hpp>
00022 #include <fei_VectorTraits_Epetra.hpp>
00023 #include <fei_Include_Trilinos.hpp>
00024 #include <fei_Vector_Impl.hpp>
00025 
00026 namespace fei {
00033   template<>
00034   struct MatrixTraits<Epetra_CrsMatrix> {
00035     static const char* typeName()
00036       { return("Epetra_CrsMatrix"); }
00037 
00038     static int setValues(Epetra_CrsMatrix* mat, double scalar)
00039       {
00040         return( mat->PutScalar(scalar) );
00041       }
00042 
00043     static int getNumLocalRows(Epetra_CrsMatrix* mat, int& numRows)
00044     {
00045       numRows = mat->NumMyRows();
00046       return(0);
00047     }
00048 
00049     static int getRowLength(Epetra_CrsMatrix* mat, int row, int& length)
00050       {
00051         length = mat->NumGlobalEntries(row);
00052         if (length < 0) return(-1);
00053         return( 0 );
00054       }
00055 
00056     static int copyOutRow(Epetra_CrsMatrix* mat,
00057                       int row, int len, double* coefs, int* indices)
00058       {
00059         int dummy;
00060         return(mat->ExtractGlobalRowCopy(row, len, dummy, coefs, indices));
00061       }
00062 
00063     static int putValuesIn(Epetra_CrsMatrix* mat,
00064                      int numRows, const int* rows,
00065                      int numCols, const int* cols,
00066                      const double* const* values,
00067                            bool sum_into)
00068       {
00070         static std::vector<int> idx;
00071         idx.resize(numRows+numCols);
00072         int* idx_row = &idx[0];
00073         int* idx_col = idx_row+numRows;
00074         for(int i=0; i<numRows; ++i) {
00075           idx_row[i] = mat->RowMap().LID(rows[i]);
00076         }
00077         for(int i=0; i<numCols; ++i) {
00078           idx_col[i] = mat->ColMap().LID(cols[i]);
00079         }
00080         if (sum_into) {
00081           for(int i=0; i<numRows; ++i) {
00082             int err = mat->SumIntoMyValues(idx_row[i], numCols,
00083                                                (double*)values[i],
00084                                                idx_col);
00085             if (err != 0) {
00086               return(err);
00087             }
00088           }
00089         }
00090         else {
00091           for(int i=0; i<numRows; ++i) {
00092             int err = mat->ReplaceMyValues(idx_row[i], numCols,
00093                                                (double*)values[i],
00094                                                idx_col);
00095             if (err != 0) {
00096               return(err);
00097             }
00098           }
00099         }
00100         return(0);
00101       }
00102 
00103     static int globalAssemble(Epetra_CrsMatrix* mat)
00104     {
00105       if (!mat->Filled()) {
00106         int err = mat->FillComplete();
00107         if (err != 0) {
00108           fei::console_out() << "MatrixTraits<Epetra_CrsMatrix>::globalAssemble"
00109                    << " ERROR in mat->FillComplete" << FEI_ENDL;
00110           return(-1);
00111         }
00112       }
00113 
00114       if (!mat->StorageOptimized()) {
00115         mat->OptimizeStorage();
00116       }
00117 
00118       return( 0 );
00119     }
00120 
00121     static int matvec(Epetra_CrsMatrix* mat,
00122                       fei::Vector* x,
00123                       fei::Vector* y)
00124     {
00125       fei::Vector_Impl<Epetra_MultiVector>* evx =
00126         dynamic_cast<fei::Vector_Impl<Epetra_MultiVector>* >(x);
00127       fei::Vector_Impl<Epetra_MultiVector>* evy =
00128         dynamic_cast<fei::Vector_Impl<Epetra_MultiVector>* >(y);
00129 
00130       if (evx == NULL || evy == NULL) {
00131         return(-1);
00132       }
00133 
00134       Epetra_MultiVector* ex = evx->getUnderlyingVector();
00135       Epetra_MultiVector* ey = evy->getUnderlyingVector();
00136 
00137       return( mat->Multiply(false, *ex, *ey) );
00138     }
00139 
00140   };//struct MatrixTraits<Epetra_CrsMatrix>
00141 }//namespace fei
00142 
00143 namespace snl_fei {
00150   template<>
00151   struct BlockMatrixTraits<Epetra_VbrMatrix> {
00152     static const char* typeName()
00153       { return("Epetra_VbrMatrix"); }
00154 
00155     static int putScalar(Epetra_VbrMatrix* mat, double scalar)
00156       {
00157         return( mat->PutScalar(scalar) );
00158       }
00159 
00160     static int getRowLength(Epetra_VbrMatrix* mat, int row, int& length)
00161       {
00162         length = mat->NumGlobalBlockEntries(row);
00163         return(0);
00164       }
00165 
00166     static int getPointRowLength(Epetra_VbrMatrix* mat, int row, int& length)
00167     {
00168       const Epetra_Map& map = mat->RowMatrixRowMap();
00169       int minLocalRow = map.MinMyGID();
00170       int localRow = row - minLocalRow;
00171       int error = mat->NumMyRowEntries(localRow, length);
00172       return(error);
00173     }
00174 
00175     static int copyOutRow(Epetra_VbrMatrix* mat,
00176                           int row, int numBlkCols,
00177                           int rowDim,
00178                           int* blkCols,
00179                           int* colDims,
00180                           double* coefs,
00181                           int coefsLen,
00182                           int& blkRowLength)
00183       {
00184         int checkRowDim;
00185         int error = mat->BeginExtractGlobalBlockRowCopy(row, numBlkCols,
00186                                                         checkRowDim,
00187                                                         blkRowLength,
00188                                                         blkCols, colDims);
00189         if (error != 0 || checkRowDim != rowDim || blkRowLength != numBlkCols) {
00190           return(error);
00191         }
00192 
00193         int offset = 0;
00194         for(int i=0; i<numBlkCols; ++i) {
00195           if (offset >= coefsLen) {
00196             cerr << "BlockMatrixTraits::copyOutRow ran off end of coefs array."
00197                  << endl;
00198             return(-2);
00199           }
00200           int numValues = rowDim*colDims[i];
00201           error = mat->ExtractEntryCopy(numValues, &(coefs[offset]),
00202                                         rowDim, false);
00203           if (error != 0) {
00204             return(error);
00205           }
00206           offset += numValues;
00207         }
00208 
00209         return(0);
00210       }
00211 
00212     static int copyOutPointRow(Epetra_VbrMatrix* mat,
00213                                int firstLocalOffset,
00214                                int row,
00215                                int len,
00216                                double* coefs,
00217                                int* indices,
00218                                int& rowLength)
00219       {
00220         int error = mat->ExtractMyRowCopy(row-firstLocalOffset,
00221                                           len, rowLength,
00222                                           coefs, indices);
00223 
00224         const Epetra_Map& colmap = mat->RowMatrixColMap();
00225         for(int i=0; i<len; ++i) {
00226           indices[i] = colmap.GID(indices[i]);
00227         }
00228 
00229         return(error);
00230       }
00231 
00232     static int sumIn(Epetra_VbrMatrix* mat,
00233                      int blockRow,
00234                      int rowDim,
00235                      int numBlockCols,
00236                      const int* blockCols,
00237                      const int* colDims,
00238                      int LDA,
00239                      const double* values)
00240     {
00241       int err, voffset = 0;
00242       for(int j=0; j<numBlockCols; ++j) {
00243         err = mat->DirectSubmitBlockEntry(blockRow, blockCols[j],
00244                                           &(values[voffset]), LDA,
00245                                          rowDim, colDims[j], true/*sum_into*/);
00246         if (err != 0) return(err);
00247 
00248         voffset += colDims[j]*LDA;
00249       }
00250 
00251       return(0);
00252     }
00253 
00254     static int copyIn(Epetra_VbrMatrix* mat,
00255                       int blockRow,
00256                       int rowDim,
00257                       int numBlockCols,
00258                       const int* blockCols,
00259                       const int* colDims,
00260                       int LDA,
00261                       const double* values)
00262     {
00263       int err, voffset = 0;
00264       for(int j=0; j<numBlockCols; ++j) {
00265         err = mat->DirectSubmitBlockEntry(blockRow, blockCols[j],
00266                                          &(values[voffset]), LDA,
00267                                     rowDim, colDims[j], false/*replace*/);
00268         if (err != 0) return(err);
00269 
00270         voffset += colDims[j]*LDA;
00271       }
00272 
00273       return(0);
00274     }
00275 
00276     static int sumIn(Epetra_VbrMatrix* mat,
00277                      int row,
00278                      int rowDim,
00279                      int numCols,
00280                      const int* cols,
00281                      const int* LDAs,
00282                      const int* colDims,
00283                      const double* const* values)
00284       {
00285         int err = 0;
00286         for(int i=0; i<numCols; ++i) {
00287           err = mat->DirectSubmitBlockEntry(row, cols[i], values[i],
00288                                             LDAs[i], rowDim, colDims[i], true);
00289           if (err != 0) return(err);
00290         }
00291 
00292         return(err);
00293       }
00294 
00295     static int copyIn(Epetra_VbrMatrix* mat,
00296                       int row,
00297                       int rowDim,
00298                       int numCols,
00299                       const int* cols,
00300                       const int* LDAs,
00301                       const int* colDims,
00302                       const double* const* values)
00303       {
00304         int err = 0;
00305         for(int i=0; i<numCols; ++i) {
00306           err = mat->DirectSubmitBlockEntry(row, cols[i], values[i],
00307                                           LDAs[i], rowDim, colDims[i], false);
00308           if (err != 0) return(err);
00309         }
00310 
00311         return(err);
00312       }
00313 
00314     static int globalAssemble(Epetra_VbrMatrix* mat)
00315     {
00316       return( mat->FillComplete() );
00317     }
00318   };//struct BlockMatrixTraits<Epetra_VbrMatrix>
00319 }//namespace snl_fei
00320 #endif //HAVE_FEI_EPETRA
00321 #endif // _fei_MatrixTraits_Epetra_hpp_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends