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 // mfh 03 Feb 2013: In a definition outside the class like this, the
00075 // return value is considered outside the class scope (for things like
00076 // resolving typedefs), but the arguments are considered inside the
00077 // class scope.
00078 template<class Scalar, 
00079      class LocalOrdinal,
00080      class GlobalOrdinal, 
00081      class Node, 
00082      class SpMatOps>
00083 Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps> >
00084 RowMatrixTransposer<Scalar, LocalOrdinal, GlobalOrdinal, Node, SpMatOps>::
00085 createTranspose()
00086 {
00087   using Teuchos::Array;
00088   using Teuchos::ArrayView;
00089   using Teuchos::ParameterList;
00090   using Teuchos::parameterList;
00091   using Teuchos::RCP;
00092   using Teuchos::rcp;
00093   using Tpetra::Import;
00094   using Tpetra::Export;
00095   typedef LocalOrdinal LO;
00096 
00097   //
00098   // This transpose is based upon the approach in EpetraExt.
00099   //
00100   size_t numLocalCols = origMatrix_->getNodeNumCols();
00101   size_t numLocalRows = origMatrix_->getNodeNumRows();
00102   size_t numLocalNnz  = origMatrix_->getNodeNumEntries();
00103 
00104   // Determine how many nonzeros there are per row in the transpose.
00105   Array<size_t> CurrentStart(numLocalCols,0);
00106   ArrayView<const LO> localIndices;
00107   ArrayView<const Scalar> localValues;
00108   RCP<const crs_matrix_type> crsMatrix = Teuchos::rcp_dynamic_cast<const crs_matrix_type>(origMatrix_);
00109   if (crsMatrix == Teuchos::null) {
00110     for (size_t i=0; i<numLocalRows; ++i) {
00111       const size_t numEntriesInRow = origMatrix_->getNumEntriesInLocalRow(i);
00112       origMatrix_->getLocalRowView(i, localIndices, localValues);
00113       for (size_t j=0; j<numEntriesInRow; ++j) {
00114         ++CurrentStart[ localIndices[j] ];
00115       }
00116     }
00117   } else {
00118     ArrayRCP<const size_t> origRowPtr_rcp;
00119     ArrayRCP<const LO>     origColInd_rcp;
00120     ArrayRCP<const Scalar> origValues_rcp;
00121     crsMatrix->getAllValues(origRowPtr_rcp, origColInd_rcp, origValues_rcp);
00122     ArrayView<const LO> origColInd = origColInd_rcp();
00123     for (LO j=0; j<origColInd.size(); ++j) {
00124       ++CurrentStart[ origColInd[j] ];
00125     }
00126   }
00127 
00128   //create temporary row-major storage for the transposed matrix
00129 
00130   ArrayRCP<size_t> rowptr_rcp(numLocalCols+1);
00131   ArrayRCP<LO>     colind_rcp(numLocalNnz);
00132   ArrayRCP<Scalar> values_rcp(numLocalNnz);
00133 
00134   // Since ArrayRCP's are slow...
00135   ArrayView<size_t> TransRowptr = rowptr_rcp();
00136   ArrayView<LO>     TransColind = colind_rcp();
00137   ArrayView<Scalar> TransValues = values_rcp();
00138 
00139   // Scansum the TransRowptr; reset CurrentStart
00140   TransRowptr[0]=0;
00141   for (size_t i=1; i<numLocalCols+1; ++i) TransRowptr[i]  = CurrentStart[i-1] + TransRowptr[i-1];
00142   for (size_t i=0; i<numLocalCols;   ++i) CurrentStart[i] = TransRowptr[i];
00143 
00144   //populate the row-major storage so that the data for the transposed matrix is easy to access
00145   if (crsMatrix == Teuchos::null) {
00146     for (size_t i=0; i<numLocalRows; ++i) {
00147       const size_t numEntriesInRow = origMatrix_->getNumEntriesInLocalRow(i);
00148       origMatrix_->getLocalRowView(i, localIndices, localValues);
00149 
00150       for (size_t j=0; j<numEntriesInRow; ++j) {
00151         size_t idx = CurrentStart[localIndices[j]];
00152         TransColind[idx] = Teuchos::as<LO>(i);
00153         TransValues[idx] = localValues[j];
00154         ++CurrentStart[localIndices[j]];
00155       }    
00156     } //for (size_t i=0; i<numLocalRows; ++i)
00157   } else {
00158     ArrayRCP<const size_t> origRowPtr_rcp;
00159     ArrayRCP<const LO>     origColInd_rcp;
00160     ArrayRCP<const Scalar> origValues_rcp;
00161     crsMatrix->getAllValues(origRowPtr_rcp, origColInd_rcp, origValues_rcp);
00162     ArrayView<const size_t>   origRowPtr = origRowPtr_rcp();
00163     ArrayView<const LO> origColInd = origColInd_rcp();
00164     ArrayView<const Scalar>   origValues = origValues_rcp();
00165     size_t k=0;
00166     for (LO i=0; i<origRowPtr.size()-1; ++i) {
00167       const LO rowIndex = Teuchos::as<LO>(i);
00168       for (size_t j=origRowPtr[i]; j<origRowPtr[i+1]; ++j) {
00169         size_t idx = CurrentStart[origColInd[k]];
00170         TransColind[idx] = rowIndex;
00171         TransValues[idx] = origValues[k];
00172         ++CurrentStart[origColInd[k++]];
00173       }
00174     }
00175   }
00176 
00177   //Allocate and populate temporary matrix with rows not uniquely owned
00178   RCP<crs_matrix_type> transMatrixWithSharedRows(new crs_matrix_type(origMatrix_->getColMap(),origMatrix_->getRowMap(),0));  
00179   transMatrixWithSharedRows->setAllValues(rowptr_rcp,colind_rcp,values_rcp);
00180 
00181 
00182   // Prebuild the importers and exporters the no-communication way, flipping the importers
00183   // and exporters around.
00184   RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> > myImport;
00185   RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > myExport;
00186   if(!origMatrix_->getGraph()->getImporter().is_null()) 
00187     myExport = rcp(new Export<LocalOrdinal,GlobalOrdinal,Node>(*origMatrix_->getGraph()->getImporter()));
00188   if(!origMatrix_->getGraph()->getExporter().is_null()) 
00189     myImport = rcp(new Import<LocalOrdinal,GlobalOrdinal,Node>(*origMatrix_->getGraph()->getExporter()));
00190 
00191   // Call ESFC
00192   transMatrixWithSharedRows->expertStaticFillComplete(origMatrix_->getRangeMap(),origMatrix_->getDomainMap(),myImport,myExport);
00193 
00194   // If transMatrixWithSharedRows has an exporter, that's what we want.  If it doesn't, the rows aren't actually shared,
00195   // and we're done!
00196   RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > exporter = transMatrixWithSharedRows->getGraph()->getExporter();
00197   if(exporter == Teuchos::null) {
00198     return transMatrixWithSharedRows;
00199   }
00200 
00201   // Finish using fusedexport
00202   return exportAndFillCompleteCrsMatrix<crs_matrix_type>(transMatrixWithSharedRows,*exporter);
00203 
00204 }
00205 //
00206 // Explicit instantiation macro
00207 //
00208 // Must be expanded from within the Tpetra namespace!
00209 //
00210 
00211 #define TPETRA_ROWMATRIXTRANSPOSER_INSTANT(SCALAR,LO,GO,NODE) \
00212   \
00213   template class RowMatrixTransposer< SCALAR, LO , GO , NODE >;
00214 
00215 
00216 }
00217 
00218 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines