Xpetra_StridedMapFactory.hpp 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 //                    Jeremie Gaidamour (jngaida@sandia.gov)
00040 //                    Jonathan Hu       (jhu@sandia.gov)
00041 //                    Ray Tuminaro      (rstumin@sandia.gov)
00042 //
00043 // ***********************************************************************
00044 //
00045 // @HEADER
00046 
00047 // WARNING: This code is experimental. Backwards compatibility should not be expected.
00048 
00049 #ifndef XPETRA_STRIDEDMAPFACTORY_HPP
00050 #define XPETRA_STRIDEDMAPFACTORY_HPP
00051 
00052 #include "Xpetra_ConfigDefs.hpp"
00053 
00054 #include "Xpetra_StridedMap.hpp"
00055 
00056 #ifdef HAVE_XPETRA_TPETRA
00057 #include "Xpetra_StridedTpetraMap.hpp"
00058 #endif
00059 #ifdef HAVE_XPETRA_EPETRA
00060 #include "Xpetra_StridedEpetraMap.hpp"
00061 #endif
00062 
00063 #include "Xpetra_Exceptions.hpp"
00064 
00065 // This factory creates Xpetra::Map. User have to specify the exact class of object that he want to create (ie: a Xpetra::TpetraMap or a Xpetra::EpetraMap).
00066 
00067 namespace Xpetra {
00068 
00069   template <class LocalOrdinal, class GlobalOrdinal = LocalOrdinal, class Node = Kokkos::DefaultNode::DefaultNodeType>
00070   class StridedMapFactory {
00071     
00072   private:
00074     StridedMapFactory() {}
00075     
00076   public:
00077     
00079     static Teuchos::RCP<StridedMap<LocalOrdinal,GlobalOrdinal, Node> > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, std::vector<size_t>& stridingInfo, const Teuchos::RCP<const Teuchos::Comm<int> > &comm, LocalOrdinal stridedBlockId=-1, GlobalOrdinal offset = 0, LocalGlobal lg=Xpetra::GloballyDistributed, const Teuchos::RCP<Node> &node = Kokkos::DefaultNode::getDefaultNode()) {
00080 
00081 #ifdef HAVE_XPETRA_TPETRA
00082       if (lib == UseTpetra)
00083         return Teuchos::rcp( new StridedTpetraMap<LocalOrdinal,GlobalOrdinal, Node> (numGlobalElements, indexBase, stridingInfo, comm, stridedBlockId, offset, lg, node) );
00084 #endif
00085 
00086       XPETRA_FACTORY_ERROR_IF_EPETRA(lib);
00087       XPETRA_FACTORY_END;
00088     }
00089 
00091     static Teuchos::RCP<StridedMap<LocalOrdinal,GlobalOrdinal, Node> > Build(UnderlyingLib lib, global_size_t numGlobalElements, size_t numLocalElements, GlobalOrdinal indexBase, std::vector<size_t>& stridingInfo, const Teuchos::RCP<const Teuchos::Comm<int> > &comm, LocalOrdinal stridedBlockId=-1, GlobalOrdinal offset = 0, const Teuchos::RCP<Node> &node = Kokkos::DefaultNode::getDefaultNode()) {
00092 
00093 #ifdef HAVE_XPETRA_TPETRA
00094       if (lib == UseTpetra)       
00095         return rcp( new StridedTpetraMap<LocalOrdinal,GlobalOrdinal, Node> (numGlobalElements, numLocalElements, indexBase, stridingInfo, comm, stridedBlockId, offset, node) );
00096 #endif
00097 
00098       XPETRA_FACTORY_ERROR_IF_EPETRA(lib);
00099       XPETRA_FACTORY_END;
00100     }
00101 
00102     static RCP<StridedMap<LocalOrdinal,GlobalOrdinal, Node> > Build(UnderlyingLib lib, const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >& map, std::vector<size_t>& stridingInfo, LocalOrdinal stridedBlockId=-1, GlobalOrdinal offset = 0) {
00103 
00104 #ifdef HAVE_XPETRA_TPETRA
00105       if (lib == UseTpetra) {
00106         Teuchos::RCP<const TpetraMap<LocalOrdinal,GlobalOrdinal,Node> > tmap = Teuchos::rcp_dynamic_cast<const TpetraMap<LocalOrdinal, GlobalOrdinal, Node> >(map);
00107         TEUCHOS_TEST_FOR_EXCEPTION(tmap == Teuchos::null, Exceptions::RuntimeError,"Xpetra::StridedMapFactory::Build: bad cast. map is not a TpetraMap.");
00108         return rcp( new StridedTpetraMap<LocalOrdinal,GlobalOrdinal, Node> (tmap->getTpetra_Map(), stridingInfo, stridedBlockId, offset) );
00109       }
00110 #endif
00111       XPETRA_FACTORY_ERROR_IF_EPETRA(lib);
00112       XPETRA_FACTORY_END;
00113     }
00114 
00115     // special constructor for generating a given subblock of a strided map
00116     static RCP<StridedMap<LocalOrdinal,GlobalOrdinal, Node> > Build(const RCP<const StridedMap<LocalOrdinal, GlobalOrdinal, Node> >& map, LocalOrdinal stridedBlockId) {
00117       TEUCHOS_TEST_FOR_EXCEPTION(stridedBlockId < 0, Exceptions::RuntimeError,"Xpetra::StridedMapFactory::Build: constructor expects stridedBlockId > -1.");
00118 #ifdef HAVE_XPETRA_TPETRA
00119       if (map->lib() == UseTpetra) {
00120         TEUCHOS_TEST_FOR_EXCEPTION(map->getStridedBlockId() != -1, Exceptions::RuntimeError,"Xpetra::StridedMapFactory::Build: constructor expects a full map (stridedBlockId == -1).");
00121         std::vector<size_t> stridingInfo = map->getStridingData();
00122 
00124         Teuchos::ArrayView< const GlobalOrdinal > dofGids = map->getNodeElementList();
00125         //std::sort(dofGids.begin(),dofGids.end()); // TODO: do i need this?
00126 
00127         // determine nStridedOffset
00128         size_t nStridedOffset = 0;
00129         for(int j=0; j<map->getStridedBlockId(); j++) {
00130           nStridedOffset += stridingInfo[j];
00131         }
00132 
00133         size_t numMyBlockDofs = stridingInfo[stridedBlockId] / map->getFixedBlockSize() * map->getNodeNumElements();
00134         std::vector<GlobalOrdinal> subBlockDofGids(numMyBlockDofs);
00135 
00136         // TODO fill vector with dofs
00137         typename Teuchos::ArrayView< const GlobalOrdinal >::iterator it;
00138         for(it = dofGids.begin(); it!=dofGids.end(); ++it) {
00139           if(map->GID2StridingBlockId( *it ) == stridedBlockId) {
00140             subBlockDofGids.push_back( *it );
00141           }
00142         }
00143 
00144         const Teuchos::ArrayView<const LocalOrdinal> subBlockDofGids_view(&subBlockDofGids[0],subBlockDofGids.size());
00145 
00146         // call constructor for TpetraMap
00147         return rcp( new StridedTpetraMap<LocalOrdinal,GlobalOrdinal, Node> (subBlockDofGids.size(), subBlockDofGids_view, map->getIndexBase(), stridingInfo, map->getComm(), stridedBlockId, map->getNode()) );
00149 
00150       }
00151 #endif
00152       XPETRA_FACTORY_ERROR_IF_EPETRA(map->lib());
00153       XPETRA_FACTORY_END;
00154     }
00155 
00156 #if 0  // TODO
00158     static Teuchos::RCP<StridedMap<LocalOrdinal,GlobalOrdinal, Node> > Build(UnderlyingLib lib, global_size_t numGlobalElements, const Teuchos::ArrayView<const GlobalOrdinal> &elementList, GlobalOrdinal indexBase, const Teuchos::RCP<const Teuchos::Comm<int> > &comm, const Teuchos::RCP<Node> &node = Kokkos::DefaultNode::getDefaultNode()) {
00159 
00160 #ifdef HAVE_XPETRA_TPETRA
00161       if (lib == UseTpetra) 
00162         return rcp( new StridedTpetraMap<LocalOrdinal,GlobalOrdinal, Node> (numGlobalElements, elementList, indexBase, comm, node) );
00163 #endif
00164 
00165       XPETRA_FACTORY_ERROR_IF_EPETRA(lib);
00166       XPETRA_FACTORY_END;
00167     }
00168 #endif
00169 
00170 
00171   };
00172 
00173   template <>
00174   class StridedMapFactory<int, int> {
00175 
00176     typedef int LocalOrdinal;
00177     typedef int GlobalOrdinal;
00178     typedef Kokkos::DefaultNode::DefaultNodeType Node;
00179     
00180   private:
00182     StridedMapFactory() {}
00183     
00184   public:
00185     
00186     static RCP<StridedMap<LocalOrdinal,GlobalOrdinal, Node> > Build(UnderlyingLib lib, global_size_t numGlobalElements, int indexBase, std::vector<size_t>& stridingInfo, const Teuchos::RCP<const Teuchos::Comm<int> > &comm, LocalOrdinal stridedBlockId=-1, GlobalOrdinal offset = 0, LocalGlobal lg=GloballyDistributed, const Teuchos::RCP<Kokkos::DefaultNode::DefaultNodeType> &node = Kokkos::DefaultNode::getDefaultNode()) {
00187 
00188 #ifdef HAVE_XPETRA_TPETRA
00189       if (lib == UseTpetra)
00190         return rcp( new StridedTpetraMap<LocalOrdinal,GlobalOrdinal, Node> (numGlobalElements, indexBase, stridingInfo, comm, stridedBlockId, offset, lg, node) );
00191 #endif
00192 
00193 #ifdef HAVE_XPETRA_EPETRA
00194       if (lib == UseEpetra)
00195         return rcp( new StridedEpetraMap(numGlobalElements, indexBase, stridingInfo, comm, stridedBlockId, offset, lg, node) );
00196 #endif
00197 
00198       XPETRA_FACTORY_END;
00199     }
00200 
00201     static RCP<StridedMap<LocalOrdinal,GlobalOrdinal, Node> > Build(UnderlyingLib lib, global_size_t numGlobalElements, size_t numLocalElements, int indexBase, std::vector<size_t>& stridingInfo, const Teuchos::RCP<const Teuchos::Comm<int> > &comm, LocalOrdinal stridedBlockId=-1, GlobalOrdinal offset = 0, const Teuchos::RCP<Kokkos::DefaultNode::DefaultNodeType> &node = Kokkos::DefaultNode::getDefaultNode()) {
00202 
00203 #ifdef HAVE_XPETRA_TPETRA
00204       if (lib == UseTpetra)       
00205         return rcp( new StridedTpetraMap<LocalOrdinal,GlobalOrdinal, Node> (numGlobalElements, numLocalElements, indexBase, stridingInfo, comm, stridedBlockId, offset, node) );
00206 #endif
00207 
00208 #ifdef HAVE_XPETRA_EPETRA
00209       if (lib == UseEpetra)
00210         return rcp( new StridedEpetraMap(numGlobalElements, numLocalElements, indexBase, stridingInfo, comm, stridedBlockId, offset, node) );
00211 #endif
00212 
00213       XPETRA_FACTORY_END;
00214     }
00215 
00216     static RCP<StridedMap<LocalOrdinal,GlobalOrdinal, Node> > Build(UnderlyingLib lib, const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >& map, std::vector<size_t>& stridingInfo, LocalOrdinal stridedBlockId=-1, GlobalOrdinal offset = 0) {
00217 
00218 #ifdef HAVE_XPETRA_TPETRA
00219       if (lib == UseTpetra) {
00220   Teuchos::RCP<const TpetraMap<LocalOrdinal,GlobalOrdinal,Node> > tmap = Teuchos::rcp_dynamic_cast<const TpetraMap<LocalOrdinal, GlobalOrdinal, Node> >(map);
00221   TEUCHOS_TEST_FOR_EXCEPTION(tmap == Teuchos::null, Exceptions::RuntimeError,"Xpetra::StridedMapFactory::Build: bad cast. map is not a TpetraMap.");
00222         return rcp( new StridedTpetraMap<LocalOrdinal,GlobalOrdinal, Node> (tmap->getTpetra_Map(), stridingInfo, stridedBlockId, offset) );
00223       }
00224 #endif
00225 
00226 #ifdef HAVE_XPETRA_EPETRA
00227       if (lib == UseEpetra) {
00228   Teuchos::RCP<const EpetraMap> emap = Teuchos::rcp_dynamic_cast<const EpetraMap>(map);
00229   TEUCHOS_TEST_FOR_EXCEPTION(emap == Teuchos::null, Exceptions::RuntimeError,"Xpetra::StridedMapFactory::Build: bad cast. map is not a EpetraMap.");  
00230         return rcp( new StridedEpetraMap(Teuchos::rcpFromRef(emap->getEpetra_Map()), stridingInfo,stridedBlockId, offset) ); // ugly: EpetraMap and TpetraMap differ here
00231       }
00232 #endif
00233 
00234       XPETRA_FACTORY_END;
00235     }
00236 
00237     // special constructor for generating a given subblock of a strided map
00238     static RCP<StridedMap<LocalOrdinal,GlobalOrdinal, Node> > Build(const RCP<const StridedMap<LocalOrdinal, GlobalOrdinal, Node> >& map, LocalOrdinal stridedBlockId) {
00239       TEUCHOS_TEST_FOR_EXCEPTION(stridedBlockId < 0, Exceptions::RuntimeError,"Xpetra::StridedMapFactory::Build: constructor expects stridedBlockId > -1.");
00240       typedef Xpetra::StridedMap<LocalOrdinal,GlobalOrdinal,Node> StridedMapClass;
00241       //typename Teuchos::ArrayView< const GlobalOrdinal >::iterator it;
00242 #ifdef HAVE_XPETRA_TPETRA
00243       if (map->lib() == UseTpetra) {
00244         TEUCHOS_TEST_FOR_EXCEPTION(map->getStridedBlockId() != -1, Exceptions::RuntimeError,"Xpetra::StridedMapFactory::Build: constructor expects a full map (stridedBlockId == -1).");
00245         std::vector<size_t> stridingInfo = map->getStridingData();
00246 
00248         Teuchos::ArrayView< const GlobalOrdinal > dofGids = map->getNodeElementList();
00249         //std::sort(dofGids.begin(),dofGids.end()); // TODO: do i need this?
00250 
00251         // determine nStridedOffset
00252         size_t nStridedOffset = 0;
00253         for(int j=0; j<map->getStridedBlockId(); j++) {
00254           nStridedOffset += stridingInfo[j];
00255         }
00256 
00257         size_t numMyBlockDofs = stridingInfo[stridedBlockId] / map->getFixedBlockSize() * map->getNodeNumElements();
00258         std::vector<GlobalOrdinal> subBlockDofGids(numMyBlockDofs);
00259 
00260         // TODO fill vector with dofs
00261 
00262         for(Teuchos::ArrayView< const GlobalOrdinal >::iterator it = dofGids.begin(); it!=dofGids.end(); ++it) {
00263           if(map->GID2StridingBlockId( *it ) == Teuchos::as<size_t>(stridedBlockId)) {
00264             subBlockDofGids.push_back( *it );
00265           }
00266         }
00267 
00268         const Teuchos::ArrayView<const LocalOrdinal> subBlockDofGids_view(&subBlockDofGids[0],subBlockDofGids.size());
00269 
00270         // call constructor for TpetraMap
00271         return rcp( new StridedTpetraMap<LocalOrdinal,GlobalOrdinal, Node> (/*subBlockDofGids.size()*/Teuchos::OrdinalTraits<global_size_t>::invalid(), subBlockDofGids_view, map->getIndexBase(), stridingInfo, map->getComm(), stridedBlockId, map->getNode()) );
00273 
00274       }
00275 #endif
00276 #ifdef HAVE_XPETRA_EPETRA
00277       if (map->lib() == UseEpetra) {
00278         TEUCHOS_TEST_FOR_EXCEPTION(map->getStridedBlockId() != -1, Exceptions::RuntimeError,"Xpetra::StridedMapFactory::Build: constructor expects a full map (stridedBlockId == -1).");
00279         std::vector<size_t> stridingInfo = map->getStridingData();
00280 
00282         Teuchos::ArrayView< const GlobalOrdinal > dofGids = map->getNodeElementList();
00283         //std::sort(dofGids.begin(),dofGids.end()); // TODO: do i need this?
00284 
00285         // determine nStridedOffset
00286         size_t nStridedOffset = 0;
00287         for(int j=0; j<map->getStridedBlockId(); j++) {
00288           nStridedOffset += stridingInfo[j];
00289         }
00290 
00291         size_t numMyBlockDofs = stridingInfo[stridedBlockId] / map->getFixedBlockSize() * map->getNodeNumElements();
00292         std::vector<GlobalOrdinal> subBlockDofGids(numMyBlockDofs);
00293 
00294         // TODO fill vector with dofs
00295         //Teuchos::ArrayView< const GlobalOrdinal >::iterator it;
00296         for(Teuchos::ArrayView< const GlobalOrdinal >::iterator it = dofGids.begin(); it!=dofGids.end(); ++it) {
00297           if(map->GID2StridingBlockId( *it ) == Teuchos::as<size_t>(stridedBlockId)) {
00298             subBlockDofGids.push_back( *it );
00299           }
00300         }
00301 
00302         const Teuchos::ArrayView<const LocalOrdinal> subBlockDofGids_view(&subBlockDofGids[0],subBlockDofGids.size());
00303 
00304         // call constructor for TpetraMap
00305         return rcp( new StridedEpetraMap(Teuchos::OrdinalTraits<global_size_t>::invalid(), subBlockDofGids_view, map->getIndexBase(), stridingInfo, map->getComm(), stridedBlockId, map->getNode()) );
00307 
00308       }
00309 #endif
00310       XPETRA_FACTORY_END;
00311     }
00312 
00313 #if 0 // TODO
00314     static RCP<StridedMap<LocalOrdinal,GlobalOrdinal, Node> > Build(UnderlyingLib lib, global_size_t numGlobalElements, const Teuchos::ArrayView<const int> &elementList, int indexBase, const Teuchos::RCP<const Teuchos::Comm<int> > &comm, const Teuchos::RCP<Kokkos::DefaultNode::DefaultNodeType> &node = Kokkos::DefaultNode::getDefaultNode()) {
00315 #ifdef HAVE_XPETRA_TPETRA
00316       if (lib == UseTpetra) 
00317         return rcp( new StridedTpetraMap<LocalOrdinal,GlobalOrdinal, Node> (numGlobalElements, elementList, indexBase, comm, node) );
00318 #endif
00319 
00320 #ifdef HAVE_XPETRA_EPETRA
00321       if (lib == UseEpetra)
00322         return rcp( new StridedEpetraMap(numGlobalElements, elementList, indexBase, comm, node) );
00323 #endif
00324       XPETRA_FACTORY_END;
00325     }
00326 #endif
00327 
00328  
00329   };
00330 
00331 }
00332 
00333 #define XPETRA_STRIDEDMAPFACTORY_SHORT
00334 #endif
00335 //TODO: removed unused methods