Tpetra Matrix/Vector Services Version of the Day
Tpetra_RowMatrixTransposer_def.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2008) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the 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 Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 
00042 #ifndef TPETRA_ROWMATRIXTRANSPOSER_DEF_HPP
00043 #define TPETRA_ROWMATRIXTRANSPOSER_DEF_HPP
00044 
00045 #include "Tpetra_Export.hpp"
00046 #include "Tpetra_Import.hpp"
00047 #include "Tpetra_Map.hpp"
00048 #include "Teuchos_DefaultSerialComm.hpp"
00049 #ifdef DOXYGEN_USE_ONLY
00050   // #include "Tpetra_RowMatrixtransposer_decl.hpp"
00051 #endif
00052 
00053 namespace Tpetra {
00054 
00055 template<class Scalar,
00056      class LocalOrdinal,
00057      class GlobalOrdinal,
00058      class Node,
00059      class SpMatOps>
00060 RowMatrixTransposer<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps>::
00061 RowMatrixTransposer (const Teuchos::RCP<const crs_matrix_type>& origMatrix)
00062   : origMatrix_ (origMatrix) {}
00063 
00064 template<class Scalar,
00065      class LocalOrdinal,
00066      class GlobalOrdinal,
00067      class Node,
00068      class SpMatOps>
00069 TEUCHOS_DEPRECATED
00070 RowMatrixTransposer<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps>::
00071 RowMatrixTransposer (const crs_matrix_type& origMatrix)
00072   : origMatrix_ (Teuchos::rcpFromRef (origMatrix)) {}
00073 
00074 template<class Scalar,
00075      class LocalOrdinal,
00076      class GlobalOrdinal,
00077      class Node,
00078      class SpMatOps>
00079 Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps> >
00080 RowMatrixTransposer<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps>::
00081 createTranspose()
00082 {
00083   using Teuchos::RCP;
00084 
00085   // Do the local transpose
00086   RCP<crs_matrix_type> transMatrixWithSharedRows = createTransposeLocal ();
00087 
00088   // If transMatrixWithSharedRows has an exporter, that's what we
00089   // want.  If it doesn't, the rows aren't actually shared, and we're
00090   // done!
00091   RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > exporter =
00092     transMatrixWithSharedRows->getGraph ()->getExporter ();
00093   if (exporter.is_null ()) {
00094     return transMatrixWithSharedRows;
00095   }
00096   else {
00097     // Use the Export object to do a fused Export and fillComplete.
00098     return exportAndFillCompleteCrsMatrix<crs_matrix_type> (transMatrixWithSharedRows, *exporter);
00099   }
00100 }
00101 
00102 
00103 // mfh 03 Feb 2013: In a definition outside the class like this, the
00104 // return value is considered outside the class scope (for things like
00105 // resolving typedefs), but the arguments are considered inside the
00106 // class scope.
00107 template<class Scalar,
00108          class LocalOrdinal,
00109          class GlobalOrdinal,
00110          class Node,
00111          class SpMatOps>
00112 Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps> >
00113 RowMatrixTransposer<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps>::
00114 createTransposeLocal ()
00115 {
00116   using Teuchos::Array;
00117   using Teuchos::ArrayRCP;
00118   using Teuchos::ArrayView;
00119   using Teuchos::RCP;
00120   using Teuchos::rcp;
00121   using Teuchos::rcp_dynamic_cast;
00122   typedef LocalOrdinal LO;
00123   typedef GlobalOrdinal GO;
00124   typedef Tpetra::Import<LO, GO, Node> import_type;
00125   typedef Tpetra::Export<LO, GO, Node> export_type;
00126 
00127   //
00128   // This transpose is based upon the approach in EpetraExt.
00129   //
00130   size_t numLocalCols = origMatrix_->getNodeNumCols();
00131   size_t numLocalRows = origMatrix_->getNodeNumRows();
00132   size_t numLocalNnz  = origMatrix_->getNodeNumEntries();
00133 
00134   // Determine how many nonzeros there are per row in the transpose.
00135   Array<size_t> CurrentStart(numLocalCols,0);
00136   ArrayView<const LO> localIndices;
00137   ArrayView<const Scalar> localValues;
00138   RCP<const crs_matrix_type> crsMatrix =
00139     rcp_dynamic_cast<const crs_matrix_type> (origMatrix_);
00140   if (crsMatrix == Teuchos::null) {
00141     for (size_t i=0; i<numLocalRows; ++i) {
00142       const size_t numEntriesInRow = origMatrix_->getNumEntriesInLocalRow(i);
00143       origMatrix_->getLocalRowView(i, localIndices, localValues);
00144       for (size_t j=0; j<numEntriesInRow; ++j) {
00145         ++CurrentStart[ localIndices[j] ];
00146       }
00147     }
00148   } else {
00149     ArrayRCP<const size_t> origRowPtr_rcp;
00150     ArrayRCP<const LO>     origColInd_rcp;
00151     ArrayRCP<const Scalar> origValues_rcp;
00152     crsMatrix->getAllValues(origRowPtr_rcp, origColInd_rcp, origValues_rcp);
00153     ArrayView<const LO> origColInd = origColInd_rcp();
00154     for (LO j=0; j<origColInd.size(); ++j) {
00155       ++CurrentStart[ origColInd[j] ];
00156     }
00157   }
00158 
00159   // create temporary row-major storage for the transposed matrix
00160 
00161   ArrayRCP<size_t> rowptr_rcp(numLocalCols+1);
00162   ArrayRCP<LO>     colind_rcp(numLocalNnz);
00163   ArrayRCP<Scalar> values_rcp(numLocalNnz);
00164 
00165   // Since ArrayRCP's are slow...
00166   ArrayView<size_t> TransRowptr = rowptr_rcp();
00167   ArrayView<LO>     TransColind = colind_rcp();
00168   ArrayView<Scalar> TransValues = values_rcp();
00169 
00170   // Scansum the TransRowptr; reset CurrentStart
00171   TransRowptr[0]=0;
00172   for (size_t i=1; i<numLocalCols+1; ++i) TransRowptr[i]  = CurrentStart[i-1] + TransRowptr[i-1];
00173   for (size_t i=0; i<numLocalCols;   ++i) CurrentStart[i] = TransRowptr[i];
00174 
00175   // populate the row-major storage so that the data for the transposed
00176   // matrix is easy to access
00177   if (crsMatrix == Teuchos::null) {
00178     for (size_t i=0; i<numLocalRows; ++i) {
00179       const size_t numEntriesInRow = origMatrix_->getNumEntriesInLocalRow (i);
00180       origMatrix_->getLocalRowView(i, localIndices, localValues);
00181 
00182       for (size_t j=0; j<numEntriesInRow; ++j) {
00183         size_t idx = CurrentStart[localIndices[j]];
00184         TransColind[idx] = Teuchos::as<LO>(i);
00185         TransValues[idx] = localValues[j];
00186         ++CurrentStart[localIndices[j]];
00187       }
00188     } //for (size_t i=0; i<numLocalRows; ++i)
00189   } else {
00190     ArrayRCP<const size_t> origRowPtr_rcp;
00191     ArrayRCP<const LO>     origColInd_rcp;
00192     ArrayRCP<const Scalar> origValues_rcp;
00193     crsMatrix->getAllValues(origRowPtr_rcp, origColInd_rcp, origValues_rcp);
00194     ArrayView<const size_t>   origRowPtr = origRowPtr_rcp();
00195     ArrayView<const LO> origColInd = origColInd_rcp();
00196     ArrayView<const Scalar>   origValues = origValues_rcp();
00197     size_t k=0;
00198     for (LO i=0; i<origRowPtr.size()-1; ++i) {
00199       const LO rowIndex = Teuchos::as<LO>(i);
00200       for (size_t j=origRowPtr[i]; j<origRowPtr[i+1]; ++j) {
00201         size_t idx = CurrentStart[origColInd[k]];
00202         TransColind[idx] = rowIndex;
00203         TransValues[idx] = origValues[k];
00204         ++CurrentStart[origColInd[k++]];
00205       }
00206     }
00207   }
00208 
00209   //Allocate and populate temporary matrix with rows not uniquely owned
00210   RCP<crs_matrix_type> transMatrixWithSharedRows =
00211     rcp (new crs_matrix_type (origMatrix_->getColMap (),
00212                               origMatrix_->getRowMap (), 0));
00213   transMatrixWithSharedRows->setAllValues (rowptr_rcp, colind_rcp, values_rcp);
00214 
00215   // Prebuild the importers and exporters the no-communication way,
00216   // flipping the importers and exporters around.
00217   RCP<const import_type> myImport;
00218   RCP<const export_type> myExport;
00219   if (! origMatrix_->getGraph ()->getImporter ().is_null ()) {
00220     myExport = rcp (new export_type (*origMatrix_->getGraph ()->getImporter ()));
00221   }
00222   if (! origMatrix_->getGraph ()->getExporter ().is_null ()) {
00223     myImport = rcp (new import_type (*origMatrix_->getGraph ()->getExporter ()));
00224   }
00225 
00226   // Call ESFC & return
00227   transMatrixWithSharedRows->expertStaticFillComplete (origMatrix_->getRangeMap (),
00228                                                        origMatrix_->getDomainMap (),
00229                                                        myImport, myExport);
00230   return transMatrixWithSharedRows;
00231 }
00232 //
00233 // Explicit instantiation macro
00234 //
00235 // Must be expanded from within the Tpetra namespace!
00236 //
00237 
00238 #define TPETRA_ROWMATRIXTRANSPOSER_INSTANT(SCALAR,LO,GO,NODE) \
00239   \
00240   template class RowMatrixTransposer< SCALAR, LO , GO , NODE >;
00241 
00242 
00243 }
00244 
00245 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines