C++ Level::Release方法代码示例

  TEUCHOS_UNIT_TEST(CoarseMap, NonStandardCaseA )
  {
    out << "version: " << MueLu::Version() << std::endl;
    Level myLevel;
    myLevel.SetLevelID(0);
    RCP<Matrix> A = TestHelpers::TestFactory<SC, LO, GO, NO, LMO>::Build1DPoisson(15);
    myLevel.Set("A", A);

    // build dummy aggregate structure
    Teuchos::RCP<Aggregates> aggs = Teuchos::rcp(new Aggregates(A->getRowMap()));
    aggs->SetNumAggregates(10); // set (local!) number of aggregates
    myLevel.Set("Aggregates", aggs);

    // build dummy nullspace vector
    Teuchos::RCP<MultiVector> nsp = MultiVectorFactory::Build(A->getRowMap(),1);
    nsp->putScalar(1.0);
    myLevel.Set("Nullspace", nsp);

    RCP<CoarseMapFactory> myCMF = Teuchos::rcp(new CoarseMapFactory());
    myLevel.Request("CoarseMap",myCMF.get());
    myCMF->SetParameter("Domain GID offsets",Teuchos::ParameterEntry(std::string("{100,50}")));
    myCMF->SetFactory("Aggregates",MueLu::NoFactory::getRCP());
    myCMF->SetFactory("Nullspace",MueLu::NoFactory::getRCP());
    myCMF->Build(myLevel);
    Teuchos::RCP<const Map> myCoarseMap = myLevel.Get<Teuchos::RCP<const Map> >("CoarseMap",myCMF.get());

    TEST_EQUALITY(myCoarseMap->getMinAllGlobalIndex() == 100, true);
    TEST_EQUALITY(myCoarseMap->getMaxLocalIndex()==9,true);

    myLevel.Release("CoarseMap",myCMF.get());
    myLevel.SetLevelID(1);
    myLevel.Request("CoarseMap",myCMF.get());
    myCMF->SetParameter("Domain GID offsets",Teuchos::ParameterEntry(std::string("{100,50}")));
    myCMF->SetFactory("Aggregates",MueLu::NoFactory::getRCP());
    myCMF->SetFactory("Nullspace",MueLu::NoFactory::getRCP());
    myCMF->Build(myLevel);
    myCoarseMap = myLevel.Get<Teuchos::RCP<const Map> >("CoarseMap",myCMF.get());

    TEST_EQUALITY(myCoarseMap->getMinAllGlobalIndex() == 50, true);
    TEST_EQUALITY(myCoarseMap->getMaxLocalIndex()==9,true);

    myLevel.Release("CoarseMap",myCMF.get());
    myLevel.SetLevelID(2);
    myLevel.Request("CoarseMap",myCMF.get());
    myCMF->SetParameter("Domain GID offsets",Teuchos::ParameterEntry(std::string("{100,50}")));
    myCMF->SetFactory("Aggregates",MueLu::NoFactory::getRCP());
    myCMF->SetFactory("Nullspace",MueLu::NoFactory::getRCP());
    myCMF->Build(myLevel);
    myCoarseMap = myLevel.Get<Teuchos::RCP<const Map> >("CoarseMap",myCMF.get());

    TEST_EQUALITY(myCoarseMap->getMinAllGlobalIndex() == 0, true);
    TEST_EQUALITY(myCoarseMap->getMaxLocalIndex()==9,true);
  }

  void TogglePFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level& fineLevel, Level &coarseLevel) const {
    FactoryMonitor m(*this, "Prolongator toggle", coarseLevel);
    std::ostringstream levelstr;
    levelstr << coarseLevel.GetLevelID();

    typedef typename Teuchos::ScalarTraits<SC>::magnitudeType Magnitude;

    TEUCHOS_TEST_FOR_EXCEPTION(nspFacts_.size() != prolongatorFacts_.size(), Exceptions::RuntimeError, "MueLu::TogglePFactory::Build: The number of provided prolongator factories and coarse nullspace factories must be identical.");
    TEUCHOS_TEST_FOR_EXCEPTION(nspFacts_.size() != 2, Exceptions::RuntimeError, "MueLu::TogglePFactory::Build: TogglePFactory needs two different transfer operator strategies for toggling."); // TODO adapt this/weaken this as soon as other toggling strategies are introduced.

    // decision routine which prolongator factory to be used
    int nProlongatorFactory = 0; // default behavior: use first prolongator in list

    // extract user parameters
    const Teuchos::ParameterList & pL = GetParameterList();
    std::string mode = Teuchos::as<std::string>(pL.get<std::string>("toggle: mode"));
    int semicoarsen_levels = Teuchos::as<int>(pL.get<int>("semicoarsen: number of levels"));

    TEUCHOS_TEST_FOR_EXCEPTION(mode!="semicoarsen", Exceptions::RuntimeError, "MueLu::TogglePFactory::Build: The 'toggle: mode' parameter must be set to 'semicoarsen'. No other mode supported, yet.");

    LO NumZDir = -1;
    if(fineLevel.IsAvailable("NumZLayers", NoFactory::get())) {
      NumZDir = fineLevel.Get<LO>("NumZLayers", NoFactory::get()); //obtain info
      GetOStream(Runtime1) << "Number of layers for semicoarsening: " << NumZDir << std::endl;
    }

    // Make a decision which prolongator to be used.
    if(fineLevel.GetLevelID() >= semicoarsen_levels || NumZDir == 1) {
      nProlongatorFactory = 1;
    } else {
      nProlongatorFactory = 0;
    }

    RCP<Matrix> P = Teuchos::null;
    RCP<MultiVector> coarseNullspace = Teuchos::null;

    // call Build for selected transfer operator
    GetOStream(Runtime0) << "TogglePFactory: call transfer factory: " << (prolongatorFacts_[nProlongatorFactory])->description() << std::endl;
    prolongatorFacts_[nProlongatorFactory]->CallBuild(coarseLevel);
    P = coarseLevel.Get< RCP<Matrix> >("P", (prolongatorFacts_[nProlongatorFactory]).get());
    coarseNullspace = coarseLevel.Get< RCP<MultiVector> >("Nullspace", (nspFacts_[nProlongatorFactory]).get());

    // Release dependencies of all prolongator and coarse level null spaces
    for(size_t t=0; t<nspFacts_.size(); ++t) {
      coarseLevel.Release(*(prolongatorFacts_[t]));
      coarseLevel.Release(*(nspFacts_[t]));
    }

    // store prolongator with this factory identification.
    Set(coarseLevel, "P", P);
    Set(coarseLevel, "Nullspace", coarseNullspace);

  } //Build()

  void ToggleCoordinatesTransferFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level & fineLevel, Level &coarseLevel) const {
    FactoryMonitor m(*this, "Coordinate transfer toggle", coarseLevel);

    typedef Xpetra::MultiVector<double,LO,GO,NO> xdMV;

    TEUCHOS_TEST_FOR_EXCEPTION(coordFacts_.size() != 2, Exceptions::RuntimeError, "MueLu::TogglePFactory::Build: ToggleCoordinatesTransferFactory needs two different transfer operator strategies for toggling.");

    int chosenP = Get< int >      (coarseLevel, "Chosen P");
    GetOStream(Runtime1) << "Transfer Coordinates" << chosenP << " to coarse level" << std::endl;
    RCP<xdMV> coarseCoords = coarseLevel.Get< RCP<xdMV> >("Coordinates",(coordFacts_[chosenP]).get());
    Set(coarseLevel, "Coordinates", coarseCoords);

    // loop through all coord facts and check whether the coarse coordinates are available.
    // This is the coarse coordinate transfer factory which belongs to the execution path
    // chosen by the TogglePFactory
    /*RCP<xdMV> coarseCoords = Teuchos::null;
    for(size_t t=0; t<coordFacts_.size(); ++t) {
      bool bIsAv = coarseLevel.IsAvailable("Coordinates",(coordFacts_[t]).get());
      std::cout << "Coordinates generated by " << (coordFacts_[t]).get() << " available? " << bIsAv << std::endl;
      if ( coarseLevel.IsAvailable("Coordinates",(coordFacts_[t]).get()) ) {
        GetOStream(Runtime1) << "Choose factory " << t << " (" << (coordFacts_[t]).get() << ")" << std::endl;
        coarseCoords = coarseLevel.Get< RCP<xdMV> >("Coordinates",(coordFacts_[t]).get());
        Set(coarseLevel, "Coordinates", coarseCoords);
      }
    }*/

    // Release dependencies of all coordinate transfer factories
    for(size_t t=0; t<coordFacts_.size(); ++t) {
      coarseLevel.Release(*(coordFacts_[t]));
    }

    //TODO: exception if coarseCoords == Teuchos::null
  }

  TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(ThresholdAFilterFactory, Basic, Scalar, LocalOrdinal, GlobalOrdinal, Node)
  {
#   include <MueLu_UseShortNames.hpp>
    MUELU_TESTING_SET_OSTREAM;
    MUELU_TESTING_LIMIT_EPETRA_SCOPE(Scalar,GlobalOrdinal,Node);
    out << "version: " << MueLu::Version() << std::endl;

    Level aLevel;
    TestHelpers::TestFactory<SC, LO, GO, NO>::createSingleLevelHierarchy(aLevel);

    RCP<Matrix> A = TestHelpers::TestFactory<SC, LO, GO, NO>::Build1DPoisson(20); //can be an empty operator

    RCP<ThresholdAFilterFactory> AfilterFactory0 = rcp(new ThresholdAFilterFactory("A",0.1)); // keep all
    RCP<ThresholdAFilterFactory> AfilterFactory1 = rcp(new ThresholdAFilterFactory("A",1.1)); // keep only diagonal
    RCP<ThresholdAFilterFactory> AfilterFactory2 = rcp(new ThresholdAFilterFactory("A",3));   // keep only diagonal

    aLevel.Set("A",A);

    aLevel.Request("A",AfilterFactory0.get());
    AfilterFactory0->Build(aLevel);
    TEST_EQUALITY(aLevel.IsAvailable("A",AfilterFactory0.get()), true);
    RCP<Matrix> A0 = aLevel.Get< RCP<Matrix> >("A",AfilterFactory0.get());
    aLevel.Release("A",AfilterFactory0.get());
    TEST_EQUALITY(aLevel.IsAvailable("A",AfilterFactory0.get()), false);
    TEST_EQUALITY(A0->getNodeNumEntries(), A->getNodeNumEntries());
    TEST_EQUALITY(A0->getGlobalNumEntries(), A->getGlobalNumEntries());

    aLevel.Request("A",AfilterFactory1.get());
    AfilterFactory1->Build(aLevel);
    TEST_EQUALITY(aLevel.IsAvailable("A",AfilterFactory1.get()), true);
    RCP<Matrix> A1 = aLevel.Get< RCP<Matrix> >("A",AfilterFactory1.get());
    aLevel.Release("A",AfilterFactory1.get());
    TEST_EQUALITY(aLevel.IsAvailable("A",AfilterFactory1.get()), false);
    TEST_EQUALITY(A1->getGlobalNumEntries(), A1->getGlobalNumRows());

    aLevel.Request("A",AfilterFactory2.get());
    AfilterFactory2->Build(aLevel);
    TEST_EQUALITY(aLevel.IsAvailable("A",AfilterFactory2.get()), true);
    RCP<Matrix> A2 = aLevel.Get< RCP<Matrix> >("A",AfilterFactory2.get());
    aLevel.Release("A",AfilterFactory2.get());
    TEST_EQUALITY(aLevel.IsAvailable("A",AfilterFactory2.get()), false);
    TEST_EQUALITY(A2->getGlobalNumEntries(), A2->getGlobalNumRows());


  }

  Teuchos::RCP<MueLu::Aggregates_kokkos<LocalOrdinal, GlobalOrdinal, Node>>
  gimmeUncoupledAggregates(const Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>>& A,
                           Teuchos::RCP<MueLu::AmalgamationInfo<LocalOrdinal, GlobalOrdinal, Node>>& amalgInfo,
                           bool bPhase1 = true, bool bPhase2a = true, bool bPhase2b = true, bool bPhase3 = true) {
#   include "MueLu_UseShortNames.hpp"

    Level level;
    TestHelpers_kokkos::TestFactory<SC,LO,GO,NO>::createSingleLevelHierarchy(level);
    level.Set("A", A);

    RCP<AmalgamationFactory>        amalgFact = rcp(new AmalgamationFactory());
    RCP<CoalesceDropFactory_kokkos> dropFact  = rcp(new CoalesceDropFactory_kokkos());
    dropFact->SetFactory("UnAmalgamationInfo", amalgFact);

    using Teuchos::ParameterEntry;

    // Setup aggregation factory (use default factory for graph)
    RCP<UncoupledAggregationFactory_kokkos> aggFact = rcp(new UncoupledAggregationFactory_kokkos());
    aggFact->SetFactory("Graph", dropFact);
    aggFact->SetParameter("aggregation: max agg size",           ParameterEntry(3));
    aggFact->SetParameter("aggregation: min agg size",           ParameterEntry(3));
    aggFact->SetParameter("aggregation: max selected neighbors", ParameterEntry(0));
    aggFact->SetParameter("aggregation: ordering",               ParameterEntry(std::string("natural")));
    aggFact->SetParameter("aggregation: enable phase 1",         ParameterEntry(bPhase1));
    aggFact->SetParameter("aggregation: enable phase 2a",        ParameterEntry(bPhase2a));
    aggFact->SetParameter("aggregation: enable phase 2b",        ParameterEntry(bPhase2b));
    aggFact->SetParameter("aggregation: enable phase 3",         ParameterEntry(bPhase3));

    level.Request("Aggregates",         aggFact.get());
    level.Request("UnAmalgamationInfo", amalgFact.get());

    level.Request(*aggFact);
    aggFact->Build(level);

    auto aggregates = level.Get<RCP<Aggregates_kokkos> >("Aggregates", aggFact.get());
    amalgInfo = level.Get<RCP<AmalgamationInfo> >("UnAmalgamationInfo", amalgFact.get());

    level.Release("UnAmalgamationInfo", amalgFact.get());
    level.Release("Aggregates",         aggFact.get());

    return aggregates;
  }

  void BlockedRAPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level &fineLevel, Level &coarseLevel) const {  //FIXME make fineLevel const!!
    FactoryMonitor m(*this, "Computing Ac (block)", coarseLevel);

    const Teuchos::ParameterList& pL = GetParameterList();

    RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,   "A");
    RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel, "P");

    RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A);
    RCP<BlockedCrsMatrix> bP = rcp_dynamic_cast<BlockedCrsMatrix>(P);
    TEUCHOS_TEST_FOR_EXCEPTION(bA.is_null() || bP.is_null(), Exceptions::BadCast, "Matrices R, A and P must be of type BlockedCrsMatrix.");

    RCP<BlockedCrsMatrix> bAP;
    RCP<BlockedCrsMatrix> bAc;
    {
      SubFactoryMonitor subM(*this, "MxM: A x P", coarseLevel);

      // Triple matrix product for BlockedCrsMatrixClass
      TEUCHOS_TEST_FOR_EXCEPTION((bA->Cols() != bP->Rows()), Exceptions::BadCast,
                               "Block matrix dimensions do not match: "
                               "A is " << bA->Rows() << "x" << bA->Cols() <<
                               "P is " << bP->Rows() << "x" << bP->Cols());

      bAP = Utils::TwoMatrixMultiplyBlock(*bA, false, *bP,  false, GetOStream(Statistics2), true, true);
    }


    // If we do not modify matrix later, allow optimization of storage.
    // This is necessary for new faster Epetra MM kernels.
    bool doOptimizeStorage = !checkAc_;

    const bool doTranspose    = true;
    const bool doFillComplete = true;
    if (pL.get<bool>("transpose: use implicit") == true) {
      SubFactoryMonitor m2(*this, "MxM: P' x (AP) (implicit)", coarseLevel);
      bAc = Utils::TwoMatrixMultiplyBlock(*bP,  doTranspose, *bAP, !doTranspose, GetOStream(Statistics2), doFillComplete, doOptimizeStorage);

    } else {
      RCP<Matrix>           R  = Get< RCP<Matrix> >(coarseLevel, "R");
      RCP<BlockedCrsMatrix> bR = rcp_dynamic_cast<BlockedCrsMatrix>(R);
      TEUCHOS_TEST_FOR_EXCEPTION(bR.is_null(), Exceptions::BadCast, "Matrix R must be of type BlockedCrsMatrix.");

      TEUCHOS_TEST_FOR_EXCEPTION(bA->Rows() != bR->Cols(), Exceptions::BadCast,
                                 "Block matrix dimensions do not match: "
                                 "R is " << bR->Rows() << "x" << bR->Cols() <<
                                 "A is " << bA->Rows() << "x" << bA->Cols());

      SubFactoryMonitor m2(*this, "MxM: R x (AP) (explicit)", coarseLevel);
      bAc = Utils::TwoMatrixMultiplyBlock(*bR, !doTranspose, *bAP, !doTranspose, GetOStream(Statistics2), doFillComplete, doOptimizeStorage);
    }


    if (checkAc_)
      CheckMainDiagonal(bAc);

    GetOStream(Statistics1) << PerfUtils::PrintMatrixInfo(*bAc, "Ac (blocked)");

    // static int run = 1;
    // RCP<CrsMatrixWrap> A11 = rcp(new CrsMatrixWrap(bAc->getMatrix(0,0)));
    // Utils::Write(toString(run) + "_A_11.mm", *A11);
    // if (!bAc->getMatrix(1,1).is_null()) {
      // RCP<CrsMatrixWrap> A22 = rcp(new CrsMatrixWrap(bAc->getMatrix(1,1)));
      // Utils::Write(toString(run) + "_A_22.mm", *A22);
    // }
    // RCP<CrsMatrixWrap> Am = rcp(new CrsMatrixWrap(bAc->Merge()));
    // Utils::Write(toString(run) + "_A.mm", *Am);
    // run++;

    Set<RCP <Matrix> >(coarseLevel, "A", bAc);

    if (transferFacts_.begin() != transferFacts_.end()) {
      SubFactoryMonitor m1(*this, "Projections", coarseLevel);

      // call Build of all user-given transfer factories
      for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
        RCP<const FactoryBase> fac = *it;
        GetOStream(Runtime0) << "BlockRAPFactory: call transfer factory: " << fac->description() << std::endl;
        fac->CallBuild(coarseLevel);
        // AP (11/11/13): I am not sure exactly why we need to call Release, but we do need it to get rid
        // of dangling data for CoordinatesTransferFactory
        coarseLevel.Release(*fac);
      }
    }
  }

//.........这里部分代码省略.........
      RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,   "A");
      RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel, "P"), AP, Ac;

      // Reuse pattern if available (multiple solve)
      if (coarseLevel.IsAvailable("AP Pattern", this)) {
        GetOStream(Runtime0) << "Ac: Using previous AP pattern" << std::endl;

        AP = Get< RCP<Matrix> >(coarseLevel, "AP Pattern");
      }

      {
        SubFactoryMonitor subM(*this, "MxM: A x P", coarseLevel);

        AP = Utils::Multiply(*A, false, *P, false, AP, GetOStream(Statistics2),true,true,std::string("MueLu::A*P-")+levelstr.str());
      }
      if (pL.get<bool>("Keep AP Pattern"))
        Set(coarseLevel, "AP Pattern", AP);

      // Reuse coarse matrix memory if available (multiple solve)
      if (coarseLevel.IsAvailable("RAP Pattern", this)) {
        GetOStream(Runtime0) << "Ac: Using previous RAP pattern" << std::endl;

        Ac = Get< RCP<Matrix> >(coarseLevel, "RAP Pattern");
        // Some eigenvalue may have been cached with the matrix in the previous run.
        // As the matrix values will be updated, we need to reset the eigenvalue.
        Ac->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
      }

      // If we do not modify matrix later, allow optimization of storage.
      // This is necessary for new faster Epetra MM kernels.
      bool doOptimizeStorage = !pL.get<bool>("RepairMainDiagonal");

      const bool doTranspose    = true;
      const bool doFillComplete = true;
      if (pL.get<bool>("transpose: use implicit") == true) {
        SubFactoryMonitor m2(*this, "MxM: P' x (AP) (implicit)", coarseLevel);
        Ac = Utils::Multiply(*P,  doTranspose, *AP, !doTranspose, Ac, GetOStream(Statistics2), doFillComplete, doOptimizeStorage,std::string("MueLu::R*(AP)-implicit-")+levelstr.str());

      } else {
        RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel, "R");

        SubFactoryMonitor m2(*this, "MxM: R x (AP) (explicit)", coarseLevel);
        Ac = Utils::Multiply(*R, !doTranspose, *AP, !doTranspose, Ac, GetOStream(Statistics2), doFillComplete, doOptimizeStorage,std::string("MueLu::R*(AP)-explicit-")+levelstr.str());
      }

      CheckRepairMainDiagonal(Ac);

      if (IsPrint(Statistics1)) {
        RCP<ParameterList> params = rcp(new ParameterList());;
        params->set("printLoadBalancingInfo", true);
        params->set("printCommInfo",          true);
        GetOStream(Statistics1) << PerfUtils::PrintMatrixInfo(*Ac, "Ac", params);
      }

      Set(coarseLevel, "A",           Ac);
      if (pL.get<bool>("Keep RAP Pattern"))
        Set(coarseLevel, "RAP Pattern", Ac);
    }

    if (transferFacts_.begin() != transferFacts_.end()) {
      SubFactoryMonitor m(*this, "Projections", coarseLevel);

      // call Build of all user-given transfer factories
      for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
        RCP<const FactoryBase> fac = *it;
        GetOStream(Runtime0) << "RAPFactory: call transfer factory: " << fac->description() << std::endl;
        fac->CallBuild(coarseLevel);
        // Coordinates transfer is marginally different from all other operations
        // because it is *optional*, and not required. For instance, we may need
        // coordinates only on level 4 if we start repartitioning from that level,
        // but we don't need them on level 1,2,3. As our current Hierarchy setup
        // assumes propagation of dependencies only through three levels, this
        // means that we need to rely on other methods to propagate optional data.
        //
        // The method currently used is through RAP transfer factories, which are
        // simply factories which are called at the end of RAP with a single goal:
        // transfer some fine data to coarser level. Because these factories are
        // kind of outside of the mainline factories, they behave different. In
        // particular, we call their Build method explicitly, rather than through
        // Get calls. This difference is significant, as the Get call is smart
        // enough to know when to release all factory dependencies, and Build is
        // dumb. This led to the following CoordinatesTransferFactory sequence:
        // 1. Request level 0
        // 2. Request level 1
        // 3. Request level 0
        // 4. Release level 0
        // 5. Release level 1
        //
        // The problem is missing "6. Release level 0". Because it was missing,
        // we had outstanding request on "Coordinates", "Aggregates" and
        // "CoarseMap" on level 0.
        //
        // This was fixed by explicitly calling Release on transfer factories in
        // RAPFactory. I am still unsure how exactly it works, but now we have
        // clear data requests for all levels.
        coarseLevel.Release(*fac);
      }
    }

  }

void BlockedRAPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node, LocalMatOps>::Build(Level &fineLevel, Level &coarseLevel) const {  //FIXME make fineLevel const!!
    typedef Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, LocalMatOps> BlockedCrsMatrixClass; // TODO move me

    FactoryMonitor m(*this, "Computing Ac (block)", coarseLevel);

    //
    // Inputs: R, A, P
    //

    RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel, "R");
    RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,   "A");
    RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel, "P");

    //
    // Dynamic casts
    //

    RCP<BlockedCrsMatrixClass> bR, bA, bP;

    try {
        /* using rcp_dynamic_cast with throw_on_fail = true */
        bR = Teuchos::rcp_dynamic_cast<BlockedCrsMatrixClass>(R, true);
        bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrixClass>(A, true);
        bP = Teuchos::rcp_dynamic_cast<BlockedCrsMatrixClass>(P, true);
    } catch(std::bad_cast e) {
        TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::BadCast, "MueLu::BlockedRAPFactory::Build(): matrices R, A and P must be of type BlockedCrsMatrix. " << e.what());
    }

    /*Utils::Write( "A00.m", CrsMatrixWrap(bA->getMatrix(0,0)) );
    Utils::Write( "A11.m", CrsMatrixWrap(bA->getMatrix(1,1)) );
    Utils::Write( "A01.m", CrsMatrixWrap(bA->getMatrix(0,1)) );
    Utils::Write( "A10.m", CrsMatrixWrap(bA->getMatrix(1,0)) );

    Utils::Write( "P00.m", CrsMatrixWrap(bP->getMatrix(0,0)) );
    Utils::Write( "P11.m", CrsMatrixWrap(bP->getMatrix(1,1)) );*/

    //
    // Build Ac = RAP
    //

    // Triple matrix product for BlockedCrsMatrixClass
    TEUCHOS_TEST_FOR_EXCEPTION((bA->Cols() != bP->Rows()) || (bA->Rows() != bR->Cols()), Exceptions::BadCast, "MueLu::BlockedRAPFactory::Build(): block matrix dimensions do not match.");
    RCP<BlockedCrsMatrixClass> bAP = Utils::TwoMatrixMultiplyBlock(*bA, false, *bP,  false, true, true);
    RCP<BlockedCrsMatrixClass> bAc = Utils::TwoMatrixMultiplyBlock(*bR, false, *bAP, false, true, true);

    if (checkAc_)
        CheckMainDiagonal(bAc);

    GetOStream(Statistics1, 0) << Utils::PrintMatrixInfo(*bAc, "Ac (blocked)");

    Set<RCP <Matrix> >(coarseLevel, "A", bAc);

    if (transferFacts_.begin() != transferFacts_.end()) {
        SubFactoryMonitor m1(*this, "Projections", coarseLevel);

        // call Build of all user-given transfer factories
        for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
            RCP<const FactoryBase> fac = *it;
            GetOStream(Runtime0, 0) << "BlockRAPFactory: call transfer factory: " << fac->description() << std::endl;
            fac->CallBuild(coarseLevel);
            // AP (11/11/13): I am not sure exactly why we need to call Release, but we do need it to get rid
            // of dangling data for CoordinatesTransferFactory
            coarseLevel.Release(*fac);
        }
    }
}