C++ Teuchos::rcpFromRef方法代码示例

    /*! \brief This routine takes the Epetra_MultiVector \c x and applies the operator
      to it resulting in the Epetra_MultiVector \c y, which is returned.
      \note It is expected that any problem with applying this operator to \c x will be
      indicated by an std::exception being thrown.
    */
    void Apply ( const Epetra_MultiVector& x, Epetra_MultiVector& y, ETrans trans=NOTRANS ) const {
      TEUCHOS_TEST_FOR_EXCEPTION(trans!=NOTRANS, XpetraOpFailure,
                         "Belos::MueLuTpetraOp::Apply, transpose mode != NOTRANS not supported.");

      Epetra_MultiVector & temp_x = const_cast<Epetra_MultiVector &>(x);

      const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> tX(rcpFromRef(temp_x));
      Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node>       tY(rcpFromRef(y));

      //FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate().
      tY.putScalar(0.0);

      Op_->apply(tX,tY);
    }

void DefaultBlockedLinearOp<Scalar>::describe(
  Teuchos::FancyOStream &out_arg
  ,const Teuchos::EVerbosityLevel verbLevel
  ) const
{
  typedef Teuchos::ScalarTraits<Scalar> ST;
  using Teuchos::rcpFromRef;
  using Teuchos::FancyOStream;
  using Teuchos::OSTab;
  assertBlockFillIsActive(false);
  RCP<FancyOStream> out = rcpFromRef(out_arg);
  OSTab tab1(out);
  switch(verbLevel) {
    case Teuchos::VERB_DEFAULT:
    case Teuchos::VERB_LOW:
      *out << this->description() << std::endl;
      break;
    case Teuchos::VERB_MEDIUM:
    case Teuchos::VERB_HIGH:
    case Teuchos::VERB_EXTREME:
    {
      *out
        << Teuchos::Describable::description() << "{"
        << "rangeDim=" << this->range()->dim()
        << ",domainDim=" << this->domain()->dim()
        << ",numRowBlocks=" << numRowBlocks_
        << ",numColBlocks=" << numColBlocks_
        << "}\n";
      OSTab tab2(out);
      *out
        << "Constituent LinearOpBase objects for M = [ Op[0,0] ..."
        << " ; ... ; ... Op[numRowBlocks-1,numColBlocks-1] ]:\n";
      tab2.incrTab();
      for( int i = 0; i < numRowBlocks_; ++i ) {
        for( int j = 0; j < numColBlocks_; ++j ) {
          *out << "Op["<<i<<","<<j<<"] = ";
          RCP<const LinearOpBase<Scalar> >
            block_i_j = getBlock(i,j);
          if(block_i_j.get())
            *out << Teuchos::describe(*getBlock(i,j),verbLevel);
          else
            *out << "NULL\n";
        }
      }
      break;
    }
    default:
      TEST_FOR_EXCEPT(true); // Should never get here!
  }
}

    /*! \brief This routine takes the Tpetra::MultiVector \c x and applies the operator
      to it resulting in the Tpetra::MultiVector \c y, which is returned.
      \note It is expected that any problem with applying this operator to \c x will be
      indicated by an std::exception being thrown.
    */
    void Apply ( const Epetra_MultiVector& x, Epetra_MultiVector& y, ETrans trans=NOTRANS ) const {

      TEUCHOS_TEST_FOR_EXCEPTION(trans!=NOTRANS, MueLuOpFailure,
                         "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");

      Epetra_MultiVector & temp_x = const_cast<Epetra_MultiVector &>(x);

      const Xpetra::EpetraMultiVector tX(rcpFromRef(temp_x));
      Xpetra::EpetraMultiVector       tY(rcpFromRef(y));

      //FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate().
      tY.putScalar(0.0);

      Hierarchy_->Iterate( tX, 1, tY , true);

    }

  Teuchos::RCP<const Teuchos::ParameterList>
  MinresSolMgr<ScalarType, MV, OP>::defaultParameters()
  {
    using Teuchos::ParameterList;
    using Teuchos::parameterList;
    using Teuchos::RCP;
    using Teuchos::rcp;
    using Teuchos::rcpFromRef;
    using Teuchos::EnhancedNumberValidator;
    typedef MagnitudeType MT;
    typedef Teuchos::ScalarTraits<MT> MST;

    // List of parameters accepted by MINRES, and their default values.
    RCP<ParameterList> pl = parameterList ("MINRES");

    pl->set ("Convergence Tolerance", MST::squareroot (MST::eps()),
	     "Relative residual tolerance that needs to be achieved by "
	     "the iterative solver, in order for the linear system to be "
	     "declared converged.",
	     rcp (new EnhancedNumberValidator<MT> (MST::zero(), MST::rmax())));
    pl->set ("Maximum Iterations", static_cast<int>(1000),
	     "Maximum number of iterations allowed for each right-hand "
	     "side solved.",
	     rcp (new EnhancedNumberValidator<int> (0, INT_MAX)));
    pl->set ("Block Size", static_cast<int>(1),
	     "Number of vectors in each block.  WARNING: The current "
	     "implementation of MINRES only accepts a block size of 1, "
	     "since it can only solve for 1 right-hand side at a time.",
	     rcp (new EnhancedNumberValidator<int> (1, 1)));
    pl->set ("Verbosity", (int) Belos::Errors,
	     "The type(s) of solver information that should "
	     "be written to the output stream.");
    pl->set ("Output Style", (int) Belos::General,
	     "What style is used for the solver information written "
	     "to the output stream.");
    pl->set ("Output Frequency", static_cast<int>(-1),
	     "How often (in terms of number of iterations) intermediate "
	     "convergence information should be written to the output stream."
	     "  -1 means never.");
    pl->set ("Output Stream", rcpFromRef(std::cout),
	     "A reference-counted pointer to the output stream where all "
	     "solver output is sent.  The output stream defaults to stdout.");
    pl->set ("Timer Label", std::string("Belos"),
	     "The string to use as a prefix for the timer labels.");
    return pl;
  }

  void
  DistObjectKA<Packet,LocalOrdinal,GlobalOrdinal,Node>::
  describe (Teuchos::FancyOStream &out,
            const Teuchos::EVerbosityLevel verbLevel) const
  {
    using Teuchos::rcpFromRef;
    using std::endl;

    const Teuchos::EVerbosityLevel vl = (verbLevel == Teuchos::VERB_DEFAULT) ?
      Teuchos::VERB_LOW : verbLevel;

    if (vl != Teuchos::VERB_NONE) {
      out << this->description () << endl;
      Teuchos::OSTab tab (rcpFromRef (out));
      out << "Export buffer size (in packets): " << exports_.size() << endl
          << "Import buffer size (in packets): " << imports_.size() << endl
          << "Map over which this object is distributed:" << endl;
      map_->describe (out, vl);
    }
  }

    /*! \brief This routine takes the Tpetra::MultiVector \c x and applies the operator
      to it resulting in the Tpetra::MultiVector \c y, which is returned.
      \note It is expected that any problem with applying this operator to \c x will be
      indicated by an std::exception being thrown.
    */
    void Apply(const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x, Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& y, ETrans trans = NOTRANS ) const {

      TEUCHOS_TEST_FOR_EXCEPTION(trans!=NOTRANS, MueLuOpFailure,
                         "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");

      //FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate(), but it matters for AMGX
      y.putScalar(0.0);

#ifdef HAVE_MUELU_AMGX
      if (!AMGX_.is_null())
        AMGX_->apply(x, y);
#endif

      if (!Hierarchy_.is_null()) {
        Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & temp_x = const_cast<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &>(x);

        const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX(rcpFromRef(temp_x));
        Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tY(rcpFromRef(y));
        Hierarchy_->Iterate(tX, tY, 1, true);
      }
    }

bool tNeumannSeries::test_scaledOp(int verbosity,std::ostream & os)
{
   bool status = false;
   bool allPassed = true;

   // perform actual test
   Thyra::LinearOpTester<double> tester;
   tester.show_all_tests(true);
   std::stringstream ss;
   Teuchos::FancyOStream fos(rcpFromRef(ss),"      |||");

   {
      // build library parameter list
      Teuchos::RCP<Teuchos::ParameterList> pl = buildLibPL(2,"Diagonal");
      if(verbosity>=10) {
         os << "   tNeumannSeries::test_scaledOp :"
            << " printing library parameter list" << std::endl;
         pl->print(os);
      }

      RCP<Teko::InverseLibrary> invLib = Teko::InverseLibrary::buildFromParameterList(*pl);
      RCP<Teko::InverseFactory> neumann = invLib->getInverseFactory("Neumann");
      RCP<Teko::InverseFactory> direct = invLib->getInverseFactory("Amesos");

      Teko::LinearOp op = buildExampleOp(2,*GetComm());
   
      Teko::LinearOp neuInv = Teko::buildInverse(*neumann,op);
      Teko::LinearOp dirInv = Teko::buildInverse(*direct,op);

      const bool result = tester.compare( *neuInv, *dirInv, Teuchos::ptrFromRef(fos) );
      TEST_ASSERT(result,
             std::endl << "   tNeumannSeries::test_scaledOp "
             << ": Comparing factory generated operator to correct operator");
      if(not result || verbosity>=10) 
         os << ss.str(); 
   }
 
   return allPassed;
}

void BlockImplicitLinearOp::applyImpl(
  const Thyra::EOpTransp M_trans,
  const Thyra::MultiVectorBase<double> & x,
  const Teuchos::Ptr<Thyra::MultiVectorBase<double> > & y,
  const double alpha,
  const double beta
  ) const
{
   TEST_FOR_EXCEPTION(M_trans!=Thyra::NOTRANS, std::runtime_error,
     "Linear operators of inherited type BlockImplicitLinearOp "
     "cannot handle conjugation (yet!)");

   // cast source vector
   RCP<const ProductMultiVectorBase<double> > src =
     rcp_dynamic_cast<const ProductMultiVectorBase<double> >(rcpFromRef(x));
   BlockedMultiVector srcX = rcp_const_cast<ProductMultiVectorBase<double> >(src);

   // cast destination vector
   BlockedMultiVector destY =
     rcp_dynamic_cast<ProductMultiVectorBase<double> >(rcpFromPtr(y));

   // call apply
   implicitApply(srcX,destY,alpha,beta);
}

  void Export<LocalOrdinal,GlobalOrdinal,Node>::
  print (std::ostream& os) const
  {
    using Teuchos::Comm;
    using Teuchos::getFancyOStream;
    using Teuchos::RCP;
    using Teuchos::rcpFromRef;
    using Teuchos::toString;
    using std::endl;

    RCP<const Comm<int> > comm = getSourceMap ()->getComm ();
    const int myImageID = comm->getRank ();
    const int numImages = comm->getSize ();
    for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
      if (myImageID == imageCtr) {
        os << endl;
        if (myImageID == 0) { // I'm the root node (only output this info once)
          os << "Export Data Members:" << endl;
        }
        os << "Image ID       : " << myImageID << endl;

        os << "permuteFromLIDs: " << toString (getPermuteFromLIDs ()) << endl;
        os << "permuteToLIDs  : " << toString (getPermuteToLIDs ()) << endl;
        os << "remoteLIDs     : " << toString (getRemoteLIDs ()) << endl;
        os << "exportLIDs     : " << toString (getExportLIDs ()) << endl;
        os << "exportPIDs     : " << toString (getExportPIDs ()) << endl;

        os << "numSameIDs     : " << getNumSameIDs () << endl;
        os << "numPermuteIDs  : " << getNumPermuteIDs () << endl;
        os << "numRemoteIDs   : " << getNumRemoteIDs () << endl;
        os << "numExportIDs   : " << getNumExportIDs () << endl;
      }
      // A few global barriers give output a chance to complete.
      comm->barrier();
      comm->barrier();
      comm->barrier();
    }
    if (myImageID == 0) {
      os << endl << endl << "Source Map:" << endl << std::flush;
    }
    comm->barrier();
    os << *getSourceMap();
    comm->barrier();

    if (myImageID == 0) {
      os << endl << endl << "Target Map:" << endl << std::flush;
    }
    comm->barrier();
    os << *getTargetMap();
    comm->barrier();

    // It's also helpful for debugging to print the Distributor
    // object.  Epetra_Export::Print() does this, so we can do a
    // side-by-side comparison.
    if (myImageID == 0) {
      os << endl << endl << "Distributor:" << endl << std::flush;
    }
    comm->barrier();
    getDistributor().describe (*(getFancyOStream (rcpFromRef (os))),
                               Teuchos::VERB_EXTREME);
  }

  int
main (int argc, char *argv[])
{
  using Teuchos::inOutArg;
  using Teuchos::ParameterList;
  using Teuchos::parameterList;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::rcpFromRef;
  using std::endl;
  typedef double                          ST;
  typedef Epetra_Operator                 OP;
  typedef Epetra_MultiVector              MV;
  typedef Belos::OperatorTraits<ST,MV,OP> OPT;
  typedef Belos::MultiVecTraits<ST,MV>    MVT;

  // This calls MPI_Init and MPI_Finalize as necessary.
  Belos::Test::MPISession session (inOutArg (argc), inOutArg (argv));
  RCP<const Epetra_Comm> comm = session.getComm ();

  bool success = false;
  bool verbose = false;
  try {
    int MyPID = comm->MyPID ();

    //
    // Parameters to read from command-line processor
    //
    int frequency = -1;  // how often residuals are printed by solver
    int numRHS = 1;  // total number of right-hand sides to solve for
    int maxIters = 13000;  // maximum number of iterations for solver to use
    std::string filename ("bcsstk14.hb");
    double tol = 1.0e-5; // relative residual tolerance

    //
    // Read in command-line arguments
    //
    Teuchos::CommandLineProcessor cmdp (false, true);
    cmdp.setOption ("verbose", "quiet", &verbose, "Print messages and results.");
    cmdp.setOption ("frequency", &frequency, "Solvers frequency for printing "
        "residuals (#iters).");
    cmdp.setOption ("tol", &tol, "Relative residual tolerance used by MINRES "
        "solver.");
    cmdp.setOption ("filename", &filename, "Filename for Harwell-Boeing test "
        "matrix.");
    cmdp.setOption ("num-rhs", &numRHS, "Number of right-hand sides to solve.");
    cmdp.setOption ("max-iters", &maxIters, "Maximum number of iterations per "
        "linear system (-1 means \"adapt to problem/block size\").");
    if (cmdp.parse (argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
      return EXIT_FAILURE;
    }
    Teuchos::oblackholestream blackHole;
    std::ostream& verbOut = (verbose && MyPID == 0) ? std::cout : blackHole;

    //
    // Generate the linear system(s) to solve.
    //
    verbOut << "Generating the linear system(s) to solve" << endl << endl;
    RCP<Epetra_CrsMatrix> A;
    RCP<Epetra_MultiVector> B, X;
    RCP<Epetra_Map> rowMap;
    try {
      // This might change the number of right-hand sides, if we read in
      // a right-hand side from the Harwell-Boeing file.
      Belos::Util::createEpetraProblem (filename, &rowMap, &A, &B, &X, &MyPID, numRHS);
    } catch (std::exception& e) {
      TEUCHOS_TEST_FOR_EXCEPTION (true, std::runtime_error,
          "Failed to create Epetra problem for matrix "
          "filename \"" << filename << "\".  "
          "createEpetraProblem() reports the following "
          "error: " << e.what());
    }
    //
    // Compute the initial residual norm of the problem, so we can see
    // by how much it improved after the solve.
    //
    std::vector<double> initialResidualNorms (numRHS);
    std::vector<double> initialResidualInfNorms (numRHS);
    Epetra_MultiVector R (*rowMap, numRHS);
    OPT::Apply (*A, *X, R);
    MVT::MvAddMv (-1.0, R, 1.0, *B, R); // R := -(A*X) + B.
    MVT::MvNorm (R, initialResidualNorms);
    MVT::MvNorm (R, initialResidualInfNorms, Belos::InfNorm);
    if (verbose) {
      verbOut << "Initial residual 2-norms:            \t";
      for (int i = 0; i < numRHS; ++i) {
        verbOut << initialResidualNorms[i];
        if (i < numRHS-1) {
          verbOut << ", ";
        }
      }
      verbOut << endl << "Initial residual Inf-norms:          \t";
      for (int i = 0; i < numRHS; ++i) {
        verbOut << initialResidualInfNorms[i];
        if (i < numRHS-1) {
          verbOut << ", ";
        }
      }
      verbOut << endl;
    }
//.........这里部分代码省略.........

int main(int argc, char *argv[]) {
#include "MueLu_UseShortNames.hpp"

  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::rcpFromRef;
  using namespace MueLuTests;

#ifdef __GNUC__
#warning Navier2DBlocked_test based tests are disabled on 12/11/2013 due to some thrown exception
#endif
  return EXIT_SUCCESS;


  Teuchos::oblackholestream blackhole;
  Teuchos::GlobalMPISession mpiSession(&argc,&argv,&blackhole);
  //
  RCP<const Teuchos::Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();
  RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
  out->setOutputToRootOnly(0);
  *out << MueLu::MemUtils::PrintMemoryUsage() << std::endl;

  // Timing
  Teuchos::Time myTime("global");
  Teuchos::TimeMonitor MM(myTime);

  // read in some command line parameters
  Teuchos::CommandLineProcessor clp(false);

  int rebalanceBlock0 = 1;      clp.setOption("rebalanceBlock0",       &rebalanceBlock0,     "rebalance block 0 (1=yes, else=no)");
  int rebalanceBlock1 = 1;      clp.setOption("rebalanceBlock1",       &rebalanceBlock1,     "rebalance block 1 (1=yes, else=no)");

  switch (clp.parse(argc,argv)) {
    case Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED:        return EXIT_SUCCESS; break;
    case Teuchos::CommandLineProcessor::PARSE_ERROR:
    case Teuchos::CommandLineProcessor::PARSE_UNRECOGNIZED_OPTION: return EXIT_FAILURE; break;
    case Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL:                               break;
  }


#if defined(HAVE_MPI) && defined(HAVE_MUELU_ZOLTAN) && defined(HAVE_MUELU_ISORROPIA)
#ifndef HAVE_TEUCHOS_LONG_LONG_INT
  *out << "Warning: scaling test was not compiled with long long int support" << std::endl;


  // custom parameters
  LocalOrdinal maxLevels = 3;

  GlobalOrdinal maxCoarseSize=1; //FIXME clp doesn't like long long int

  int globalNumDofs = 8898;  // used for the maps
  int nDofsPerNode = 3;      // used for generating the fine level null-space

  // build strided maps
  // striding information: 2 velocity dofs and 1 pressure dof = 3 dofs per node
  std::vector<size_t> stridingInfo;
  stridingInfo.push_back(2);
  stridingInfo.push_back(1);

  /////////////////////////////////////// build strided maps
  // build strided maps:
  // xstridedfullmap: full map (velocity and pressure dof gids), continous
  // xstridedvelmap: only velocity dof gid maps (i.e. 0,1,3,4,6,7...)
  // xstridedpremap: only pressure dof gid maps (i.e. 2,5,8,...)
  Xpetra::UnderlyingLib lib = Xpetra::UseEpetra;
  RCP<StridedMap> xstridedfullmap = StridedMapFactory::Build(lib,globalNumDofs,0,stridingInfo,comm,-1);
  RCP<StridedMap> xstridedvelmap  = StridedMapFactory::Build(xstridedfullmap,0);
  RCP<StridedMap> xstridedpremap  = StridedMapFactory::Build(xstridedfullmap,1);

  /////////////////////////////////////// transform Xpetra::Map objects to Epetra
  // this is needed for AztecOO
  const RCP<const Epetra_Map> fullmap = rcpFromRef(Xpetra::toEpetra(*xstridedfullmap));
  RCP<const Epetra_Map>       velmap  = rcpFromRef(Xpetra::toEpetra(*xstridedvelmap));
  RCP<const Epetra_Map>       premap  = rcpFromRef(Xpetra::toEpetra(*xstridedpremap));

  /////////////////////////////////////// import problem matrix and RHS from files (-> Epetra)

  // read in problem
  Epetra_CrsMatrix * ptrA = 0;
  Epetra_Vector * ptrf = 0;
  Epetra_MultiVector* ptrNS = 0;

  *out << "Reading matrix market file" << std::endl;

  EpetraExt::MatrixMarketFileToCrsMatrix("A5932_re1000.txt",*fullmap,*fullmap,*fullmap,ptrA);
  EpetraExt::MatrixMarketFileToVector("b5932_re1000.txt",*fullmap,ptrf);
  //EpetraExt::MatrixMarketFileToCrsMatrix("/home/tobias/promotion/trilinos/fc17-dyn/packages/muelu/test/navierstokes/A5932_re1000.txt",*fullmap,*fullmap,*fullmap,ptrA);
  //EpetraExt::MatrixMarketFileToVector("/home/tobias/promotion/trilinos/fc17-dyn/packages/muelu/test/navierstokes/b5932_re1000.txt",*fullmap,ptrf);

  RCP<Epetra_CrsMatrix> epA = Teuchos::rcp(ptrA);
  RCP<Epetra_Vector> epv = Teuchos::rcp(ptrf);
  RCP<Epetra_MultiVector> epNS = Teuchos::rcp(ptrNS);


  /////////////////////////////////////// split system into 2x2 block system

  *out << "Split matrix into 2x2 block matrix" << std::endl;

  // split fullA into A11,..., A22
  Teuchos::RCP<Epetra_CrsMatrix> A11;
//.........这里部分代码省略.........

int 
main (int argc, char *argv[]) 
{
  using Teuchos::Comm;
  using Teuchos::FancyOStream;
  using Teuchos::getFancyOStream;
  using Teuchos::oblackholestream;
  using Teuchos::OSTab;
  using Teuchos::ParameterList;
  using Teuchos::parameterList;
  using Teuchos::RCP;
  using Teuchos::rcpFromRef;
  using std::cout;
  using std::endl;
  //
  // Typedefs for Tpetra template arguments.
  //
  typedef double scalar_type;
  typedef long int global_ordinal_type;
  typedef int local_ordinal_type;
  typedef Kokkos::DefaultNode::DefaultNodeType node_type;
  //
  // Tpetra objects which are the MV and OP template parameters of the
  // Belos specialization which we are testing.
  //
  typedef Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> MV;
  typedef Tpetra::Operator<scalar_type, local_ordinal_type, global_ordinal_type, node_type> OP;
  // 
  // Other typedefs.
  // 
  typedef Teuchos::ScalarTraits<scalar_type> STS;
  typedef Tpetra::CrsMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> sparse_matrix_type;

  Teuchos::GlobalMPISession mpiSession (&argc, &argv, &cout);
  RCP<const Comm<int> > comm = Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
  RCP<node_type> node = Tpetra::DefaultPlatform::getDefaultPlatform().getNode();
  RCP<oblackholestream> blackHole (new oblackholestream);
  const int myRank = comm->getRank();

  // Output stream that prints only on Rank 0.
  RCP<FancyOStream> out;
  if (myRank == 0) {
    out = Teuchos::getFancyOStream (rcpFromRef (cout));
  } else {
    out = Teuchos::getFancyOStream (blackHole);
  }

  //
  // Get test parameters from command-line processor.
  //  
  // CommandLineProcessor always understands int, but may not
  // understand global_ordinal_type.  We convert to the latter below.
  int numRows = comm->getSize() * 100;
  bool tolerant = false;
  bool verbose = false;
  bool debug = false;
  Teuchos::CommandLineProcessor cmdp (false, true);
  cmdp.setOption("numRows", &numRows,
		 "Global number of rows (and columns) in the sparse matrix to generate.");
  cmdp.setOption("tolerant", "intolerant", &tolerant,
		 "Whether to parse files tolerantly.");
  cmdp.setOption("verbose", "quiet", &verbose, 
		 "Print messages and results.");
  cmdp.setOption("debug", "release", &debug, 
		 "Run debugging checks and print copious debugging output.");
  if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
    *out << "\nEnd Result: TEST FAILED" << endl;
    return EXIT_FAILURE;
  }
  // Output stream for verbose output.
  RCP<FancyOStream> verbOut = verbose ? out : getFancyOStream (blackHole);

  const bool success = true;

  // Test whether it's possible to instantiate the solver.
  // This is a minimal compilation test.
  *verbOut << "Instantiating Block GCRODR solver" << endl;
  Belos::BlockGCRODRSolMgr<scalar_type, MV, OP> solver;
  //
  // Test setting solver parameters.  For now, we just use an empty
  // (but non-null) parameter list, which the solver should fill in
  // with defaults.
  //
  *verbOut << "Setting solver parameters" << endl;
  RCP<ParameterList> solverParams = parameterList ();
  solver.setParameters (solverParams);
  //
  // Create a linear system to solve.
  //
  *verbOut << "Creating linear system" << endl;
  RCP<sparse_matrix_type> A;
  RCP<MV> X_guess, X_exact, B;
  {
    typedef Belos::Tpetra::ProblemMaker<sparse_matrix_type> factory_type;
    factory_type factory (comm, node, out, tolerant, debug);
    
    RCP<ParameterList> problemParams = parameterList ();
    problemParams->set ("Global number of rows", 
			static_cast<global_ordinal_type> (numRows));
    problemParams->set ("Problem type", std::string ("Nonsymmetric"));
//.........这里部分代码省略.........

Teuchos::RCP<const Teuchos::ParameterList>
LSQRSolMgr<ScalarType,MV,OP,false>::getValidParameters() const
{
    using Teuchos::ParameterList;
    using Teuchos::parameterList;
    using Teuchos::RCP;
    using Teuchos::rcp;
    using Teuchos::rcpFromRef;
    typedef Teuchos::ScalarTraits<MagnitudeType> STM;

    // Set all the valid parameters and their default values.
    if (is_null (validParams_)) {
        // We use Teuchos::as just in case MagnitudeType doesn't have a
        // constructor that takes an int.  Otherwise, we could just write
        // "MagnitudeType(10)".
        const MagnitudeType ten = Teuchos::as<MagnitudeType> (10);
        const MagnitudeType sqrtEps = STM::squareroot (STM::eps());

        const MagnitudeType lambda = STM::zero();
        RCP<std::ostream> outputStream = rcpFromRef (std::cout);
        const MagnitudeType relRhsErr = ten * sqrtEps;
        const MagnitudeType relMatErr = ten * sqrtEps;
        const MagnitudeType condMax = STM::one() / STM::eps();
        const int maxIters = 1000;
        const int termIterMax = 1;
        const std::string orthoType ("DGKS");
        const MagnitudeType orthoKappa = Teuchos::as<MagnitudeType> (-1.0);
        const int verbosity = Belos::Errors;
        const int outputStyle = Belos::General;
        const int outputFreq = -1;
        const std::string label ("Belos");

        RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
        pl->set("Output Stream", outputStream,
                "is a reference-counted pointer to the output stream receiving\n"
                "all solver output.");
        pl->set("Lambda", lambda, "is the damping parameter.");
        pl->set("Rel RHS Err", relRhsErr,
                "estimates the error in the data defining the right-\n"
                "hand side.");
        pl->set("Rel Mat Err", relMatErr,
                "estimates the error in the data defining the matrix.");
        pl->set("Condition Limit", condMax,
                "bounds the estimated condition number of Abar.");
        pl->set("Maximum Iterations", maxIters,
                "allows at most the maximum number of iterations.");
        pl->set("Term Iter Max", termIterMax,
                "consecutive iterations meeting thresholds are necessary for\n"
                "for convergence.");
        pl->set("Orthogonalization", orthoType,
                "uses orthogonalization of either DGKS, ICGS, IMGS, or TSQR.");
        {
            OrthoManagerFactory<ScalarType, MV, OP> factory;
            pl->set("Orthogonalization", orthoType,
                    "refers to the orthogonalization method to use.  Valid "
                    "options: " + factory.validNamesString());
            RCP<const ParameterList> orthoParams =
                factory.getDefaultParameters (orthoType);
            pl->set ("Orthogonalization Parameters", *orthoParams,
                     "Parameters specific to the type of orthogonalization used.");
        }
        pl->set("Orthogonalization Constant", orthoKappa,
                "is the threshold used by DGKS orthogonalization to determine\n"
                "whether or not to repeat classical Gram-Schmidt.  This parameter\n"
                "is ignored if \"Orthogonalization\" is not \"DGKS\".");
        pl->set("Verbosity", verbosity,
                "type(s) of solver information are outputted to the output\n"
                "stream.");
        pl->set("Output Style", outputStyle,
                "the style used for the solver information outputted to the\n"
                "output stream.");
        pl->set("Output Frequency", outputFreq,
                "is the frequency at which information is written to the\n"
                "output stream.");
        pl->set("Timer Label", label,
                "is the string to use as a prefix for the timer labels.");
        //  pl->set("Restart Timers", restartTimers_);
        pl->set("Block Size", 1, "Block size parameter (currently, this must always be 1).");

        validParams_ = pl;
    }
    return validParams_;
}

void DenseContainer<MatrixType, LocalScalarType>::
applyImpl (const local_mv_type& X,
           local_mv_type& Y,
           Teuchos::ETransp mode,
           LocalScalarType alpha,
           LocalScalarType beta) const
{
    using Teuchos::ArrayRCP;
    using Teuchos::RCP;
    using Teuchos::rcp;
    using Teuchos::rcpFromRef;

    TEUCHOS_TEST_FOR_EXCEPTION(
        X.getLocalLength () != Y.getLocalLength (),
        std::logic_error, "Ifpack2::DenseContainer::applyImpl: X and Y have "
        "incompatible dimensions (" << X.getLocalLength () << " resp. "
        << Y.getLocalLength () << ").  Please report this bug to "
        "the Ifpack2 developers.");
    TEUCHOS_TEST_FOR_EXCEPTION(
        localMap_->getNodeNumElements () != X.getLocalLength (),
        std::logic_error, "Ifpack2::DenseContainer::applyImpl: The inverse "
        "operator and X have incompatible dimensions (" <<
        localMap_->getNodeNumElements () << " resp. "
        << X.getLocalLength () << ").  Please report this bug to "
        "the Ifpack2 developers.");
    TEUCHOS_TEST_FOR_EXCEPTION(
        localMap_->getNodeNumElements () != Y.getLocalLength (),
        std::logic_error, "Ifpack2::DenseContainer::applyImpl: The inverse "
        "operator and Y have incompatible dimensions (" <<
        localMap_->getNodeNumElements () << " resp. "
        << Y.getLocalLength () << ").  Please report this bug to "
        "the Ifpack2 developers.");
    TEUCHOS_TEST_FOR_EXCEPTION(
        X.getLocalLength () != static_cast<size_t> (diagBlock_.numRows ()),
        std::logic_error, "Ifpack2::DenseContainer::applyImpl: The input "
        "multivector X has incompatible dimensions from those of the "
        "inverse operator (" << X.getLocalLength () << " vs. "
        << (mode == Teuchos::NO_TRANS ? diagBlock_.numCols () : diagBlock_.numRows ())
        << ").  Please report this bug to the Ifpack2 developers.");
    TEUCHOS_TEST_FOR_EXCEPTION(
        X.getLocalLength () != static_cast<size_t> (diagBlock_.numRows ()),
        std::logic_error, "Ifpack2::DenseContainer::applyImpl: The output "
        "multivector Y has incompatible dimensions from those of the "
        "inverse operator (" << Y.getLocalLength () << " vs. "
        << (mode == Teuchos::NO_TRANS ? diagBlock_.numRows () : diagBlock_.numCols ())
        << ").  Please report this bug to the Ifpack2 developers.");

    typedef Teuchos::ScalarTraits<local_scalar_type> STS;
    const int numVecs = static_cast<int> (X.getNumVectors ());
    if (alpha == STS::zero ()) { // don't need to solve the linear system
        if (beta == STS::zero ()) {
            // Use BLAS AXPY semantics for beta == 0: overwrite, clobbering
            // any Inf or NaN values in Y (rather than multiplying them by
            // zero, resulting in NaN values).
            Y.putScalar (STS::zero ());
        }
        else { // beta != 0
            Y.scale (beta);
        }
    }
    else { // alpha != 0; must solve the linear system
        Teuchos::LAPACK<int, local_scalar_type> lapack;
        // If beta is nonzero or Y is not constant stride, we have to use
        // a temporary output multivector.  It gets a (deep) copy of X,
        // since GETRS overwrites its (multi)vector input with its output.
        RCP<local_mv_type> Y_tmp;
        if (beta == STS::zero () ) {
            Tpetra::deep_copy (Y, X);
            Y_tmp = rcpFromRef (Y);
        }
        else {
            Y_tmp = rcp (new local_mv_type (X, Teuchos::Copy));
        }
        const int Y_stride = static_cast<int> (Y_tmp->getStride ());
        ArrayRCP<local_scalar_type> Y_view = Y_tmp->get1dViewNonConst ();
        local_scalar_type* const Y_ptr = Y_view.getRawPtr ();

        int INFO = 0;
        const char trans =
            (mode == Teuchos::CONJ_TRANS ? 'C' : (mode == Teuchos::TRANS ? 'T' : 'N'));
        lapack.GETRS (trans, diagBlock_.numRows (), numVecs,
                      diagBlock_.values (), diagBlock_.stride (),
                      ipiv_.getRawPtr (), Y_ptr, Y_stride, &INFO);
        TEUCHOS_TEST_FOR_EXCEPTION(
            INFO != 0, std::runtime_error, "Ifpack2::DenseContainer::applyImpl: "
            "LAPACK's _GETRS (solve using LU factorization with partial pivoting) "
            "failed with INFO = " << INFO << " != 0.");

        if (beta != STS::zero ()) {
            Y.update (alpha, *Y_tmp, beta);
        }
        else if (! Y.isConstantStride ()) {
            Tpetra::deep_copy (Y, *Y_tmp);
        }
    }
}

int main(int argc, char *argv[]){
  Teuchos::GlobalMPISession mpiSession(&argc,&argv);
  typedef float SCALAR;
  typedef int LO;
  typedef int GO;
  typedef std::complex<SCALAR> cmplx;
  typedef CrsMatrix<cmplx,LO,GO,Node> MAT;
  typedef ScalarTraits<cmplx> ST;
  typedef MultiVector<cmplx,LO,GO,Node> MV;
  typedef ST::magnitudeType Mag;
  typedef ScalarTraits<Mag> MT;
  const size_t numVecs = 1;

  std::ostream &out = std::cout;
  RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(out));

  Platform &platform = Tpetra::DefaultPlatform::getDefaultPlatform();
  RCP<const Comm<int> > comm = platform.getComm();
  RCP<Node>             node = platform.getNode();

  RCP<MAT> A =
    Tpetra::MatrixMarket::Reader<MAT>::readSparseFile("../test/matrices/amesos2_test_mat3.mtx",comm,node);

  RCP<const Map<LO,GO,Node> > dmnmap = A->getDomainMap();
  RCP<const Map<LO,GO,Node> > rngmap = A->getRangeMap();

  // Create the know-solution vector
  std::map<GO,cmplx> xValues;
  xValues[0] = cmplx(-0.58657, 0.10646);
  xValues[1] = cmplx(0.86716, 0.75421); 
  xValues[2] = cmplx(0.58970, 0.29876);

  std::map<GO,cmplx>::iterator it;
  RCP<MV> X = rcp(new MV(dmnmap, numVecs));
  X->setObjectLabel("X");

  for( it = xValues.begin(); it != xValues.end(); ++it ){
    if( rngmap->isNodeGlobalElement( (*it).first ) ){
      X->replaceGlobalValue( (*it).first, 0, (*it).second );
    }
  }

  RCP<MV> Xhat = rcp(new MV(dmnmap, numVecs));
  RCP<MV> B    = rcp(new MV(rngmap, numVecs));
  X->setObjectLabel("X");
  B->setObjectLabel("B");
  Xhat->setObjectLabel("Xhat");

  A->apply(*X,*B,Teuchos::CONJ_TRANS); // use conjugate-transpose

  Xhat->randomize();

  // Solve A*Xhat = B for Xhat using the Superlu solver
  RCP<Amesos2::Solver<MAT,MV> > solver
    = Amesos2::create<MAT,MV>("SuperLU", A, Xhat, B );

  Teuchos::ParameterList amesos2_params;
  // Setting the following will cause Amesos2 to complain
  amesos2_params.sublist("SuperLU").set("Trans","CONJ","Solve with conjugate-transpose");

  solver->setParameters( rcpFromRef(amesos2_params) );
  solver->symbolicFactorization().numericFactorization().solve();

  B->describe(*fos, Teuchos::VERB_EXTREME);
  Xhat->describe(*fos, Teuchos::VERB_EXTREME);
  X->describe(*fos, Teuchos::VERB_EXTREME);

  // Check result of solve
  Array<Mag> xhatnorms(numVecs), xnorms(numVecs);
  Xhat->norm2(xhatnorms());
  X->norm2(xnorms());
}