Xpetra_StridedEpetraMap.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 //                    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 #include "Xpetra_ConfigDefs.hpp"
00050 
00051 #ifdef HAVE_XPETRA_EPETRA
00052 
00053 #include "Xpetra_StridedEpetraMap.hpp"
00054 #include "Xpetra_EpetraUtils.hpp"
00055 
00056 #include "Xpetra_EpetraExceptions.hpp"
00057 
00058 // MPI helper
00059 #define sumAll(rcpComm, in, out)                                        \
00060   Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_SUM, in, Teuchos::outArg(out));
00061 #define minAll(rcpComm, in, out)                                        \
00062   Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MIN, in, Teuchos::outArg(out));
00063 #define maxAll(rcpComm, in, out)                                        \
00064   Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MAX, in, Teuchos::outArg(out));
00065 
00066 namespace Xpetra {
00067 
00068   // Implementation note for constructors: the Epetra_Comm is cloned in the constructor of Epetra_BlockMap. We don't need to keep a reference on it.
00069   // TODO: use toEpetra() function here.
00070   StridedEpetraMap::StridedEpetraMap(global_size_t numGlobalElements, int indexBase, std::vector<size_t>& stridingInfo, const Teuchos::RCP<const Teuchos::Comm<int> > &comm,
00071                         LocalOrdinal stridedBlockId, GlobalOrdinal offset, LocalGlobal lg, const Teuchos::RCP<Node> &node)
00072   : EpetraMap(numGlobalElements, indexBase, comm, lg, node), StridedMap<int, int>(numGlobalElements, indexBase, stridingInfo, comm, stridedBlockId, offset)
00073   {    
00074     // check input data and reorganize map
00075     global_size_t numGlobalNodes = Teuchos::OrdinalTraits<global_size_t>::invalid();
00076     if(numGlobalElements != Teuchos::OrdinalTraits<global_size_t>::invalid())
00077       numGlobalNodes = numGlobalElements / getFixedBlockSize(); // number of nodes (over all processors)
00078     
00079     // build an equally distributed node map
00080     RCP<Epetra_Map> nodeMap = Teuchos::null;
00081     IF_EPETRA_EXCEPTION_THEN_THROW_GLOBAL_INVALID_ARG((nodeMap = (rcp(new Epetra_Map(static_cast<int>(numGlobalNodes), indexBase, *toEpetra(comm))))));
00082     
00083      // translate local node ids to local dofs
00084     int nStridedOffset = 0;
00085     int nDofsPerNode = Teuchos::as<int>(getFixedBlockSize()); // dofs per node for local striding block
00086     if(stridedBlockId > -1) {
00087       // determine nStridedOffset
00088       for(int j=0; j<stridedBlockId; j++) {
00089         nStridedOffset += stridingInfo[j];
00090       }
00091       nDofsPerNode = stridingInfo[stridedBlockId];
00092       
00093       numGlobalElements = nodeMap->NumGlobalElements()*nDofsPerNode;
00094     }
00095     std::vector<int> dofgids;
00096     for(int i = 0; i<nodeMap->NumMyElements(); i++) {
00097       int gid = nodeMap->GID(i);
00098       for(int dof = 0; dof < nDofsPerNode; ++dof) {
00099         // dofs are calculated by
00100         // global offset + node_GID * full size of strided map + striding offset of current striding block + dof id of current striding block
00101         dofgids.push_back(offset_ + gid*getFixedBlockSize() + nStridedOffset + dof);
00102       }
00103     }
00104 
00105     if (numGlobalElements == Teuchos::OrdinalTraits<global_size_t>::invalid()) {
00106       IF_EPETRA_EXCEPTION_THEN_THROW_GLOBAL_INVALID_ARG((map_ = (rcp(new Epetra_BlockMap(-1, dofgids.size(), &dofgids[0], 1, indexBase, *toEpetra(comm))))));
00107     } else {
00108       IF_EPETRA_EXCEPTION_THEN_THROW_GLOBAL_INVALID_ARG((map_ = (rcp(new Epetra_BlockMap(numGlobalElements, dofgids.size(), &dofgids[0], 1, indexBase, *toEpetra(comm))))));
00109     }
00110     
00111     TEUCHOS_TEST_FOR_EXCEPTION(map_->NumMyPoints() % nDofsPerNode != 0, Exceptions::RuntimeError, "StridedEpetraMap::StridedEpetraMap: wrong distribution of dofs among processors.");
00112     if(stridedBlockId == -1) {
00113       TEUCHOS_TEST_FOR_EXCEPTION(getNodeNumElements() != Teuchos::as<size_t>(nodeMap->NumMyElements()*nDofsPerNode), Exceptions::RuntimeError, "StridedEpetraMap::StridedEpetraMap: wrong distribution of dofs among processors.");
00114       TEUCHOS_TEST_FOR_EXCEPTION(getGlobalNumElements() != Teuchos::as<size_t>(nodeMap->NumGlobalElements()*nDofsPerNode), Exceptions::RuntimeError, "StridedEpetraMap::StridedEpetraMap: wrong distribution of dofs among processors.");
00115     } else {
00116       TEUCHOS_TEST_FOR_EXCEPTION(stridingInfo.size() < Teuchos::as<size_t>(stridedBlockId), Exceptions::RuntimeError, "StridedEpetraMap::StridedEpetraMap: stridedBlockId > stridingInfo.size()");      
00117       int nDofsInStridedBlock = stridingInfo[stridedBlockId];
00118       TEUCHOS_TEST_FOR_EXCEPTION(getNodeNumElements() != Teuchos::as<size_t>(nodeMap->NumMyElements()*nDofsInStridedBlock), Exceptions::RuntimeError, "StridedEpetraMap::StridedEpetraMap: wrong distribution of dofs among processors.");
00119       TEUCHOS_TEST_FOR_EXCEPTION(getGlobalNumElements() != Teuchos::as<size_t>(nodeMap->NumGlobalElements()*nDofsInStridedBlock), Exceptions::RuntimeError, "StridedEpetraMap::StridedEpetraMap: wrong distribution of dofs among processors.");
00120     }  
00121     
00122     TEUCHOS_TEST_FOR_EXCEPTION(CheckConsistency() == false, Exceptions::RuntimeError, "StridedEpetraMap::StridedEpetraMap: CheckConsistency() == false");      
00123   }
00124 
00125   StridedEpetraMap::StridedEpetraMap(global_size_t numGlobalElements, size_t numLocalElements, int indexBase,
00126                        std::vector<size_t>& stridingInfo, const Teuchos::RCP<const Teuchos::Comm<int> > &comm, LocalOrdinal stridedBlockId, GlobalOrdinal offset, const Teuchos::RCP<Node> &node)
00127   : EpetraMap(numGlobalElements, numLocalElements, indexBase, comm, node), StridedMap<int, int>(numGlobalElements, numLocalElements, indexBase, stridingInfo, comm, stridedBlockId, offset)
00128   {
00129     // check input data and reorganize map
00130 
00131     global_size_t numGlobalNodes = Teuchos::OrdinalTraits<global_size_t>::invalid();
00132     if(numGlobalElements != Teuchos::OrdinalTraits<global_size_t>::invalid())
00133       numGlobalNodes = numGlobalElements / getFixedBlockSize(); // number of nodes (over all processors)
00134     size_t blockSize = getFixedBlockSize();
00135     //if(stridedBlockId > -1) {
00136     //  blockSize = stridingInfo[stridedBlockId];
00137     //}
00138     size_t        numLocalNodes  = numLocalElements / blockSize;      // number of nodes (on each processor)
00139     
00140     // build an equally distributed node map
00141     RCP<Epetra_Map> nodeMap = Teuchos::null;
00142     IF_EPETRA_EXCEPTION_THEN_THROW_GLOBAL_INVALID_ARG((nodeMap = (rcp(new Epetra_Map(static_cast<int>(numGlobalNodes), numLocalNodes, indexBase, *toEpetra(comm))))));
00143     
00144     // translate local node ids to local dofs
00145     int nStridedOffset = 0;
00146     int nDofsPerNode = Teuchos::as<int>(getFixedBlockSize()); // dofs per node for local striding block
00147     if(stridedBlockId > -1) {
00148       // determine nStridedOffset
00149       for(int j=0; j<stridedBlockId; j++) {
00150         nStridedOffset += stridingInfo[j];
00151       }
00152       nDofsPerNode = stridingInfo[stridedBlockId];
00153       
00154       numGlobalElements = nodeMap->NumGlobalElements()*nDofsPerNode;
00155     }
00156     std::vector<int> dofgids;
00157     for(int i = 0; i<nodeMap->NumMyElements(); i++) {
00158       int gid = nodeMap->GID(i);
00159       for(int dof = 0; dof < nDofsPerNode; ++dof) {
00160   // dofs are calculated by
00161   // global offset + node_GID * full size of strided map + striding offset of current striding block + dof id of current striding block
00162         dofgids.push_back(offset_ + gid*getFixedBlockSize() + nStridedOffset + dof);
00163       }
00164     }
00165     
00166     if (numGlobalElements == Teuchos::OrdinalTraits<global_size_t>::invalid()) {
00167       IF_EPETRA_EXCEPTION_THEN_THROW_GLOBAL_INVALID_ARG((map_ = (rcp(new Epetra_BlockMap(-1, dofgids.size(), &dofgids[0], 1, indexBase, *toEpetra(comm))))));
00168     } else {
00169       IF_EPETRA_EXCEPTION_THEN_THROW_GLOBAL_INVALID_ARG((map_ = (rcp(new Epetra_BlockMap(numGlobalElements, dofgids.size(), &dofgids[0], 1, indexBase, *toEpetra(comm))))));
00170     }
00171     
00172     TEUCHOS_TEST_FOR_EXCEPTION(map_->NumMyPoints() % nDofsPerNode != 0, Exceptions::RuntimeError, "StridedEpetraMap::StridedEpetraMap: wrong distribution of dofs among processors.");
00173     if(stridedBlockId == -1) {
00174       TEUCHOS_TEST_FOR_EXCEPTION(getNodeNumElements() != Teuchos::as<size_t>(nodeMap->NumMyElements()*nDofsPerNode), Exceptions::RuntimeError, "StridedEpetraMap::StridedEpetraMap: wrong distribution of dofs among processors.");
00175       TEUCHOS_TEST_FOR_EXCEPTION(getGlobalNumElements() != Teuchos::as<size_t>(nodeMap->NumGlobalElements()*nDofsPerNode), Exceptions::RuntimeError, "StridedEpetraMap::StridedEpetraMap: wrong distribution of dofs among processors.");
00176     } else {
00177       TEUCHOS_TEST_FOR_EXCEPTION(stridingInfo.size() < Teuchos::as<size_t>(stridedBlockId), Exceptions::RuntimeError, "StridedEpetraMap::StridedEpetraMap: stridedBlockId > stridingInfo.size()");      
00178       int nDofsInStridedBlock = stridingInfo[stridedBlockId];
00179       TEUCHOS_TEST_FOR_EXCEPTION(getNodeNumElements() != Teuchos::as<size_t>(nodeMap->NumMyElements()*nDofsInStridedBlock), Exceptions::RuntimeError, "StridedEpetraMap::StridedEpetraMap: wrong distribution of dofs among processors.");
00180       TEUCHOS_TEST_FOR_EXCEPTION(getGlobalNumElements() != Teuchos::as<size_t>(nodeMap->NumGlobalElements()*nDofsInStridedBlock), Exceptions::RuntimeError, "StridedEpetraMap::StridedEpetraMap: wrong distribution of dofs among processors.");
00181     }    
00182     
00183     TEUCHOS_TEST_FOR_EXCEPTION(CheckConsistency() == false, Exceptions::RuntimeError, "StridedEpetraMap::StridedEpetraMap: CheckConsistency() == false");      
00184   }
00185     
00186 
00187   StridedEpetraMap::StridedEpetraMap(global_size_t numGlobalElements, const Teuchos::ArrayView<const int> &elementList, int indexBase,
00188                        std::vector<size_t>& stridingInfo, const Teuchos::RCP<const Teuchos::Comm<int> > &comm, LocalOrdinal stridedBlockId, const Teuchos::RCP<Node> &node)
00189   : EpetraMap(numGlobalElements, elementList, indexBase, comm, node), Xpetra::StridedMap<int,int>(numGlobalElements, elementList, indexBase, stridingInfo, comm, stridedBlockId)
00190   {
00191     int nDofsPerNode = Teuchos::as<int>(getFixedBlockSize());
00192     if(stridedBlockId != -1) {
00193       TEUCHOS_TEST_FOR_EXCEPTION(stridingInfo.size() < Teuchos::as<size_t>(stridedBlockId), Exceptions::RuntimeError, "StridedEpetraMap::StridedEpetraMap: stridedBlockId > stridingInfo.size()");
00194       nDofsPerNode = Teuchos::as<int>(stridingInfo[stridedBlockId]);
00195     }
00196     TEUCHOS_TEST_FOR_EXCEPTION(map_->NumMyPoints() % nDofsPerNode != 0, Exceptions::RuntimeError, "StridedEpetraMap::StridedEpetraMap: wrong distribution of dofs among processors.");
00197 
00198     // create EpetraMap using the dofs from ElementList
00199     IF_EPETRA_EXCEPTION_THEN_THROW_GLOBAL_INVALID_ARG((map_ = (rcp(new Epetra_BlockMap(numGlobalElements, elementList.size(), &elementList[0], 1, indexBase, *toEpetra(comm))))));
00200 
00201     // set parameters for striding information
00202     TEUCHOS_TEST_FOR_EXCEPTION(CheckConsistency() == false, Exceptions::RuntimeError, "StridedEpetraMap::StridedEpetraMap: CheckConsistency() == false");
00203   }
00204 
00205 
00206   std::string StridedEpetraMap::description() const {
00207     std::ostringstream oss;
00208     oss << EpetraMap::description();
00209     oss << "{getGlobalNumElements() = " << EpetraMap::getGlobalNumElements()
00210         << ", getNodeNumElements() = " << EpetraMap::getNodeNumElements()
00211         << ", isContiguous() = " << EpetraMap::isContiguous()
00212         << ", isDistributed() = " << EpetraMap::isDistributed()
00213         << "}";
00214     return oss.str();
00215   }
00216 
00217   void StridedEpetraMap::describe( Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const {
00218       
00219     const Teuchos::RCP<const Teuchos::Comm<int> > comm_ = EpetraMap::getComm();
00220 
00221     // This implementation come from Tpetra_Map_def.hpp (without modification)
00222     using std::endl;
00223     using std::setw;
00224     using Teuchos::VERB_DEFAULT;
00225     using Teuchos::VERB_NONE;
00226     using Teuchos::VERB_LOW;
00227     using Teuchos::VERB_MEDIUM;
00228     using Teuchos::VERB_HIGH;
00229     using Teuchos::VERB_EXTREME;
00230      
00231     const size_t nME = EpetraMap::getNodeNumElements();
00232     Teuchos::ArrayView<const int> myEntries = EpetraMap::getNodeElementList();
00233     int myImageID = comm_->getRank();
00234     int numImages = comm_->getSize();
00235      
00236     Teuchos::EVerbosityLevel vl = verbLevel;
00237     if (vl == VERB_DEFAULT) vl = VERB_LOW;
00238      
00239     size_t width = 1;
00240     for (size_t dec=10; dec<EpetraMap::getGlobalNumElements(); dec *= 10) {
00241       ++width;
00242     }
00243     width = std::max<size_t>(width,12) + 2;
00244      
00245     Teuchos::OSTab tab(out);
00246      
00247     if (vl == VERB_NONE) {
00248       // do nothing
00249     }
00250     else if (vl == VERB_LOW) {
00251       out << this->description() << endl;
00252     }
00253     else {  // MEDIUM, HIGH or EXTREME
00254       for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
00255         if (myImageID == imageCtr) {
00256           if (myImageID == 0) { // this is the root node (only output this info once)
00257             out << endl
00258                 << "Number of Global Entries = " << EpetraMap::getGlobalNumElements()  << endl
00259                 << "Maximum of all GIDs      = " << EpetraMap::getMaxAllGlobalIndex() << endl
00260                 << "Minimum of all GIDs      = " << EpetraMap::getMinAllGlobalIndex() << endl
00261                 << "Index Base               = " << EpetraMap::getIndexBase()         << endl;
00262           }
00263           out << endl;
00264           if (vl == VERB_HIGH || vl == VERB_EXTREME) {
00265             out << "Number of Local Elements   = " << nME           << endl
00266                 << "Maximum of my GIDs         = " << EpetraMap::getMaxGlobalIndex() << endl
00267                 << "Minimum of my GIDs         = " << EpetraMap::getMinGlobalIndex() << endl;
00268             out << endl;
00269           }
00270           if (vl == VERB_EXTREME) {
00271             out << std::setw(width) << "Node ID"
00272                 << std::setw(width) << "Local Index"
00273                 << std::setw(width) << "Global Index"
00274                 << endl;
00275             for (size_t i=0; i < nME; i++) {
00276               out << std::setw(width) << myImageID
00277                   << std::setw(width) << i
00278                   << std::setw(width) << myEntries[i]
00279                   << endl;
00280             }
00281             out << std::flush;
00282           }
00283         }
00284         // Do a few global ops to give I/O a chance to complete
00285         comm_->barrier();
00286         comm_->barrier();
00287         comm_->barrier();
00288       }
00289     }
00290   }
00291 
00292 }
00293 
00294 #endif