Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_OverlappingRowMatrix_def.hpp
00001 /*@HEADER
00002 // ***********************************************************************
00003 // 
00004 //       Ifpack2: Tempated Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2009) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #ifndef IFPACK2_OVERLAPPINGROWMATRIX_DEF_HPP
00031 #define IFPACK2_OVERLAPPINGROWMATRIX_DEF_HPP
00032 #include "Ifpack2_OverlappingRowMatrix_decl.hpp"
00033 #include "Tpetra_CrsMatrix.hpp"
00034 
00035 #include "Teuchos_CommHelpers.hpp"
00036 
00037 namespace Ifpack2 {
00038 //==========================================================================
00039 // Standard constructor
00040 template<class MatrixType>
00041 OverlappingRowMatrix<MatrixType>::OverlappingRowMatrix(const Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& Matrix_in,
00042                    int OverlapLevel_in) :
00043   A_(Matrix_in),
00044   OverlapLevel_(OverlapLevel_in),
00045   UseSubComm_(false)
00046 {  
00047   using Teuchos::RCP;
00048   using Teuchos::rcp;
00049   using Teuchos::Array;
00050 
00051   // Short form names
00052   typedef typename Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node>        MapType;
00053   typedef typename Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node>           ImportType;
00054   typedef typename Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> CrsMatrixType;
00055   typedef typename Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> RowMatrixType;
00056 
00057   TEUCHOS_TEST_FOR_EXCEPTION(OverlapLevel_ <= 0, std::runtime_error, "Ifpack2::OverlappingRowMatrix: OverlapLevel must be > 0.");
00058   TEUCHOS_TEST_FOR_EXCEPTION(A_->getComm()->getSize() == 1, std::runtime_error, "Ifpack2::OverlappingRowMatrix: Matrix must be parallel.");
00059 
00060   RCP<const CrsMatrixType> ACRS = Teuchos::rcp_dynamic_cast<const CrsMatrixType,const RowMatrixType>(A_);
00061   TEUCHOS_TEST_FOR_EXCEPTION(ACRS == Teuchos::null, std::runtime_error, "Ifpack2::OverlappingRowMatrix: Matrix must be Tpetra::CrsMatrix.");
00062   
00063   NumMyRowsA_ = A_->getNodeNumRows();
00064 
00065   GlobalOrdinal global_invalid = Teuchos::OrdinalTraits<global_size_t>::invalid();
00066 
00067   // Temp arrays
00068   Array<GlobalOrdinal> ExtElements;  
00069   RCP<MapType>         TmpMap;
00070   RCP<CrsMatrixType>   TmpMatrix; 
00071   RCP<ImportType>      TmpImporter;
00072   RCP<const MapType>   RowMap, ColMap;
00073 
00074   // The big import loop
00075   for (int overlap = 0 ; overlap < OverlapLevel_ ; ++overlap) {
00076     // Get the current maps
00077     if(overlap==0){
00078       RowMap = A_->getRowMap();
00079       ColMap = A_->getColMap(); 
00080     }
00081     else {
00082       RowMap = TmpMatrix->getRowMap();
00083       ColMap = TmpMatrix->getColMap(); 
00084     }
00085 
00086     size_t size = ColMap->getNodeNumElements() - RowMap->getNodeNumElements(); 
00087     Array<GlobalOrdinal> mylist(size); 
00088     size_t count = 0; 
00089 
00090     // define the set of rows that are in ColMap but not in RowMap 
00091     for (LocalOrdinal i = 0 ; (size_t) i < ColMap->getNodeNumElements() ; ++i) { 
00092       GlobalOrdinal GID = ColMap->getGlobalElement(i); 
00093       if (A_->getRowMap()->getLocalElement(GID) ==  global_invalid) { 
00094   typename Array<GlobalOrdinal>::iterator pos 
00095     = std::find(ExtElements.begin(),ExtElements.end(),GID); 
00096         if (pos == ExtElements.end()) { 
00097           ExtElements.push_back(GID);
00098           mylist[count] = GID; 
00099           ++count; 
00100         } 
00101       } 
00102     }
00103     
00104     // Allocate & import new matrices, maps, etc.
00105     TmpMap      = rcp(new MapType(global_invalid,mylist(0,count),Teuchos::OrdinalTraits<GlobalOrdinal>::zero(),A_->getComm(),A_->getNode()));
00106     TmpMatrix   = rcp(new CrsMatrixType(TmpMap,0));
00107     TmpImporter = rcp(new ImportType(A_->getRowMap(),TmpMap));
00108 
00109     TmpMatrix->doImport(*ACRS,*TmpImporter,Tpetra::INSERT);
00110     TmpMatrix->fillComplete(A_->getDomainMap(),TmpMap);
00111   }
00112 
00113   // build the map containing all the nodes (original
00114   // matrix + extended matrix)
00115   Array<GlobalOrdinal> mylist(NumMyRowsA_ + ExtElements.size());
00116   for (LocalOrdinal i = 0 ; (size_t)i < NumMyRowsA_ ; ++i)
00117     mylist[i] = A_->getRowMap()->getGlobalElement(i);
00118   for (LocalOrdinal i = 0 ; i < ExtElements.size() ; ++i)
00119     mylist[i + NumMyRowsA_] = ExtElements[i];
00120 
00121   RowMap_= rcp(new MapType(global_invalid,mylist(),Teuchos::OrdinalTraits<GlobalOrdinal>::zero(),A_->getComm(),A_->getNode()));
00122   ColMap_= RowMap_;
00123 
00124   // now build the map corresponding to all the external nodes
00125   // (with respect to A().RowMatrixRowMap().
00126   ExtMap_      = rcp(new MapType(global_invalid,ExtElements(),Teuchos::OrdinalTraits<GlobalOrdinal>::zero(),A_->getComm(),A_->getNode()));
00127   ExtMatrix_   = rcp(new CrsMatrixType(ExtMap_,ColMap_,0));
00128   ExtImporter_ = rcp(new ImportType(A_->getRowMap(),ExtMap_));
00129 
00130   RCP<CrsMatrixType> ExtMatrixCRS = Teuchos::rcp_dynamic_cast<CrsMatrixType,RowMatrixType>(ExtMatrix_);
00131   ExtMatrixCRS->doImport(*ACRS,*ExtImporter_,Tpetra::INSERT);
00132   ExtMatrixCRS->fillComplete(A_->getDomainMap(),RowMap_);
00133 
00134   Importer_ = rcp(new ImportType(A_->getRowMap(),RowMap_));
00135 
00136   // fix indices for overlapping matrix
00137   NumMyRowsB_ = ExtMatrix_->getNodeNumRows();
00138   NumMyRows_ = NumMyRowsA_ + NumMyRowsB_;
00139   NumMyCols_ = NumMyRows_;
00140 
00141   NumMyDiagonals_ = A_->getNodeNumDiags() + ExtMatrix_->getNodeNumDiags(); 
00142   NumMyNonzeros_  = A_->getNodeNumEntries() + ExtMatrix_->getNodeNumEntries();
00143 
00144   // FIXME: Fix this later when Teuchos::Comm gets redone
00145   Tpetra::global_size_t NumMyNonzeros_tmp = NumMyNonzeros_;
00146   Teuchos::reduceAll<int,Tpetra::global_size_t>(*A_->getComm(),Teuchos::REDUCE_SUM,NumMyNonzeros_tmp,Teuchos::outArg(NumGlobalNonzeros_));  
00147   Tpetra::global_size_t NumMyRows_tmp = NumMyRows_;
00148   Teuchos::reduceAll<int,Tpetra::global_size_t>(*A_->getComm(),Teuchos::REDUCE_SUM,NumMyRows_tmp,Teuchos::outArg(NumGlobalRows_));  
00149 
00150   MaxNumEntries_ = A_->getNodeMaxNumRowEntries();  
00151   if (MaxNumEntries_ < ExtMatrix_->getNodeMaxNumRowEntries())
00152     MaxNumEntries_ = ExtMatrix_->getNodeMaxNumRowEntries();
00153 
00154   // Resize temp arrays
00155   Indices_.resize(MaxNumEntries_);
00156   Values_.resize(MaxNumEntries_);
00157 }
00158 
00159 
00160 
00161 //==========================================================================
00162 // Sub-communicator constructor
00163 template<class MatrixType>
00164 OverlappingRowMatrix<MatrixType>::OverlappingRowMatrix(const Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& Matrix_in,
00165                     int OverlapLevel_in, int subdomainID)
00166 {
00167   //FIXME
00168   throw std::runtime_error("Ifpack2::OverlappingRowMatrix - subdomain code not implemented yet.");
00169 }
00170 
00171 //==========================================================================
00172 // Destructor
00173 template<class MatrixType>
00174 OverlappingRowMatrix<MatrixType>::~OverlappingRowMatrix()
00175 {
00176 
00177 }
00178 //==========================================================================
00179 // Returns the communicator.
00180 template<class MatrixType>
00181 const Teuchos::RCP<const Teuchos::Comm<int> > & OverlappingRowMatrix<MatrixType>::getComm() const
00182 {
00183   return A_->getComm();
00184 }
00185   
00186 //==========================================================================
00187 // Returns the underlying node.
00188 template<class MatrixType>
00189 Teuchos::RCP<typename MatrixType::node_type> OverlappingRowMatrix<MatrixType>::getNode() const
00190 {
00191   return A_->getNode();
00192 }
00193   
00194 //==========================================================================
00195 // Returns the Map that describes the row distribution in this matrix.
00196 template<class MatrixType>
00197 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> > & OverlappingRowMatrix<MatrixType>::getRowMap() const
00198 {
00199   return RowMap_;
00200 }
00201   
00202 //==========================================================================
00203 //  Returns the Map that describes the column distribution in this matrix.
00204 template<class MatrixType>
00205 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> > & OverlappingRowMatrix<MatrixType>::getColMap() const
00206 {
00207   return ColMap_;
00208 }
00209 
00210 //==========================================================================
00211 //  Returns the Map that describes the column distribution in this matrix.
00212 template<class MatrixType>
00213 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> > & OverlappingRowMatrix<MatrixType>::getDomainMap() const
00214 {
00215   return A_->getDomainMap();
00216 }
00217 
00218 //==========================================================================
00219 //  Returns the Map that describes the column distribution in this matrix.
00220 template<class MatrixType>
00221 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> > & OverlappingRowMatrix<MatrixType>::getRangeMap() const
00222 {
00223   return A_->getRangeMap();
00224 }
00225   
00226 //==========================================================================
00227   // Returns the RowGraph associated with this matrix. 
00228 template<class MatrixType>
00229 Teuchos::RCP<const Tpetra::RowGraph<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> > OverlappingRowMatrix<MatrixType>::getGraph() const
00230 {
00231   throw std::runtime_error("Ifpack2::OverlappingRowMatrix - ERROR getGraph() is not implemented.");
00232 }
00233   
00234 //==========================================================================
00235 // Returns the number of global rows in this matrix.
00236 template<class MatrixType>
00237 global_size_t OverlappingRowMatrix<MatrixType>::getGlobalNumRows() const
00238 {
00239   return NumGlobalRows_;
00240 }
00241   
00242 //==========================================================================
00243 //  Returns the number of global columns in this matrix.
00244 template<class MatrixType>
00245 global_size_t OverlappingRowMatrix<MatrixType>::getGlobalNumCols() const
00246 {
00247   return NumGlobalRows_;
00248 }
00249   
00250 //==========================================================================
00251 // Returns the number of rows owned on the calling node.
00252 template<class MatrixType>
00253 size_t OverlappingRowMatrix<MatrixType>::getNodeNumRows() const
00254 {
00255   return NumMyRows_;
00256 }
00257   
00258 //==========================================================================
00259 // Returns the number of columns needed to apply the forward operator on this node, i.e., the number of elements listed in the column map.
00260 template<class MatrixType>
00261 size_t OverlappingRowMatrix<MatrixType>::getNodeNumCols() const
00262 {
00263   return NumMyCols_;
00264 }
00265   
00266 //==========================================================================
00267 // Returns the index base for global indices for this matrix. 
00268 template<class MatrixType>
00269 typename MatrixType::global_ordinal_type OverlappingRowMatrix<MatrixType>::getIndexBase() const
00270 {
00271   return A_->getIndexBase();
00272 }
00273   
00274 //==========================================================================
00275 // Returns the global number of entries in this matrix.
00276 template<class MatrixType>
00277 Tpetra::global_size_t OverlappingRowMatrix<MatrixType>::getGlobalNumEntries() const
00278 {
00279   return NumGlobalNonzeros_;
00280 }
00281   
00282 //==========================================================================
00283 // Returns the local number of entries in this matrix.
00284 template<class MatrixType>
00285 size_t OverlappingRowMatrix<MatrixType>::getNodeNumEntries() const
00286 {
00287   return NumMyNonzeros_;
00288 }
00289   
00290 //==========================================================================
00291 //  Returns the current number of entries on this node in the specified global row.
00292 /* Returns Teuchos::OrdinalTraits<size_t>::invalid() if the specified global row is not valid for this graph. */
00293 template<class MatrixType>
00294 size_t OverlappingRowMatrix<MatrixType>::getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const
00295 {
00296   LocalOrdinal localRow = RowMap_->getLocalElement(globalRow);
00297   if(localRow == Teuchos::OrdinalTraits<LocalOrdinal>::invalid()) return Teuchos::OrdinalTraits<size_t>::invalid();
00298   else return getNumEntriesInLocalRow(localRow);
00299 }
00300 
00301   
00302 //==========================================================================
00303 // Returns the current number of entries on this node in the specified local row.
00304 /* Returns Teuchos::OrdinalTraits<size_t>::invalid() if the specified local row is not valid for this graph. */
00305 template<class MatrixType>
00306 size_t OverlappingRowMatrix<MatrixType>::getNumEntriesInLocalRow(LocalOrdinal localRow) const
00307 {
00308   if ((size_t)localRow < NumMyRowsA_)
00309     return A_->getNumEntriesInLocalRow(localRow);
00310   else
00311     return ExtMatrix_->getNumEntriesInLocalRow((LocalOrdinal)(localRow-NumMyRowsA_));
00312 }
00313   
00314 //==========================================================================
00315 //  Returns the number of global diagonal entries, based on global row/column index comparisons. 
00316 template<class MatrixType>
00317 global_size_t OverlappingRowMatrix<MatrixType>::getGlobalNumDiags() const
00318 {
00319   throw std::runtime_error("Ifpack2::OverlappingRowMatrix::getGlobalNumDiags() not supported.");
00320 }
00321   
00322 //==========================================================================
00323 //  Returns the number of local diagonal entries, based on global row/column index comparisons. 
00324 template<class MatrixType>
00325 size_t OverlappingRowMatrix<MatrixType>::getNodeNumDiags() const
00326 {
00327   return A_->getNodeNumDiags();
00328 }
00329   
00330 //==========================================================================
00331 //  Returns the maximum number of entries across all rows/columns on all nodes.
00332 template<class MatrixType>
00333 size_t OverlappingRowMatrix<MatrixType>::getGlobalMaxNumRowEntries() const
00334 {
00335   throw std::runtime_error("Ifpack2::OverlappingRowMatrix::getGlobalMaxNumRowEntries() not supported.");
00336 }
00337   
00338 //==========================================================================
00339 template<class MatrixType>
00340 //  Returns the maximum number of entries across all rows/columns on this node.
00341 size_t OverlappingRowMatrix<MatrixType>::getNodeMaxNumRowEntries() const
00342 {
00343   return MaxNumEntries_;
00344 }
00345   
00346 //==========================================================================
00347 //  Indicates whether this matrix has a well-defined column map
00348 template<class MatrixType>
00349 bool OverlappingRowMatrix<MatrixType>::hasColMap() const
00350 {
00351   return true;
00352 }
00353   
00354 //==========================================================================
00355 //  Indicates whether this matrix is lower triangular.
00356 template<class MatrixType>
00357 bool OverlappingRowMatrix<MatrixType>::isLowerTriangular() const
00358 {
00359   return A_->isLowerTriangular();
00360 }
00361   
00362 //==========================================================================
00363 //  Indicates whether this matrix is upper triangular.
00364 template<class MatrixType>
00365 bool OverlappingRowMatrix<MatrixType>::isUpperTriangular() const
00366 {
00367   return A_->isUpperTriangular();
00368 } 
00369 //==========================================================================
00370 //  If matrix indices are in the local range, this function returns true. Otherwise, this function returns false. */
00371 template<class MatrixType>
00372 bool OverlappingRowMatrix<MatrixType>::isLocallyIndexed() const
00373 {
00374   return true;
00375 }
00376    
00377 //==========================================================================
00378 //  If matrix indices are in the global range, this function returns true. Otherwise, this function returns false. */
00379 template<class MatrixType>
00380 bool OverlappingRowMatrix<MatrixType>::isGloballyIndexed() const
00381 {
00382   return false;
00383 }
00384   
00385 //==========================================================================
00386 // Returns \c true if fillComplete() has been called.
00387 template<class MatrixType>
00388 bool OverlappingRowMatrix<MatrixType>::isFillComplete() const
00389 {
00390   return true;
00391 }
00392   
00393 //==========================================================================
00394 // Extract a list of entries in a specified global row of this matrix. Put into pre-allocated storage.
00395   /*
00396     \param LocalRow - (In) Global row number for which indices are desired.
00397     \param Indices - (Out) Global column indices corresponding to values.
00398     \param Values - (Out) Matrix values.
00399     \param NumEntries - (Out) Number of indices.
00400     
00401     Note: A std::runtime_error exception is thrown if either \c Indices or \c Values is not large enough to hold the data associated
00402     with row \c GlobalRow. If \c GlobalRow does not belong to this node, then \c Indices and \c Values are unchanged and \c NumIndices is 
00403     returned as Teuchos::OrdinalTraits<size_t>::invalid().
00404   */
00405 template<class MatrixType>
00406 void OverlappingRowMatrix<MatrixType>::getGlobalRowCopy(GlobalOrdinal GlobalRow,
00407         const Teuchos::ArrayView<GlobalOrdinal> &Indices,
00408         const Teuchos::ArrayView<Scalar> &Values,
00409         size_t &NumEntries) const
00410 {
00411   LocalOrdinal LocalRow=RowMap_->getLocalElement(GlobalRow);
00412   if(LocalRow == Teuchos::OrdinalTraits<LocalOrdinal>::invalid()) NumEntries = Teuchos::OrdinalTraits<size_t>::invalid();
00413   else {
00414     if ((size_t)LocalRow < NumMyRowsA_) 
00415       A_->getGlobalRowCopy(GlobalRow,Indices,Values,NumEntries);
00416     else 
00417       ExtMatrix_->getGlobalRowCopy(GlobalRow,Indices,Values,NumEntries);
00418   }
00419 
00420 }
00421   
00422 //==========================================================================
00423 // Extract a list of entries in a specified local row of the graph. Put into storage allocated by calling routine.
00424   /*
00425     \param LocalRow - (In) Local row number for which indices are desired.
00426     \param Indices - (Out) Local column indices corresponding to values.
00427     \param Values - (Out) Matrix values.
00428     \param NumIndices - (Out) Number of indices.
00429     
00430     Note: A std::runtime_error exception is thrown if either \c Indices or \c Values is not large enough to hold the data associated
00431     with row \c LocalRow. If \c LocalRow is not valid for this node, then \c Indices and \c Values are unchanged and \c NumIndices is 
00432     returned as Teuchos::OrdinalTraits<size_t>::invalid().
00433   */
00434 template<class MatrixType>
00435 void OverlappingRowMatrix<MatrixType>::getLocalRowCopy(LocalOrdinal LocalRow, 
00436                    const Teuchos::ArrayView<LocalOrdinal> &Indices, 
00437                    const Teuchos::ArrayView<Scalar> &Values,
00438                    size_t &NumEntries) const
00439 {
00440   if ((size_t)LocalRow < NumMyRowsA_) 
00441     A_->getLocalRowCopy(LocalRow,Indices,Values,NumEntries);
00442   else 
00443     ExtMatrix_->getLocalRowCopy(LocalRow-(LocalOrdinal)(NumMyRowsA_),Indices,Values,NumEntries);
00444 }
00445   
00446 //==========================================================================
00447 // Extract a const, non-persisting view of global indices in a specified row of the matrix.
00448   /*
00449     \param GlobalRow - (In) Global row number for which indices are desired.
00450     \param Indices   - (Out) Global column indices corresponding to values.
00451     \param Values    - (Out) Row values
00452     \pre <tt>isLocallyIndexed() == false</tt>
00453     \post <tt>indices.size() == getNumEntriesInGlobalRow(GlobalRow)</tt>
00454     
00455     Note: If \c GlobalRow does not belong to this node, then \c indices is set to null.
00456   */
00457 template<class MatrixType>
00458 void OverlappingRowMatrix<MatrixType>::getGlobalRowView(GlobalOrdinal GlobalRow, 
00459               Teuchos::ArrayView<const GlobalOrdinal> &indices, 
00460               Teuchos::ArrayView<const Scalar> &values) const
00461 {
00462   LocalOrdinal LocalRow=RowMap_->getLocalElement(GlobalRow);
00463   if(LocalRow == Teuchos::OrdinalTraits<LocalOrdinal>::invalid())  {
00464     indices=Teuchos::null; 
00465     values =Teuchos::null;
00466   }
00467   else {
00468     if ((size_t)LocalRow < NumMyRowsA_) 
00469       A_->getGlobalRowView(GlobalRow,indices,values);
00470     else 
00471       ExtMatrix_->getGlobalRowView(GlobalRow,indices,values);
00472   }
00473 }
00474   
00475 //==========================================================================
00476 // Extract a const, non-persisting view of local indices in a specified row of the matrix.
00477   /*
00478     \param LocalRow - (In) Local row number for which indices are desired.
00479     \param Indices  - (Out) Global column indices corresponding to values.
00480     \param Values   - (Out) Row values
00481     \pre <tt>isGloballyIndexed() == false</tt>
00482     \post <tt>indices.size() == getNumEntriesInLocalRow(LocalRow)</tt>
00483     
00484     Note: If \c LocalRow does not belong to this node, then \c indices is set to null.
00485   */
00486 template<class MatrixType>
00487 void OverlappingRowMatrix<MatrixType>::getLocalRowView(LocalOrdinal LocalRow, 
00488                    Teuchos::ArrayView<const LocalOrdinal> &indices, 
00489                    Teuchos::ArrayView<const Scalar> &values) const
00490 {
00491   if ((size_t)LocalRow < NumMyRowsA_) 
00492     A_->getLocalRowView(LocalRow,indices,values);
00493   else 
00494     ExtMatrix_->getLocalRowView(LocalRow-(LocalOrdinal)(NumMyRowsA_),indices,values);
00495 }
00496   
00497 //==========================================================================
00498 //  Get a copy of the diagonal entries owned by this node, with local row indices.
00499   /* Returns a distributed Vector object partitioned according to this matrix's row map, containing the 
00500     the zero and non-zero diagonals owned by this node. */
00501 template<class MatrixType>
00502 void OverlappingRowMatrix<MatrixType>::getLocalDiagCopy(Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &diag) const
00503 {
00504   throw std::runtime_error("Ifpack2::OverlappingRowMatrix::getLocalDiagCopy not supported.");
00505 }
00506   
00507 //==========================================================================
00508 template<class MatrixType>
00509 void OverlappingRowMatrix<MatrixType>::leftScale(const Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x)
00510 {
00511   throw std::runtime_error("Ifpack2::OverlappingRowMatrix does not support leftScale.");
00512 }
00513   
00514 //==========================================================================
00515 template<class MatrixType>
00516 void OverlappingRowMatrix<MatrixType>::rightScale(const Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x)
00517 {
00518   throw std::runtime_error("Ifpack2::OverlappingRowMatrix does not support leftScale.");
00519 }
00520   
00521 //==========================================================================
00522 // Returns the Frobenius norm of the matrix. 
00523 template<class MatrixType>
00524 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType OverlappingRowMatrix<MatrixType>::getFrobeniusNorm() const
00525 {
00526   throw std::runtime_error("Ifpack2::OverlappingRowMatrix does not support getFrobeniusNorm.");
00527 }
00528 
00529 
00530 
00531 //==========================================================================
00532   // \brief Computes the operator-multivector application.
00533   /* Loosely, performs \f$Y = \alpha \cdot A^{\textrm{mode}} \cdot X + \beta \cdot Y\f$. However, the details of operation
00534     vary according to the values of \c alpha and \c beta. Specifically
00535     - if <tt>beta == 0</tt>, apply() <b>must</b> overwrite \c Y, so that any values in \c Y (including NaNs) are ignored.
00536     - if <tt>alpha == 0</tt>, apply() <b>may</b> short-circuit the operator, so that any values in \c X (including NaNs) are ignored.
00537 
00538     This is analagous to the *Multiply* function in Ifpack, not the *Apply*
00539   */
00540 template<class MatrixType>
00541 void OverlappingRowMatrix<MatrixType>::apply(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X, 
00542                Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y, 
00543                Teuchos::ETransp mode, 
00544                Scalar alpha,
00545                Scalar beta) const
00546 {
00547   // Note: This isn't AztecOO compliant.  But neither was Ifpack's version.
00548   TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
00549            "Ifpack2::OverlappingRowMatrix::apply ERROR: X.getNumVectors() != Y.getNumVectors().");
00550 
00551   Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
00552   Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar> > x_ptr = X.get2dView();
00553   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> >       y_ptr = Y.get2dViewNonConst();
00554   Y.putScalar(zero);
00555   size_t NumVectors = Y.getNumVectors();
00556 
00557   for (size_t i = 0 ; i < NumMyRowsA_ ; ++i) {
00558     size_t Nnz;
00559     // Use this class's getrow to make the below code simpler
00560     A_->getLocalRowCopy(i,Indices_(),Values_(),Nnz);
00561     if (mode==Teuchos::NO_TRANS){
00562       for (size_t j = 0 ; j < Nnz ; ++j) 
00563   for (size_t k = 0 ; k < NumVectors ; ++k)
00564     y_ptr[k][i] += Values_[j] * x_ptr[k][Indices_[j]];      
00565     }
00566     else if (mode==Teuchos::TRANS){
00567       for (size_t j = 0 ; j < Nnz ; ++j) 
00568   for (size_t k = 0 ; k < NumVectors ; ++k)
00569     y_ptr[k][Indices_[j]] += Values_[j] * x_ptr[k][i];
00570     }
00571     else { //mode==Teuchos::CONJ_TRANS
00572       for (size_t j = 0 ; j < Nnz ; ++j) 
00573   for (size_t k = 0 ; k < NumVectors ; ++k)
00574     y_ptr[k][Indices_[j]] += Teuchos::ScalarTraits<Scalar>::conjugate(Values_[j]) * x_ptr[k][i];
00575     }
00576   }
00577 
00578  for (size_t i = 0 ; i < NumMyRowsB_ ; ++i) {
00579     size_t Nnz;
00580     // Use this class's getrow to make the below code simpler
00581     ExtMatrix_->getLocalRowCopy(i,Indices_(),Values_(),Nnz);
00582     if (mode==Teuchos::NO_TRANS){
00583       for (size_t j = 0 ; j < Nnz ; ++j) 
00584   for (size_t k = 0 ; k < NumVectors ; ++k)
00585     y_ptr[k][NumMyRowsA_+i] += Values_[j] * x_ptr[k][Indices_[j]];      
00586     }
00587     else if (mode==Teuchos::TRANS){
00588       for (size_t j = 0 ; j < Nnz ; ++j) 
00589   for (size_t k = 0 ; k < NumVectors ; ++k)
00590     y_ptr[k][NumMyRowsA_+Indices_[j]] += Values_[j] * x_ptr[k][i];
00591     }
00592     else { //mode==Teuchos::CONJ_TRANS
00593       for (size_t j = 0 ; j < Nnz ; ++j) 
00594   for (size_t k = 0 ; k < NumVectors ; ++k)
00595     y_ptr[k][NumMyRowsA_+Indices_[j]] += Teuchos::ScalarTraits<Scalar>::conjugate(Values_[j]) * x_ptr[k][i];
00596     }
00597   }
00598 }
00599   
00600 
00601 
00602 // ======================================================================
00603 template<class MatrixType>
00604 void OverlappingRowMatrix<MatrixType>::importMultiVector(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X, 
00605                Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &OvX, 
00606                Tpetra::CombineMode CM)
00607 {
00608   OvX.doImport(X,*Importer_,CM);
00609 }
00610 
00611 
00612 // ======================================================================
00613 template<class MatrixType>
00614 void OverlappingRowMatrix<MatrixType>::exportMultiVector(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &OvX, 
00615                Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X, 
00616                Tpetra::CombineMode CM)
00617 {
00618   X.doExport(OvX,*Importer_,CM);
00619 }
00620 
00621 //==========================================================================
00622 // Indicates whether this operator supports applying the adjoint operator.
00623 template<class MatrixType>
00624 bool OverlappingRowMatrix<MatrixType>::hasTransposeApply() const
00625 {
00626   return true;
00627 }
00628 
00629 //==========================================================================
00630 template<class MatrixType> 
00631 bool OverlappingRowMatrix<MatrixType>::supportsRowViews() const
00632 {
00633   return false;
00634 }
00635 
00636 
00637 } // namespace Ifpack2
00638 
00639 #endif // IFPACK2_OVERLAPPINGROWMATRIX_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends