MueLu_Level.cpp

Go to the documentation of this file.

// @HEADER
//
// ***********************************************************************
//
//        MueLu: A package for multigrid based preconditioning
//                  Copyright 2012 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact
//                    Jonathan Hu       (jhu@sandia.gov)
//                    Andrey Prokopenko (aprokop@sandia.gov)
//                    Ray Tuminaro      (rstumin@sandia.gov)
//
// ***********************************************************************
//
// @HEADER
#include <Teuchos_TabularOutputter.hpp>

#include "MueLu_Level.hpp"

#include "MueLu_FactoryManagerBase.hpp"

namespace MueLu {

RCP<Level> Level::Build() {
  RCP<Level> newLevel = rcp(new Level());

  // Copy 'keep' status of variables
  for (TwoKeyMap::const_iterator kt = map_.begin(); kt != map_.end(); kt++) {
    const FactoryBase* factory = kt->first;

    for (SubMap::const_iterator it = kt->second.begin(); it != kt->second.end(); it++) {
      const std::string& ename = it->first;

      if (IsKept(ename, factory, MueLu::Keep)) {  // MueLu::Keep is the only flag propagated
        if (factory == NULL)                      // TODO: Is this possible?? Throw exception. Not supposed to use the FactoryManager here.
          newLevel->Keep(ename, NoFactory::get());
        else
          newLevel->Keep(ename, factory);
      }
    }
  }

  return newLevel;
}

int Level::GetLevelID() const { return levelID_; }

void Level::SetLevelID(int levelID) {
  if (levelID_ != -1 && levelID_ != levelID)
    GetOStream(Warnings1) << "Level::SetLevelID(): Changing an already defined LevelID (previousID=" << levelID_ << ", newID=" << levelID << ")" << std::endl;

  levelID_ = levelID;
}

void Level::SetPreviousLevel(const RCP<Level>& previousLevel) {
  if (previousLevel_ != Teuchos::null && previousLevel_ != previousLevel)
    GetOStream(Warnings1) << "Level::SetPreviousLevel(): PreviousLevel was already defined" << std::endl;

  previousLevel_ = previousLevel;
}

void Level::SetFactoryManager(const RCP<const FactoryManagerBase>& factoryManager) {
  factoryManager_ = factoryManager;
}

const RCP<const FactoryManagerBase> Level::GetFactoryManager() {
  return factoryManager_;
}

void Level::AddKeepFlag(const std::string& ename, const FactoryBase* factory, KeepType keep) {
  if (!IsKey(factory, ename)) {
    // If the entry does not exist, create it to store the keep flag
    Teuchos::RCP<MueLu::VariableContainer> newVar = Teuchos::rcp(new MueLu::VariableContainer);
    map_[factory][ename]                          = newVar;
  }
  // Set the flag
  map_[factory][ename]->AddKeepFlag(keep);
}

void Level::RemoveKeepFlag(const std::string& ename, const FactoryBase* factory, KeepType keep) {
  // No entry = nothing to do
  if (!IsKey(factory, ename))
    return;

  // Remove the flag
  Teuchos::RCP<MueLu::VariableContainer>& v = map_[factory][ename];
  v->RemoveKeepFlag(keep);

  // Remove data if no keep flag left and counter == 0
  if ((v->IsRequested() == false) && (v->GetKeepFlag() == 0)) {
    v = Teuchos::null;  // free data

    map_[factory].erase(ename);
    if (map_.count(factory) == 0)
      map_.erase(factory);
  }
}

KeepType Level::GetKeepFlag(const std::string& ename, const FactoryBase* factory) const {
  if (!IsKey(factory, ename))
    return false;

  return Get(factory, ename)->GetKeepFlag();
}

void Level::Request(const FactoryBase& factory) {
  RequestMode prev = requestMode_;
  requestMode_     = REQUEST;
  factory.CallDeclareInput(*this);
  requestMode_ = prev;
}

void Level::Release(const FactoryBase& factory) {
  RequestMode prev = requestMode_;
  requestMode_     = RELEASE;
  factory.CallDeclareInput(*this);
  requestMode_ = prev;
}

void Level::DeclareInput(const std::string& ename, const FactoryBase* factory, const FactoryBase* requestedBy) {
  if (requestMode_ == REQUEST) {
    try {
      Request(ename, factory, requestedBy);

    } catch (Exceptions::DependencyError& e) {
      std::ostringstream msg;
      msg << requestedBy->ShortClassName() << "::DeclareInput: (" << e.what() << ") unable to find or generate requested data \""
          << ename << "\" with generating factory \"" << ((factory != NULL) ? factory->ShortClassName() : "null") << "\" [" << factory << "]";
      msg << "\n    during request for data \"" << std::setw(15) << ename << "\" on level " << GetLevelID()
          << " by factory " << std::setw(25) << requestedBy->ShortClassName() << " [" << requestedBy << "]";
      throw Exceptions::RuntimeError(msg.str());

    } catch (Exceptions::RuntimeError& e) {
      std::ostringstream msg;
      msg << e.what() << "\n    during request for data \"" << std::setw(15) << ename << "\" on level " << GetLevelID()
          << " by factory " << std::setw(25) << requestedBy->ShortClassName() << " [" << requestedBy << "]";
      throw Exceptions::RuntimeError(msg.str());
    }

  } else if (requestMode_ == RELEASE) {
    Release(ename, factory, requestedBy);

  } else
    TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::Level::DeclareInput(): requestMode_ undefined.");
}

void Level::DeclareDependencies(const FactoryBase* factory, bool bRequestOnly, bool bReleaseOnly) {  // TODO: replace bReleaseOnly, bReleaseOnly by one RequestMode enum
  if (bRequestOnly && bReleaseOnly)
    TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::Level::DeclareDependencies(): Both bRequestOnly and bReleaseOnly set to true makes no sense.");

  if (requestMode_ == REQUEST) {
    if (bReleaseOnly == false) Request(*factory);

  } else if (requestMode_ == RELEASE) {
    if (bRequestOnly == false) Release(*factory);

  } else
    TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::Level::DeclareDependencies(): requestMode_ undefined.");
}

void Level::Request(const std::string& ename, const FactoryBase* factory, const FactoryBase* requestedBy) {
  const FactoryBase* fac = GetFactory(ename, factory);
  // printf("(l=%d) [%43s] requesting \"%20s\" generated by %10p  [actually, generated by %p (%43s)]\n",
  // levelID_, requestedBy->description().c_str(), ename.c_str(), factory, fac, fac->description().c_str());

  // We request factory only if necessary, i.e.
  //   - factory has not been requested before, and
  //   - we need data which is not available
  // Let us consider an RAP factory in the reuse scenario. The factory generates "A" (coarse matrix) and
  // "RAP graph" (pattern of A). Typically, "A" has keep flag Final, which is cleared during next Setup, but
  // "RAP graph" has flag Keep, which is not cleared as it is part of NextRun. Therefore, during next Setup
  // we have RAP factory with available "RAP graph" but not available "A".
  //
  // During the regular construction phase, we will do a single request: ("A", RAPFactory). Suppose, we used
  //    bool test = (IsRequestedFactory(fac) == false && IsAvailable(fac) == false);
  // This is incorrect, as IsAvailable(fac) checks whether there is *any* data generated by factory, which there is
  // ("A"), and the dependencies of the factory would not be requested despite the need for them (we need fine A, P, and
  // R to generate Ac even if we know sparsity pattern of Ac).
  //
  // On the other hand,
  //   bool test = (IsRequestedFactory(fac) == false && IsAvailable(ename, fac) == false);
  // is correct as ("A", fac) is not available (only ("RAP graph", fac) is), and dependent factories would be
  // properly requested.
  //
  // This way, factory is requested only once (because of the IsRequested(fac) check), and only when one of the needed
  // pieces of data is not availble.
  bool test = (IsRequestedFactory(fac) == false && IsAvailable(ename, fac) == false);

  // This request must be done before calling Request(*fac) to avoid circular dependency problems.
  if (!IsKey(fac, ename)) {
    Teuchos::RCP<MueLu::VariableContainer> newVar = Teuchos::rcp(new MueLu::VariableContainer);
    map_[fac][ename]                              = newVar;
  }

  Teuchos::RCP<MueLu::VariableContainer>& v = map_[fac][ename];
  v->Request(requestedBy);

  // The value of IsRequestedFactory(fac) is true, due to the above request.
  // That is why a temporary boolean "test" is used!
  TEUCHOS_TEST_FOR_EXCEPTION(IsRequestedFactory(fac) != true, Exceptions::RuntimeError, "Level::Request(ename, factory): internal logic error.");

  if (test) {
    // Call Request for factory dependencies.
    // We only do that if necessary, see comments above
    Request(*fac);
  }
}

void Level::Release(const std::string& ename, const FactoryBase* factory, const FactoryBase* requestedBy) {
  const FactoryBase* fac = GetFactory(ename, factory);
  // printf("(l=%d) [%43s] releasing  \"%20s\" generated by %10p  [actually, generated by %p (%43s)]\n",
  // levelID_, requestedBy->description().c_str(), ename.c_str(), factory, fac, fac->description().c_str());

  // Only a factory which has requested (fac,ename) is allowed to release it again.
  // Do not release data if it has not been requested by the factory "requestedBy"
  // Note: when data is released (fac,ename) depends on it often happened that some
  //       of this data has (recursively) been released too often
  if (IsRequestedBy(fac, ename, requestedBy)) {
    // In general all data (fac,ename) depends on is released when calling Get in generating factory (fac) Build method
    // Here we check the need to release the dependencies of some data that has been requested (by factory "requestedBy")
    // but the corresponding Build function of factory "fac" has never been called. Therefore the dependencies
    // have never been released. Do it now.
    if (CountRequestedFactory(fac) == 1 &&   // check if factory fac is not requested by another factory
        IsAvailableFactory(fac) == false) {  // check if Build function of factory fac has been called
      Release(*fac);
    }

    TEUCHOS_TEST_FOR_EXCEPTION(!IsKey(fac, ename), Exceptions::RuntimeError, "\"" + ename + "\" not found. Do a request first.");

    Teuchos::RCP<MueLu::VariableContainer>& v = map_[fac][ename];
    v->Release(requestedBy);

    // Remove data if no keep flag left and counter == 0
    if ((v->IsRequested() == false) && (v->GetKeepFlag() == 0)) {
      v = Teuchos::null;  // free data

      map_[fac].erase(ename);
      if (map_.count(fac) == 0)
        map_.erase(fac);
    }
  }
}

void Level::Clear() {
  // TODO: needs some love, ugly as it is
  // The ugliness is the fact that we restart both loops when we remove a single element
  bool wasRemoved;
  do {
    wasRemoved = false;
    for (TwoKeyMap::const_iterator kt = map_.begin(); kt != map_.end(); kt++) {
      const FactoryBase* factory = kt->first;

      for (SubMap::const_iterator it = kt->second.begin(); it != kt->second.end(); it++) {
        // We really want a reference here, but because later we'll need to check whether the
        // key was removed, we should copy the value
        const std::string ename = it->first;

        // We clear all the data that
        //   a) has not been requested
        //   b) is not being kept using NextRun (e.g., we clear out Final data)
        if (!IsKept(ename, factory, MueLu::NextRun)) {
          RemoveKeepFlag(ename, factory, MueLu::All);  // will delete the data if counter == 0

          // To prevent infinite looping, we need to check whether we have
          // actually removed the data. In buggy code it may happen that we
          // were unable to do that, for instance, if the data was outstanding
          // request
          if (IsKey(factory, ename)) {
            GetOStream(Errors) << "Level::Clear found Internal data inconsistency" << std::endl;
            print(GetOStream(Errors), VERB_EXTREME);

            throw Exceptions::RuntimeError("Level::Clear found Internal data inconsistency");
          }

          wasRemoved = true;
          break;
        }
      }

      if (wasRemoved)
        break;
    }

  } while (wasRemoved == true);
}

void Level::ExpertClear() {
  TwoKeyMap::const_iterator kt = map_.begin();
  while (kt != map_.end()) {
    const FactoryBase* factory = kt->first;

    SubMap::const_iterator it = kt->second.begin();
    while (it != kt->second.end()) {
      const std::string& ename = it->first;

      // obtain variable container
      Teuchos::RCP<MueLu::VariableContainer>& v = map_[factory][ename];

      if (v->GetKeepFlag() == 0 ||
          v->IsKept(MueLu::UserData) == true ||
          v->IsKept(MueLu::Final) == true) {
        it++;
        v = Teuchos::null;  // free data
        map_[factory].erase(ename);
        if (map_.count(factory) == 0) {
          break;  // last occurrence for factory has been removed. proceed with next factory
        }
      } else
        it++;
    }  // end for it

    if (map_.count(factory) == 0) {
      kt++;
      map_.erase(factory);
    } else
      kt++;
  }
}

std::string Level::description() const {
  std::ostringstream out;
  out << BaseClass::description();
  out << "{ levelID = " << levelID_ << "}";
  return out.str();
}

void Level::print(std::ostream& out, const VerbLevel verbLevel) const {
  if (!(verbLevel & Debug))
    return;

  out << "LevelID = " << GetLevelID() << std::endl;

  typedef Teuchos::TabularOutputter TTO;
  TTO outputter(out);
  outputter.pushFieldSpec("data name", TTO::STRING, TTO::LEFT, TTO::GENERAL, 20);
  outputter.pushFieldSpec("gen. factory addr.", TTO::STRING, TTO::LEFT, TTO::GENERAL, 40);
  outputter.pushFieldSpec("req", TTO::INT, TTO::LEFT, TTO::GENERAL, 3);
  outputter.pushFieldSpec("keep", TTO::STRING, TTO::LEFT, TTO::GENERAL, 5);
  outputter.pushFieldSpec("type", TTO::STRING, TTO::LEFT, TTO::GENERAL, 18);
  outputter.pushFieldSpec("data", TTO::STRING, TTO::LEFT, TTO::GENERAL, 14);
  outputter.pushFieldSpec("req'd by", TTO::STRING, TTO::LEFT, TTO::GENERAL, 20);
  outputter.outputHeader();

  for (TwoKeyMap::const_iterator kt = map_.begin(); kt != map_.end(); kt++) {
    const FactoryBase* factory = kt->first;

    for (SubMap::const_iterator it = kt->second.begin(); it != kt->second.end(); it++) {
      const std::string& ename = it->first;

      outputter.outputField(ename);  // variable name

      // NOTE: We cannot dereference the factory pointer and call factory->description() as we do not know
      // if the factory still exist (the factory pointer is a raw pointer by design). Instead, the level
      // should store the factory description internally as a string for debugging purpose (and in debug mode only).
      //         // factory name
      //         std::stringstream ss1;
      //         ss1 << (*kt)->description();
      //         outputter.outputField((ss1.str()).substr(0,30));
      if (factory == 0) {
        outputter.outputField("Null");
      } else if (factory == NoFactory::get()) {
        outputter.outputField("NoFactory");
      } else {
        std::ostringstream oss;
        oss << factory->ShortClassName() << "[" << factory->GetID() << "]";
#ifdef HAVE_MUELU_DEBUG
        oss << "(" << factory << ")";
#endif
        outputter.outputField(oss.str());
      }

      int reqcount = NumRequests(factory, ename);  // request counter
      outputter.outputField(reqcount);

      KeepType keepType = GetKeepFlag(ename, factory);
      if (keepType != 0) {
        std::stringstream ss;
        if (keepType & MueLu::UserData) {
          ss << "User";
        }
        if (keepType & MueLu::Keep) {
          ss << "Keep";
        }
        if (keepType & MueLu::Final) {
          ss << "Final";
        }
        outputter.outputField(ss.str());
      } else {
        outputter.outputField("No");
      }

      if (IsAvailable(ename, factory)) {
        std::string strType = it->second->GetTypeName();

        if (strType == "int") {
          outputter.outputField(strType);
          outputter.outputField(it->second->GetData<int>());
        } else if (strType == "double") {
          outputter.outputField(strType);
          outputter.outputField(it->second->GetData<double>());
        } else if (strType == "string") {
          outputter.outputField(strType);
          outputter.outputField(it->second->GetData<std::string>());
        } else {
          size_t npos = std::string::npos;

          if (strType.find("MueLu::Aggregates") != npos)
            outputter.outputField("Aggregates");
          else if (strType.find("MueLu::AmalgamationInfo") != npos)
            outputter.outputField("AmalgamationInfo");
          else if (strType.find("MueLu::Constraint") != npos)
            outputter.outputField("Constraint");
          else if (strType.find("MueLu::SmootherBase") != npos)
            outputter.outputField("SmootherBase");
          else if (strType.find("MueLu::SmootherPrototype") != npos)
            outputter.outputField("SmootherPrototype");
          else if (strType.find("Xpetra::Export") != npos)
            outputter.outputField("Export");
          else if (strType.find("Xpetra::Import") != npos)
            outputter.outputField("Import");
          else if (strType.find("Xpetra::Map") != npos)
            outputter.outputField("Map");
          else if (strType.find("Xpetra::Matrix") != npos)
            outputter.outputField("Matrix");
          else if (strType.find("Xpetra::MultiVector") != npos)
            outputter.outputField("Vector");
          else if (strType.find("Xpetra::Operator") != npos)
            outputter.outputField("Operator");
          else
            outputter.outputField(strType);

          outputter.outputField("available");
        }

      } else {
        outputter.outputField("unknown");
        outputter.outputField("not available");
      }

      typedef VariableContainer::request_container container_type;
      const container_type& requestedBy = it->second->Requests();
      std::ostringstream ss;
      for (container_type::const_iterator ct = requestedBy.begin(); ct != requestedBy.end(); ct++) {
        if (ct != requestedBy.begin()) ss << ",";
        ss << ct->first->ShortClassName() << "[" << ct->first->GetID() << "]";
#ifdef HAVE_MUELU_DEBUG
        ss << "(" << ct->first << ")";
#endif

        if (ct->second > 1) ss << "x" << ct->second;
      }
      outputter.outputField(ss.str());

      outputter.nextRow();
    }
  }  // for (TwoKeyMap::const_iterator kt = map_.begin(); kt != map_.end(); kt++) {
}

#if defined(HAVE_MUELU_BOOST) && defined(HAVE_MUELU_BOOST_FOR_REAL) && defined(BOOST_VERSION) && (BOOST_VERSION >= 104400)
void Level::UpdateGraph(std::map<const FactoryBase*, BoostVertex>& vindices,
                        std::map<std::pair<BoostVertex, BoostVertex>, std::string>& edges,
                        BoostProperties& dp,
                        BoostGraph& graph) const {
  size_t vind = vindices.size();

  for (TwoKeyMap::const_iterator it1 = map_.begin(); it1 != map_.end(); it1++) {
    if (vindices.find(it1->first) == vindices.end()) {
      BoostVertex boost_vertex = boost::add_vertex(graph);
      std::ostringstream oss;
      oss << it1->first->ShortClassName() << "[" << it1->first->GetID() << "]";
      boost::put("label", dp, boost_vertex, oss.str());
      vindices[it1->first] = vind++;
    }

    for (SubMap::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); it2++) {
      const VariableContainer::request_container& requests = it2->second->Requests();
      for (VariableContainer::request_container::const_iterator rit = requests.begin(); rit != requests.end(); rit++) {
        if (vindices.find(rit->first) == vindices.end()) {
          // requested by factory which is unknown
          BoostVertex boost_vertex = boost::add_vertex(graph);
          std::ostringstream oss;
          oss << rit->first->ShortClassName() << "[" << rit->first->GetID() << "]";
          boost::put("label", dp, boost_vertex, oss.str());
          vindices[rit->first] = vind++;
        }

        edges[std::pair<BoostVertex, BoostVertex>(vindices[rit->first], vindices[it1->first])] = it2->first;
      }
    }
  }
}
#endif

// JG Note: should the option IgnoreUserData() moved to the Factory interface or on the specific factories that are using this option? It would simplify the level class.
const FactoryBase* Level::GetFactory(const std::string& ename, const FactoryBase* factory) const {
  if (factory != NULL)
    return factory;

  // If IgnoreUserData == false and if variable "ename" generated by NoFactory is provided by the user (MueLu::UserData),
  // use user-provided data by default without querying the FactoryManager.
  // When FactoryManager == null, we consider that IgnoreUserData == false.
  if ((factoryManager_ == Teuchos::null || factoryManager_->IgnoreUserData() == false) &&
      (IsAvailable(ename, NoFactory::get()) && IsKept(ename, NoFactory::get(), MueLu::UserData))) {
    return NoFactory::get();
  }

  // Query factory manager
  TEUCHOS_TEST_FOR_EXCEPTION(factoryManager_ == null, Exceptions::RuntimeError, "MueLu::Level(" << levelID_ << ")::GetFactory(" << ename << ", " << factory << "): No FactoryManager");
  const FactoryBase* fac = factoryManager_->GetFactory(ename).get();
  TEUCHOS_TEST_FOR_EXCEPTION(fac == NULL, Exceptions::RuntimeError, "MueLu::Level(" << levelID_ << ")::GetFactory(" << ename << ", " << factory << "): Default factory returned by FactoryManager cannot be NULL");
  return fac;
}

Level::RequestMode Level::requestMode_ = UNDEF;

}  // namespace MueLu

// TODO: Caps should not matter