Xpetra_EpetraCrsMatrix.cpp Source File

Go to the documentation of this file.
00001 // @HEADER
00002 //
00003 // ***********************************************************************
00004 //
00005 //             Xpetra: A linear algebra interface package
00006 //                  Copyright 2012 Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact
00039 //                    Jonathan Hu       (jhu@sandia.gov)
00040 //                    Andrey Prokopenko (aprokop@sandia.gov)
00041 //                    Ray Tuminaro      (rstumin@sandia.gov)
00042 //
00043 // ***********************************************************************
00044 //
00045 // @HEADER
00046 #include <Teuchos_Array.hpp>
00047 #include "Xpetra_EpetraCrsMatrix.hpp"
00048 
00049 namespace Xpetra {
00050 
00051   template<class EpetraGlobalOrdinal>
00052   EpetraCrsMatrixT<EpetraGlobalOrdinal>::EpetraCrsMatrixT(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &rowMap, size_t maxNumEntriesPerRow, ProfileType pftype, const Teuchos::RCP< Teuchos::ParameterList > &plist)
00053     : mtx_(Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra(rowMap), maxNumEntriesPerRow, toEpetra(pftype)))), isFillResumed_(false) { }
00054 
00055   template<class EpetraGlobalOrdinal>
00056   EpetraCrsMatrixT<EpetraGlobalOrdinal>::EpetraCrsMatrixT(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &rowMap, const ArrayRCP< const size_t > &NumEntriesPerRowToAlloc, ProfileType pftype, const Teuchos::RCP< Teuchos::ParameterList > &plist)
00057     : isFillResumed_(false)
00058   {
00059     Teuchos::Array<int> numEntriesPerRowToAlloc(NumEntriesPerRowToAlloc.begin(), NumEntriesPerRowToAlloc.end()); // convert array of "size_t" to array of "int"
00060     mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra(rowMap), numEntriesPerRowToAlloc.getRawPtr(), toEpetra(pftype)));
00061   }
00062 
00063   template<class EpetraGlobalOrdinal>
00064   EpetraCrsMatrixT<EpetraGlobalOrdinal>::EpetraCrsMatrixT(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &rowMap, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &colMap, size_t maxNumEntriesPerRow, ProfileType pftype, const Teuchos::RCP< Teuchos::ParameterList > &plist)
00065     : mtx_(Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra(rowMap), toEpetra(colMap), maxNumEntriesPerRow, toEpetra(pftype)))), isFillResumed_(false) { }
00066 
00067   template<class EpetraGlobalOrdinal>
00068   EpetraCrsMatrixT<EpetraGlobalOrdinal>::EpetraCrsMatrixT(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &rowMap, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &colMap, const ArrayRCP< const size_t > &NumEntriesPerRowToAlloc, ProfileType pftype, const Teuchos::RCP< Teuchos::ParameterList > &plist)
00069     : isFillResumed_(false)
00070   {
00071     Teuchos::Array<int> numEntriesPerRowToAlloc(NumEntriesPerRowToAlloc.begin(), NumEntriesPerRowToAlloc.end()); // convert array of "size_t" to array of "int"
00072     mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra(rowMap), toEpetra(colMap), numEntriesPerRowToAlloc.getRawPtr(), toEpetra(pftype)));
00073   }
00074 
00075   template<class EpetraGlobalOrdinal>
00076   EpetraCrsMatrixT<EpetraGlobalOrdinal>::EpetraCrsMatrixT(const Teuchos::RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph, const Teuchos::RCP< Teuchos::ParameterList > &plist)
00077     : mtx_(Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra(graph)))), isFillResumed_(false) { }
00078 
00079   template<class EpetraGlobalOrdinal>
00080   EpetraCrsMatrixT<EpetraGlobalOrdinal>::EpetraCrsMatrixT(const EpetraCrsMatrixT& matrix)
00081     : mtx_(Teuchos::rcp(new Epetra_CrsMatrix(*(matrix.mtx_)))), isFillResumed_(false) { }
00082 
00083 
00084   template<class EpetraGlobalOrdinal>
00085   EpetraCrsMatrixT<EpetraGlobalOrdinal>::EpetraCrsMatrixT(const Teuchos::RCP<const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& sourceMatrix,
00086                                    const Import<LocalOrdinal,GlobalOrdinal,Node> &importer,
00087                                    const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& domainMap,
00088                                    const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rangeMap,
00089                                    const Teuchos::RCP<Teuchos::ParameterList>& params):
00090     isFillResumed_(false)
00091   {
00092     XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal>, *sourceMatrix, tSourceMatrix, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraCrsMatrixT as an input argument.");
00093     XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal>, importer, tImporter, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraImportT as an input argument.");
00094 
00095     const Epetra_Map* myDomainMap = (domainMap!=Teuchos::null)? &toEpetra(domainMap): 0;
00096     const Epetra_Map* myRangeMap  = (rangeMap !=Teuchos::null)? &toEpetra(rangeMap) : 0;
00097 
00098     // Follows the Tpetra parameters
00099     bool restrictComm=false;
00100     if(!params.is_null()) restrictComm = params->get("Restrict Communicator",restrictComm);
00101     mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(*tSourceMatrix.getEpetra_CrsMatrix(),*tImporter.getEpetra_Import(),myDomainMap,myRangeMap,restrictComm));
00102     if(restrictComm && mtx_->NumMyRows()==0)
00103       mtx_=Teuchos::null;
00104   }
00105 
00106   template<class EpetraGlobalOrdinal>
00107   EpetraCrsMatrixT<EpetraGlobalOrdinal>::EpetraCrsMatrixT(const Teuchos::RCP<const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& sourceMatrix,
00108                                    const Export<LocalOrdinal,GlobalOrdinal,Node> &exporter,
00109                                    const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& domainMap,
00110                                    const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rangeMap,
00111                                    const Teuchos::RCP<Teuchos::ParameterList>& params):
00112     isFillResumed_(false)
00113   {
00114     XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal>, *sourceMatrix, tSourceMatrix, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraCrsMatrixT as an input argument.");
00115     XPETRA_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal>, exporter, tExporter, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraExportT as an input argument.");
00116 
00117     const Epetra_Map* myDomainMap = (domainMap!=Teuchos::null)? &toEpetra(domainMap): 0;
00118     const Epetra_Map* myRangeMap  = (rangeMap !=Teuchos::null)? &toEpetra(rangeMap) : 0;
00119 
00120     // Follows the Tpetra parameters
00121     bool restrictComm=false;
00122     if(!params.is_null()) restrictComm = params->get("Restrict Communicator",restrictComm);
00123 
00124     mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(*tSourceMatrix.getEpetra_CrsMatrix(),*tExporter.getEpetra_Export(),myDomainMap,myRangeMap,restrictComm));
00125   }
00126 
00127 
00128 
00129   template<class EpetraGlobalOrdinal>
00130   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::insertGlobalValues(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal> &cols, const ArrayView<const Scalar> &vals) {
00131     XPETRA_MONITOR("EpetraCrsMatrixT::insertGlobalValues");
00132     XPETRA_ERR_CHECK(mtx_->InsertGlobalValues(globalRow, vals.size(), vals.getRawPtr(), cols.getRawPtr()));
00133   }
00134 
00135   template<class EpetraGlobalOrdinal>
00136   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::insertLocalValues(LocalOrdinal localRow, const ArrayView<const LocalOrdinal> &cols, const ArrayView<const Scalar> &vals) {
00137     XPETRA_MONITOR("EpetraCrsMatrixT::insertLocalValues");
00138     XPETRA_ERR_CHECK(mtx_->InsertMyValues(localRow, vals.size(), vals.getRawPtr(), cols.getRawPtr()));
00139   }
00140 
00141   template<class EpetraGlobalOrdinal>
00142   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::replaceGlobalValues(GlobalOrdinal globalRow, const ArrayView< const GlobalOrdinal > &indices, const ArrayView< const Scalar > &values) {
00143     XPETRA_MONITOR("EpetraCrsMatrixT::replaceGlobalValues");
00144 
00145     {
00146       const std::string tfecfFuncName("replaceGlobalValues");
00147       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(! isFillActive(), std::runtime_error,
00148                                             ": Fill must be active in order to call this method.  If you have already "
00149                                             "called fillComplete(), you need to call resumeFill() before you can "
00150                                             "replace values.");
00151 
00152       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(values.size() != indices.size(),
00153                                             std::runtime_error, ": values.size() must equal indices.size().");
00154     }
00155 
00156     XPETRA_ERR_CHECK(mtx_->ReplaceGlobalValues(globalRow, indices.size(), values.getRawPtr(), indices.getRawPtr()));
00157 
00158   }
00159 
00160   template<class EpetraGlobalOrdinal>
00161   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::replaceLocalValues(LocalOrdinal localRow, const ArrayView< const LocalOrdinal > &indices, const ArrayView< const Scalar > &values) {
00162     XPETRA_MONITOR("EpetraCrsMatrixT::replaceLocalValues");
00163 
00164     {
00165       const std::string tfecfFuncName("replaceLocalValues");
00166       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(! isFillActive(), std::runtime_error,
00167                                             ": Fill must be active in order to call this method.  If you have already "
00168                                             "called fillComplete(), you need to call resumeFill() before you can "
00169                                             "replace values.");
00170 
00171       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(values.size() != indices.size(),
00172                                             std::runtime_error, ": values.size() must equal indices.size().");
00173     }
00174 
00175     XPETRA_ERR_CHECK(mtx_->ReplaceMyValues(localRow, indices.size(), values.getRawPtr(), indices.getRawPtr()));
00176 
00177   }
00178 
00179   template<class EpetraGlobalOrdinal>
00180   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::allocateAllValues(size_t numNonZeros, ArrayRCP<size_t>& rowptr, ArrayRCP<LocalOrdinal>& colind, ArrayRCP<Scalar>& values) {
00181      XPETRA_MONITOR("EpetraCrsMatrixT::allocateAllValues");
00182 
00183     // Row offsets
00184     // Unfortunately, we cannot do this in the same manner as column indices
00185     // and values (see below).  The problem is that Tpetra insists on using
00186     // size_t, and Epetra uses int internally.  So we only resize here, and
00187     // will need to copy in setAllValues
00188     rowptr.resize(getNodeNumRows()+1);
00189 
00190     int  lowerOffset = 0;
00191     bool ownMemory   = false;
00192 
00193     // Column indices
00194     // Extract, resize, set colind
00195     Epetra_IntSerialDenseVector& myColind = mtx_->ExpertExtractIndices();
00196     myColind.Resize(numNonZeros);
00197     colind = Teuchos::arcp(myColind.Values(), lowerOffset, numNonZeros, ownMemory);
00198 
00199     // Values
00200     // Extract, reallocate, set values
00201     double *& myValues = mtx_->ExpertExtractValues();
00202     delete [] myValues;
00203     myValues = new double[numNonZeros];
00204     values = Teuchos::arcp(myValues,lowerOffset,numNonZeros,ownMemory);
00205   }
00206 
00207   template<class EpetraGlobalOrdinal>
00208   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::setAllValues(const ArrayRCP<size_t>& rowptr, const ArrayRCP<LocalOrdinal>& colind, const ArrayRCP<Scalar>& values) {
00209     XPETRA_MONITOR("EpetraCrsMatrixT::setAllValues");
00210 
00211     // Check sizes
00212     TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(rowptr.size()) != getNodeNumRows()+1, Xpetra::Exceptions::RuntimeError,
00213                                "An exception is thrown to let you know that the size of your rowptr array is incorrect.");
00214     TEUCHOS_TEST_FOR_EXCEPTION(values.size() != colind.size(), Xpetra::Exceptions::RuntimeError,
00215                                "An exception is thrown to let you know that you mismatched your pointers.");
00216 
00217     // Check pointers
00218     if (values.size() > 0) {
00219       TEUCHOS_TEST_FOR_EXCEPTION(colind.getRawPtr() != mtx_->ExpertExtractIndices().Values(), Xpetra::Exceptions::RuntimeError,
00220                                  "An exception is thrown to let you know that you mismatched your pointers.");
00221       TEUCHOS_TEST_FOR_EXCEPTION(values.getRawPtr() != mtx_->ExpertExtractValues(), Xpetra::Exceptions::RuntimeError,
00222                                  "An exception is thrown to let you know that you mismatched your pointers.");
00223     }
00224 
00225     // We have to make a copy here, it is unavoidable
00226     // See comments in allocateAllValues
00227     const size_t N = getNodeNumRows();
00228 
00229     Epetra_IntSerialDenseVector& myRowptr = mtx_->ExpertExtractIndexOffset();
00230     myRowptr.Resize(N+1);
00231     for (size_t i = 0; i < N+1; i++)
00232       myRowptr[i] = Teuchos::as<int>(rowptr[i]);
00233   }
00234 
00235   template<class EpetraGlobalOrdinal>
00236   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::getAllValues(ArrayRCP<const size_t>& rowptr, ArrayRCP<const LocalOrdinal>& colind, ArrayRCP<const Scalar>& values) const {
00237     XPETRA_MONITOR("EpetraCrsMatrixT::getAllValues");
00238 
00239     int  lowerOffset = 0;
00240     bool ownMemory   = false;
00241 
00242     const size_t n   = getNodeNumRows();
00243     const size_t nnz = getNodeNumEntries();
00244 
00245     // Row offsets
00246     // We have to make a copy here, it is unavoidable (see comments in allocateAllValues)
00247     Epetra_IntSerialDenseVector& myRowptr = mtx_->ExpertExtractIndexOffset();
00248     rowptr.resize(n+1);
00249     for (size_t i = 0; i < n+1; i++)
00250       (*const_cast<size_t*>(&rowptr[i])) = Teuchos::as<size_t>(myRowptr[i]);
00251 
00252     // Column indices
00253     colind = Teuchos::arcp(mtx_->ExpertExtractIndices().Values(), lowerOffset, nnz, ownMemory);
00254 
00255     // Values
00256     values = Teuchos::arcp(mtx_->ExpertExtractValues(), lowerOffset, nnz, ownMemory);
00257   }
00258 
00259   template<class EpetraGlobalOrdinal>
00260   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::resumeFill(const RCP< ParameterList > &params) {
00261     XPETRA_MONITOR("EpetraCrsMatrixT::resumeFill");
00262 
00263     // According to Tpetra documentation, resumeFill() may be called repeatedly.
00264     isFillResumed_ = true;
00265   }
00266 
00267   template<class EpetraGlobalOrdinal>
00268   bool EpetraCrsMatrixT<EpetraGlobalOrdinal>::isFillComplete() const { XPETRA_MONITOR("EpetraCrsMatrixT::isFillComplete"); if (isFillResumed_) return false; else return mtx_->Filled(); }
00269 
00270   template<class EpetraGlobalOrdinal>
00271   bool EpetraCrsMatrixT<EpetraGlobalOrdinal>::isFillActive() const { XPETRA_MONITOR("EpetraCrsMatrixT::isFillActive"); return !isFillComplete(); }
00272 
00273   template<class EpetraGlobalOrdinal>
00274   bool EpetraCrsMatrixT<EpetraGlobalOrdinal>::supportsRowViews() const { XPETRA_MONITOR("EpetraCrsMatrixT::supportsRowViews"); return true; }
00275 
00276   //TODO: throw same exception as Tpetra
00277   template<class EpetraGlobalOrdinal>
00278   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::getLocalRowCopy(LocalOrdinal LocalRow, const ArrayView<LocalOrdinal> &Indices, const ArrayView<Scalar> &Values, size_t &NumEntries) const {
00279     XPETRA_MONITOR("EpetraCrsMatrixT::getLocalRowCopy");
00280 
00281     int numEntries = -1;
00282     XPETRA_ERR_CHECK(mtx_->ExtractMyRowCopy(LocalRow, Indices.size(), numEntries, Values.getRawPtr(), Indices.getRawPtr()));
00283     NumEntries = numEntries;
00284   }
00285 
00286   template<class EpetraGlobalOrdinal>
00287   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::getGlobalRowCopy(GlobalOrdinal GlobalRow, const ArrayView<GlobalOrdinal> &Indices, const ArrayView<Scalar> &Values, size_t &NumEntries) const {
00288     XPETRA_MONITOR("EpetraCrsMatrixT::getGlobalRowCopy");
00289 
00290     int numEntries = -1;
00291     XPETRA_ERR_CHECK(mtx_->ExtractGlobalRowCopy(GlobalRow, Indices.size(), numEntries, Values.getRawPtr(), Indices.getRawPtr()));
00292     NumEntries = numEntries;
00293   }
00294 
00295   template<class EpetraGlobalOrdinal>
00296   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal> &indices, ArrayView<const Scalar> &values) const {
00297     XPETRA_MONITOR("EpetraCrsMatrixT::getGlobalRowView");
00298 
00299     int      numEntries;
00300     double * eValues;
00301     GlobalOrdinal * eIndices;
00302 
00303     XPETRA_ERR_CHECK(mtx_->ExtractGlobalRowView(GlobalRow, numEntries, eValues, eIndices));
00304     if (numEntries == 0) { eValues = NULL; eIndices = NULL; } // Cf. TEUCHOS_TEST_FOR_EXCEPT( p == 0 && size_in != 0 ) in Teuchos ArrayView constructor.
00305 
00306     indices = ArrayView<const GlobalOrdinal>(eIndices, numEntries);
00307     values  = ArrayView<const double>(eValues, numEntries);
00308   }
00309 
00310   template<class EpetraGlobalOrdinal>
00311   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal> &indices, ArrayView<const Scalar> &values) const {
00312     XPETRA_MONITOR("EpetraCrsMatrixT::getLocalRowView");
00313 
00314     int      numEntries;
00315     double * eValues;
00316     int    * eIndices;
00317 
00318     XPETRA_ERR_CHECK(mtx_->ExtractMyRowView(LocalRow, numEntries, eValues, eIndices));
00319     if (numEntries == 0) { eValues = NULL; eIndices = NULL; } // Cf. TEUCHOS_TEST_FOR_EXCEPT( p == 0 && size_in != 0 ) in Teuchos ArrayView constructor.
00320 
00321     indices = ArrayView<const int>(eIndices, numEntries);
00322     values  = ArrayView<const double>(eValues, numEntries);
00323   }
00324 
00325   template<class EpetraGlobalOrdinal>
00326   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::apply(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Y, Teuchos::ETransp mode, Scalar alpha, Scalar beta) const {
00327     XPETRA_MONITOR("EpetraCrsMatrixT::apply");
00328 
00329     //TEUCHOS_TEST_FOR_EXCEPTION((alpha != 1) || (beta != 0), Xpetra::Exceptions::NotImplemented, "Xpetra::EpetraCrsMatrixT.multiply() only accept alpha==1 and beta==0");
00330 
00331     XPETRA_DYNAMIC_CAST(const EpetraMultiVectorT<GlobalOrdinal>, X, eX, "Xpetra::EpetraCrsMatrixT->apply() only accept Xpetra::EpetraMultiVectorT as input arguments.");
00332     XPETRA_DYNAMIC_CAST(      EpetraMultiVectorT<GlobalOrdinal>, Y, eY, "Xpetra::EpetraCrsMatrixT->apply() only accept Xpetra::EpetraMultiVectorT as input arguments.");
00333 
00334     TEUCHOS_TEST_FOR_EXCEPTION((mode != Teuchos::NO_TRANS) && (mode != Teuchos::TRANS), Xpetra::Exceptions::NotImplemented, "Xpetra::EpetraCrsMatrixT->apply() only accept mode == NO_TRANS or mode == TRANS");
00335     bool eTrans = toEpetra(mode);
00336 
00337     // /!\ UseTranspose value
00338     TEUCHOS_TEST_FOR_EXCEPTION(mtx_->UseTranspose(), Xpetra::Exceptions::NotImplemented, "An exception is throw to let you know that Xpetra::EpetraCrsMatrixT->apply() do not take into account the UseTranspose() parameter of Epetra_CrsMatrix.");
00339 
00340     RCP<Epetra_MultiVector> epY = eY.getEpetra_MultiVector();
00341 
00342     // helper vector: tmp = A*x
00343     RCP<Epetra_MultiVector> tmp = Teuchos::rcp(new Epetra_MultiVector(*epY));
00344     tmp->PutScalar(0.0);
00345     XPETRA_ERR_CHECK(mtx_->Multiply(eTrans, *eX.getEpetra_MultiVector(), *tmp));
00346 
00347     // calculate alpha * A * x + beta * y
00348     XPETRA_ERR_CHECK(eY.getEpetra_MultiVector()->Update(alpha,*tmp,beta));
00349   }
00350 
00351   template<class EpetraGlobalOrdinal>
00352   std::string EpetraCrsMatrixT<EpetraGlobalOrdinal>::description() const {
00353     XPETRA_MONITOR("EpetraCrsMatrixT::description");
00354 
00355     // This implementation come from Tpetra_CrsMatrix_def.hpp (without modification)
00356     std::ostringstream oss;
00357     //TODO: oss << DistObject<char, LocalOrdinal,GlobalOrdinal>::description();
00358     if (isFillComplete()) {
00359       oss << "{status = fill complete"
00360           << ", global rows = " << getGlobalNumRows()
00361           << ", global cols = " << getGlobalNumCols()
00362           << ", global num entries = " << getGlobalNumEntries()
00363           << "}";
00364     }
00365     else {
00366       oss << "{status = fill not complete"
00367           << ", global rows = " << getGlobalNumRows()
00368           << "}";
00369     }
00370     return oss.str();
00371 
00372   }
00373 
00374   template<class EpetraGlobalOrdinal>
00375   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const {
00376     XPETRA_MONITOR("EpetraCrsMatrixT::describe");
00377 
00378     // This implementation come from Tpetra_CrsMatrix_def.hpp (without modification)
00379     using std::endl;
00380     using std::setw;
00381     using Teuchos::VERB_DEFAULT;
00382     using Teuchos::VERB_NONE;
00383     using Teuchos::VERB_LOW;
00384     using Teuchos::VERB_MEDIUM;
00385     using Teuchos::VERB_HIGH;
00386     using Teuchos::VERB_EXTREME;
00387     Teuchos::EVerbosityLevel vl = verbLevel;
00388     if (vl == VERB_DEFAULT) vl = VERB_LOW;
00389     RCP<const Comm<int> > comm = this->getComm();
00390     const int myImageID = comm->getRank(),
00391       numImages = comm->getSize();
00392     size_t width = 1;
00393     for (size_t dec=10; dec<getGlobalNumRows(); dec *= 10) {
00394       ++width;
00395     }
00396     width = std::max<size_t>(width,11) + 2;
00397     Teuchos::OSTab tab(out);
00398     //    none: print nothing
00399     //     low: print O(1) info from node 0
00400     //  medium: print O(P) info, num entries per node
00401     //    high: print O(N) info, num entries per row
00402     // extreme: print O(NNZ) info: print indices and values
00403     //
00404     // for medium and higher, print constituent objects at specified verbLevel
00405     if (vl != VERB_NONE) {
00406       if (myImageID == 0) out << this->description() << std::endl;
00407       // O(1) globals, minus what was already printed by description()
00408       if (isFillComplete() && myImageID == 0) {
00409         out << "Global number of diagonals = " << getGlobalNumDiags() << std::endl;
00410         out << "Global max number of entries = " << getGlobalMaxNumRowEntries() << std::endl;
00411       }
00412       // constituent objects
00413       if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
00414         if (myImageID == 0) out << "\nRow map: " << std::endl;
00415         getRowMap()->describe(out,vl);
00416         //
00417         if (getColMap() != null) {
00418           if (getColMap() == getRowMap()) {
00419             if (myImageID == 0) out << "\nColumn map is row map.";
00420           }
00421           else {
00422             if (myImageID == 0) out << "\nColumn map: " << std::endl;
00423             getColMap()->describe(out,vl);
00424           }
00425         }
00426         if (getDomainMap() != null) {
00427           if (getDomainMap() == getRowMap()) {
00428             if (myImageID == 0) out << "\nDomain map is row map.";
00429           }
00430           else if (getDomainMap() == getColMap()) {
00431             if (myImageID == 0) out << "\nDomain map is row map.";
00432           }
00433           else {
00434             if (myImageID == 0) out << "\nDomain map: " << std::endl;
00435             getDomainMap()->describe(out,vl);
00436           }
00437         }
00438         if (getRangeMap() != null) {
00439           if (getRangeMap() == getDomainMap()) {
00440             if (myImageID == 0) out << "\nRange map is domain map." << std::endl;
00441           }
00442           else if (getRangeMap() == getRowMap()) {
00443             if (myImageID == 0) out << "\nRange map is row map." << std::endl;
00444           }
00445           else {
00446             if (myImageID == 0) out << "\nRange map: " << std::endl;
00447             getRangeMap()->describe(out,vl);
00448           }
00449         }
00450         if (myImageID == 0) out << std::endl;
00451       }
00452       // O(P) data
00453       if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
00454         for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
00455           if (myImageID == imageCtr) {
00456             out << "Node ID = " << imageCtr << std::endl;
00457             // TODO: need a graph
00458             //               if (staticGraph_->indicesAreAllocated() == false) {
00459             //                 out << "Node not allocated" << std::endl;
00460             //               }
00461             //               else {
00462             //                 out << "Node number of allocated entries = " << staticGraph_->getNodeAllocationSize() << std::endl;
00463             //               }
00464 
00465             // TMP:
00466             //            const Epetra_CrsGraph & staticGraph_ = mtx_->Graph();
00467             // End of TMP
00468 
00469             out << "Node number of entries = " << getNodeNumEntries() << std::endl;
00470             if (isFillComplete()) {
00471               out << "Node number of diagonals = " << getNodeNumDiags() << std::endl;
00472             }
00473             out << "Node max number of entries = " << getNodeMaxNumRowEntries() << std::endl;
00474           }
00475           comm->barrier();
00476           comm->barrier();
00477           comm->barrier();
00478         }
00479       }
00480       // O(N) and O(NNZ) data
00481       if (vl == VERB_HIGH || vl == VERB_EXTREME) {
00482         for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
00483           if (myImageID == imageCtr) {
00484             out << std::setw(width) << "Node ID"
00485                 << std::setw(width) << "Global Row"
00486                 << std::setw(width) << "Num Entries";
00487             if (vl == VERB_EXTREME) {
00488               out << std::setw(width) << "(Index,Value)";
00489             }
00490             out << std::endl;
00491             for (size_t r=0; r < getNodeNumRows(); ++r) {
00492               const size_t nE = getNumEntriesInLocalRow(r);
00493               GlobalOrdinal gid = getRowMap()->getGlobalElement(r);
00494               out << std::setw(width) << myImageID
00495                   << std::setw(width) << gid
00496                   << std::setw(width) << nE;
00497               if (vl == VERB_EXTREME) {
00498                 if (isGloballyIndexed()) {
00499                   ArrayView<const GlobalOrdinal> rowinds;
00500                   ArrayView<const Scalar> rowvals;
00501                   getGlobalRowView(gid,rowinds,rowvals);
00502                   for (size_t j=0; j < nE; ++j) {
00503                     out << " (" << rowinds[j]
00504                         << ", " << rowvals[j]
00505                         << ") ";
00506                   }
00507                 }
00508                 else if (isLocallyIndexed()) {
00509                   ArrayView<const LocalOrdinal> rowinds;
00510                   ArrayView<const Scalar> rowvals;
00511                   getLocalRowView(r,rowinds,rowvals);
00512                   for (size_t j=0; j < nE; ++j) {
00513                     out << " (" << getColMap()->getGlobalElement(rowinds[j])
00514                         << ", " << rowvals[j]
00515                         << ") ";
00516                   }
00517                 }
00518               }
00519               out << std::endl;
00520             }
00521           }
00522           comm->barrier();
00523           comm->barrier();
00524           comm->barrier();
00525         }
00526       }
00527     }
00528 
00529   }
00530 
00531   // TODO: use toEpetra()
00532   template<class EpetraGlobalOrdinal>
00533   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::doImport(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node> &source,
00534                                  const Import<LocalOrdinal, GlobalOrdinal, Node> &importer, CombineMode CM) {
00535     XPETRA_MONITOR("EpetraCrsMatrixT::doImport");
00536 
00537     XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal>, source, tSource, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraCrsMatrixT as input arguments.");
00538     XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal>, importer, tImporter, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraImportT as input arguments.");
00539 
00540     RCP<const Epetra_CrsMatrix> v = tSource.getEpetra_CrsMatrix();
00541     int err = mtx_->Import(*v, *tImporter.getEpetra_Import(), toEpetra(CM));
00542     TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Catch error code returned by Epetra.");
00543   }
00544 
00545   template<class EpetraGlobalOrdinal>
00546   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::doExport(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node> &dest,
00547                                  const Import<LocalOrdinal, GlobalOrdinal, Node>& importer, CombineMode CM) {
00548     XPETRA_MONITOR("EpetraCrsMatrixT::doExport");
00549 
00550     XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal>, dest, tDest, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraCrsMatrixT as input arguments.");
00551     XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal>, importer, tImporter, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraImportT as input arguments.");
00552 
00553     RCP<const Epetra_CrsMatrix> v = tDest.getEpetra_CrsMatrix();
00554     int err = mtx_->Export(*v, *tImporter.getEpetra_Import(), toEpetra(CM));
00555     TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Catch error code returned by Epetra.");
00556   }
00557 
00558   template<class EpetraGlobalOrdinal>
00559   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::doImport(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node> &source,
00560                                  const Export<LocalOrdinal, GlobalOrdinal, Node>& exporter, CombineMode CM) {
00561     XPETRA_MONITOR("EpetraCrsMatrixT::doImport");
00562 
00563     XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal>, source, tSource, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraCrsMatrixT as input arguments.");
00564     XPETRA_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal>, exporter, tExporter, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraImportT as input arguments.");
00565 
00566     RCP<const Epetra_CrsMatrix> v = tSource.getEpetra_CrsMatrix();
00567     int err = mtx_->Import(*v, *tExporter.getEpetra_Export(), toEpetra(CM));
00568     TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Catch error code returned by Epetra.");
00569 
00570   }
00571 
00572   template<class EpetraGlobalOrdinal>
00573   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::doExport(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node> &dest,
00574                                  const Export<LocalOrdinal, GlobalOrdinal, Node>& exporter, CombineMode CM) {
00575     XPETRA_MONITOR("EpetraCrsMatrixT::doExport");
00576 
00577     XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal>, dest, tDest, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraCrsMatrixT as input arguments.");
00578     XPETRA_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal>, exporter, tExporter, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraImportT as input arguments.");
00579 
00580     RCP<const Epetra_CrsMatrix> v = tDest.getEpetra_CrsMatrix();
00581     int err = mtx_->Export(*v, *tExporter.getEpetra_Export(), toEpetra(CM));
00582     TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Catch error code returned by Epetra.");
00583   }
00584 
00585     template<class EpetraGlobalOrdinal>
00586     void EpetraCrsMatrixT<EpetraGlobalOrdinal>::fillComplete(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &domainMap,
00587                                        const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &rangeMap,
00588                                        const RCP< ParameterList > &params) {
00589       XPETRA_MONITOR("EpetraCrsMatrixT::fillComplete");
00590 
00591       // For Epetra matrices, resumeFill() is a fictive operation. There is no need for a fillComplete after some resumeFill() operations.
00592       if (isFillResumed_ == true) { isFillResumed_ = false; return; }
00593 
00594       bool doOptimizeStorage = true;
00595       if (params != null && params->get("Optimize Storage",true) == false) doOptimizeStorage = false;
00596       mtx_->FillComplete(toEpetra(domainMap), toEpetra(rangeMap), doOptimizeStorage);
00597     }
00598 
00599     template<class EpetraGlobalOrdinal>
00600     void EpetraCrsMatrixT<EpetraGlobalOrdinal>::fillComplete(const RCP< ParameterList > &params) {
00601       XPETRA_MONITOR("EpetraCrsMatrixT::fillComplete");
00602 
00603       // For Epetra matrices, resumeFill() is a fictive operation. There is no need for a fillComplete after some resumeFill() operations.
00604       if (isFillResumed_ == true) { isFillResumed_ = false; return; }
00605 
00606       bool doOptimizeStorage = true;
00607       if (params != null && params->get("Optimize Storage",true) == false) doOptimizeStorage = false;
00608       mtx_->FillComplete(doOptimizeStorage);
00609     }
00610 
00611 
00612   template<class EpetraGlobalOrdinal>
00613   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::replaceDomainMapAndImporter(const Teuchos::RCP< const  Map< LocalOrdinal, GlobalOrdinal, Node > >& newDomainMap, Teuchos::RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> >  & newImporter) {
00614       XPETRA_MONITOR("EpetraCrsMatrixT::replaceDomainMapAndImporter");
00615       XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal>, *newImporter, eImporter, "Xpetra::EpetraCrsMatrixT::replaceDomainMapAndImporter only accepts Xpetra::EpetraImportT.");
00616 
00617       const RCP<const Epetra_Import> & myImport = eImporter.getEpetra_Import();
00618       int rv=0;
00619       if(myImport==Teuchos::null)
00620         rv=mtx_->ReplaceDomainMapAndImporter( toEpetra(newDomainMap),0);
00621       else
00622         rv=mtx_->ReplaceDomainMapAndImporter( toEpetra(newDomainMap),&*myImport);
00623       TEUCHOS_TEST_FOR_EXCEPTION(rv != 0, std::runtime_error, "Xpetra::EpetraCrsMatrixT::replaceDomainMapAndImporter FAILED!");
00624   }
00625 
00626 
00627   template<class EpetraGlobalOrdinal>
00628   void EpetraCrsMatrixT<EpetraGlobalOrdinal>::expertStaticFillComplete(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & domainMap,
00629                                                  const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & rangeMap,
00630                                                  const RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> > &importer,
00631                                                  const RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > &exporter,
00632                                                  const RCP<ParameterList> & params) {
00633     XPETRA_MONITOR("EpetraCrsMatrixT::expertStaticFillComplete");
00634     int rv=0;
00635     const Epetra_Import * myimport =0;
00636     const Epetra_Export * myexport =0;
00637 
00638     if(!importer.is_null()) {
00639       XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal>, *importer, eImporter, "Xpetra::EpetraCrsMatrixT::expertStaticFillComplete only accepts Xpetra::EpetraImportT.");
00640       myimport = eImporter.getEpetra_Import().getRawPtr();
00641     }
00642     if(!exporter.is_null()) {
00643       XPETRA_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal>, *exporter, eExporter, "Xpetra::EpetraCrsMatrixT::expertStaticFillComplete only accepts Xpetra::EpetraImportT.");
00644       myexport = eExporter.getEpetra_Export().getRawPtr();
00645     }
00646 
00647     rv=mtx_->ExpertStaticFillComplete(toEpetra(domainMap), toEpetra(rangeMap), myimport, myexport);
00648 
00649     TEUCHOS_TEST_FOR_EXCEPTION(rv != 0, std::runtime_error, "Xpetra::EpetraCrsMatrixT::expertStaticFillComplete FAILED!");
00650   }
00651 
00652 #ifndef XPETRA_EPETRA_NO_32BIT_GLOBAL_INDICES
00653 template class EpetraCrsMatrixT<int>;
00654 #endif
00655 
00656 #ifndef XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES
00657 template class EpetraCrsMatrixT<long long>;
00658 #endif
00659 
00660 }