Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_CreateOverlapGraph.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 // 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 
00043 #ifndef IFPACK2_CREATEOVERLAPGRAPH_HPP
00044 #define IFPACK2_CREATEOVERLAPGRAPH_HPP
00045 
00046 #include "Ifpack2_ConfigDefs.hpp"
00047 #include "Tpetra_CrsGraph.hpp"
00048 #include "Tpetra_CrsMatrix.hpp"
00049 #include "Tpetra_Import.hpp"
00050 #include "Teuchos_RefCountPtr.hpp"
00051 
00052 
00053 namespace Ifpack2 {
00054 
00071 template<class GraphType>
00072 Teuchos::RCP<const GraphType>
00073 createOverlapGraph (const Teuchos::RCP<const GraphType>& inputGraph,
00074                     const int overlapLevel)
00075 {
00076   using Teuchos::RCP;
00077   using Teuchos::rcp;
00078   typedef Tpetra::Map<typename GraphType::local_ordinal_type,
00079                       typename GraphType::global_ordinal_type,
00080                       typename GraphType::node_type> map_type;
00081   typedef Tpetra::Import<typename GraphType::local_ordinal_type,
00082                          typename GraphType::global_ordinal_type,
00083                          typename GraphType::node_type> import_type;
00084   TEUCHOS_TEST_FOR_EXCEPTION(
00085     overlapLevel < 0, std::invalid_argument,
00086     "Ifpack2::createOverlapGraph: overlapLevel must be >= 0, "
00087     "but you specified overlapLevel = " << overlapLevel << ".");
00088 
00089   const int numProcs = inputGraph->getMap ()->getComm ()->getSize ();
00090   if (overlapLevel == 0 || numProcs < 2) {
00091     return inputGraph;
00092   }
00093 
00094   RCP<const map_type> overlapRowMap = inputGraph->getRowMap ();
00095   RCP<const map_type> domainMap = inputGraph->getDomainMap ();
00096   RCP<const map_type> rangeMap = inputGraph->getRangeMap ();
00097 
00098   RCP<GraphType> overlapGraph;
00099   RCP<const GraphType> oldGraph;
00100   RCP<const map_type> oldRowMap;
00101   for (int level = 0; level < overlapLevel; ++level) {
00102     oldGraph = overlapGraph;
00103     oldRowMap = overlapRowMap;
00104 
00105     RCP<const import_type> overlapImporter;
00106     if (level == 0) {
00107       overlapImporter = inputGraph->getImporter ();
00108     } else {
00109       overlapImporter = oldGraph->getImporter ();
00110     }
00111 
00112     overlapRowMap = overlapImporter->getTargetMap ();
00113     if (level < overlapLevel - 1) {
00114       overlapGraph = rcp (new GraphType (overlapRowMap, 0));
00115     }
00116     else {
00117       // On last iteration, we want to filter out all columns except those that
00118       // correspond to rows in the graph.  This ensures that our graph is square
00119       overlapGraph = rcp (new GraphType (overlapRowMap, overlapRowMap, 0));
00120     }
00121 
00122     overlapGraph->doImport (*inputGraph, *overlapImporter, Tpetra::INSERT);
00123     overlapGraph->fillComplete (domainMap, rangeMap);
00124   }
00125 
00126   return overlapGraph;
00127 }
00128 
00134 template<class GraphType>
00135 TEUCHOS_DEPRECATED Teuchos::RCP<const GraphType>
00136 CreateOverlapGraph (const Teuchos::RCP<const GraphType>& inputGraph,
00137                     int OverlapLevel)
00138 {
00139   return createOverlapGraph (inputGraph, OverlapLevel);
00140 }
00141 
00151 template<class MatrixType>
00152 Teuchos::RCP<const MatrixType>
00153 createOverlapMatrix (const Teuchos::RCP<const MatrixType>& inputMatrix,
00154                      const int overlapLevel)
00155 {
00156   using Teuchos::RCP;
00157   using Teuchos::rcp;
00158   typedef typename MatrixType::map_type map_type;
00159   typedef Tpetra::Import<typename MatrixType::local_ordinal_type,
00160     typename MatrixType::global_ordinal_type,
00161     typename MatrixType::node_type> import_type;
00162 
00163   TEUCHOS_TEST_FOR_EXCEPTION(
00164     overlapLevel < 0, std::invalid_argument,
00165     "Ifpack2::createOverlapMatrix: overlapLevel must be >= 0, "
00166     "but you specified overlapLevel = " << overlapLevel << ".");
00167 
00168   const int numProcs = inputMatrix->getMap ()->getComm ()->getSize ();
00169   if (overlapLevel == 0 || numProcs < 2) {
00170     return inputMatrix;
00171   }
00172 
00173   RCP<const map_type> overlapRowMap = inputMatrix->getRowMap ();
00174   RCP<const map_type> domainMap = inputMatrix->getDomainMap ();
00175   RCP<const map_type> rangeMap = inputMatrix->getRangeMap ();
00176 
00177   RCP<MatrixType> overlapMatrix;
00178   RCP<const MatrixType> oldMatrix;
00179   RCP<const map_type> oldRowMap;
00180   for (int level = 0; level < overlapLevel; ++level) {
00181     oldMatrix = overlapMatrix;
00182     oldRowMap = overlapRowMap;
00183 
00184     RCP<const import_type> overlapImporter;
00185     if (level == 0) {
00186       overlapImporter = inputMatrix->getGraph ()->getImporter ();
00187     } else {
00188       overlapImporter = oldMatrix->getGraph ()->getImporter ();
00189     }
00190 
00191     overlapRowMap = overlapImporter->getTargetMap ();
00192     if (level < overlapLevel - 1) {
00193       overlapMatrix = rcp (new MatrixType (overlapRowMap, 0));
00194     }
00195     else {
00196       // On last iteration, we want to filter out all columns except those that
00197       // correspond to rows in the matrix.  This assures that our matrix is square
00198       overlapMatrix = rcp (new MatrixType (overlapRowMap, overlapRowMap, 0));
00199     }
00200 
00201     overlapMatrix->doImport (*inputMatrix, *overlapImporter, Tpetra::INSERT);
00202     overlapMatrix->fillComplete (domainMap, rangeMap);
00203   }
00204 
00205   return overlapMatrix;
00206 }
00207 
00213 template<class MatrixType>
00214 TEUCHOS_DEPRECATED Teuchos::RCP<const MatrixType>
00215 CreateOverlapMatrix (const Teuchos::RCP<const MatrixType>& inputGraph,
00216                      int overlapLevel)
00217 {
00218   return createOverlapMatrix (inputGraph, overlapLevel);
00219 }
00220 
00221 } // namespace Ifpack2
00222 
00223 #endif // IFPACK2_CREATEOVERLAPGRAPH_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends