C++ SerialDenseMatrix::putScalar方法代码示例

//.........这里部分代码省略.........
  M->SetTracebackMode (1); // Shutdown Epetra Warning tracebacks

  //************************************
  // Call the LOBPCG solver manager
  //***********************************
  //
  // Variables used for the LOBPCG Method
  const int nev       = 10;
  const int blockSize = 5;
  const int maxIters  = 500;
  const double tol    = 1.0e-8;

  typedef Epetra_MultiVector MV;
  typedef Epetra_Operator OP;
  typedef MultiVecTraits<double, Epetra_MultiVector> MVT;

  // Create an Epetra_MultiVector for an initial vector to start the
  // solver.  Note: This needs to have the same number of columns as
  // the blocksize.
  RCP<Epetra_MultiVector> ivec = rcp (new Epetra_MultiVector (Map, blockSize));
  ivec->Random (); // fill the initial vector with random values

  // Create the eigenproblem.
  RCP<BasicEigenproblem<double, MV, OP> > MyProblem =
    rcp (new BasicEigenproblem<double, MV, OP> (A, ivec));

  // Inform the eigenproblem that the operator A is symmetric
  MyProblem->setHermitian (true);

  // Set the number of eigenvalues requested
  MyProblem->setNEV (nev);

  // Tell the eigenproblem that you are finishing passing it information.
  const bool success = MyProblem->setProblem ();
  if (! success) {
    printer.print (Errors, "Anasazi::BasicEigenproblem::setProblem() reported an error.\n");
#ifdef HAVE_MPI
    MPI_Finalize ();
#endif // HAVE_MPI
    return -1;
  }

  // Create parameter list to pass into the solver manager
  Teuchos::ParameterList MyPL;
  MyPL.set ("Which", which);
  MyPL.set ("Block Size", blockSize);
  MyPL.set ("Maximum Iterations", maxIters);
  MyPL.set ("Convergence Tolerance", tol);
  MyPL.set ("Full Ortho", true);
  MyPL.set ("Use Locking", true);

  // Create the solver manager
  LOBPCGSolMgr<double, MV, OP> MySolverMan (MyProblem, MyPL);

  // Solve the problem
  ReturnType returnCode = MySolverMan.solve ();

  // Get the eigenvalues and eigenvectors from the eigenproblem
  Eigensolution<double,MV> sol = MyProblem->getSolution ();
  std::vector<Value<double> > evals = sol.Evals;
  RCP<MV> evecs = sol.Evecs;

  // Compute residuals.
  std::vector<double> normR (sol.numVecs);
  if (sol.numVecs > 0) {
    Teuchos::SerialDenseMatrix<int,double> T (sol.numVecs, sol.numVecs);
    Epetra_MultiVector tempAevec (Map, sol.numVecs );
    T.putScalar (0.0);
    for (int i = 0; i < sol.numVecs; ++i) {
      T(i,i) = evals[i].realpart;
    }
    A->Apply (*evecs, tempAevec);
    MVT::MvTimesMatAddMv (-1.0, *evecs, T, 1.0, tempAevec);
    MVT::MvNorm (tempAevec, normR);
  }

  // Print the results
  std::ostringstream os;
  os.setf (std::ios_base::right, std::ios_base::adjustfield);
  os << "Solver manager returned "
     << (returnCode == Converged ? "converged." : "unconverged.") << endl;
  os << endl;
  os << "------------------------------------------------------" << endl;
  os << std::setw(16) << "Eigenvalue"
     << std::setw(18) << "Direct Residual"
     << endl;
  os << "------------------------------------------------------" << endl;
  for (int i = 0; i < sol.numVecs; ++i) {
    os << std::setw(16) << evals[i].realpart
       << std::setw(18) << normR[i] / evals[i].realpart
       << endl;
  }
  os << "------------------------------------------------------" << endl;
  printer.print (Errors, os.str ());

#ifdef HAVE_MPI
  MPI_Finalize ();
#endif // HAVE_MPI
  return 0;
}

    void
    factorExplicit (Kokkos::MultiVector<Scalar, NodeType>& A,
		    Kokkos::MultiVector<Scalar, NodeType>& Q,
		    Teuchos::SerialDenseMatrix<LocalOrdinal, Scalar>& R,
		    const bool contiguousCacheBlocks,
		    const bool forceNonnegativeDiagonal=false)
    {
      using Teuchos::asSafe;
      typedef Kokkos::MultiVector<Scalar, NodeType> KMV;

      // Tsqr currently likes LocalOrdinal ordinals, but
      // Kokkos::MultiVector has size_t ordinals.  Do conversions
      // here.  
      //
      // Teuchos::asSafe() can do safe conversion (e.g., checking for
      // overflow when casting to a narrower integer type), if a
      // custom specialization is defined for
      // Teuchos::ValueTypeConversionTraits<size_t, LocalOrdinal>.
      // Otherwise, this has the same (potentially) unsafe effect as
      // static_cast<LocalOrdinal>(...) would have.
      const LocalOrdinal A_numRows = asSafe<LocalOrdinal> (A.getNumRows());
      const LocalOrdinal A_numCols = asSafe<LocalOrdinal> (A.getNumCols());
      const LocalOrdinal A_stride = asSafe<LocalOrdinal> (A.getStride());
      const LocalOrdinal Q_numRows = asSafe<LocalOrdinal> (Q.getNumRows());
      const LocalOrdinal Q_numCols = asSafe<LocalOrdinal> (Q.getNumCols());
      const LocalOrdinal Q_stride = asSafe<LocalOrdinal> (Q.getStride());

      // Sanity checks for matrix dimensions
      if (A_numRows < A_numCols) {
	std::ostringstream os;
	os << "In Tsqr::factorExplicit: input matrix A has " << A_numRows 
	   << " local rows, and " << A_numCols << " columns.  The input "
	  "matrix must have at least as many rows on each processor as "
	  "there are columns.";
	throw std::invalid_argument(os.str());
      } else if (A_numRows != Q_numRows) {
	std::ostringstream os;
	os << "In Tsqr::factorExplicit: input matrix A and output matrix Q "
	  "must have the same number of rows.  A has " << A_numRows << " rows"
	  " and Q has " << Q_numRows << " rows.";
	throw std::invalid_argument(os.str());
      } else if (R.numRows() < R.numCols()) {
	std::ostringstream os;
	os << "In Tsqr::factorExplicit: output matrix R must have at least "
	  "as many rows as columns.  R has " << R.numRows() << " rows and "
	   << R.numCols() << " columns.";
	throw std::invalid_argument(os.str());
      } else if (A_numCols != R.numCols()) {
	std::ostringstream os;
	os << "In Tsqr::factorExplicit: input matrix A and output matrix R "
	  "must have the same number of columns.  A has " << A_numCols 
	   << " columns and R has " << R.numCols() << " columns.";
	throw std::invalid_argument(os.str());
      }

      // Check for quick exit, based on matrix dimensions
      if (Q_numCols == 0)
	return;

      // Hold on to nonconst views of A and Q.  This will make TSQR
      // correct (if perhaps inefficient) for all possible Kokkos Node
      // types, even GPU nodes.
      Teuchos::ArrayRCP<scalar_type> A_ptr = A.getValuesNonConst();
      Teuchos::ArrayRCP<scalar_type> Q_ptr = Q.getValuesNonConst();

      R.putScalar (STS::zero());
      NodeOutput nodeResults = 
	nodeTsqr_->factor (A_numRows, A_numCols, A_ptr.getRawPtr(), A_stride,
			   R.values(), R.stride(), contiguousCacheBlocks);
      // FIXME (mfh 19 Oct 2010) Replace actions on raw pointer with
      // actions on the Kokkos::MultiVector or at least the ArrayRCP.
      nodeTsqr_->fill_with_zeros (Q_numRows, Q_numCols, 
				  Q_ptr.getRawPtr(), Q_stride,
				  contiguousCacheBlocks);
      matview_type Q_rawView (Q_numRows, Q_numCols, 
			      Q_ptr.getRawPtr(), Q_stride);
      matview_type Q_top_block = 
	nodeTsqr_->top_block (Q_rawView, contiguousCacheBlocks);
      if (Q_top_block.nrows() < R.numCols()) {
	std::ostringstream os;
	os << "The top block of Q has too few rows.  This means that the "
	   << "the intranode TSQR implementation has a bug in its top_block"
	   << "() method.  The top block should have at least " << R.numCols()
	   << " rows, but instead has only " << Q_top_block.ncols() 
	   << " rows.";
	throw std::logic_error (os.str());
      }
      {
	matview_type Q_top (R.numCols(), Q_numCols, Q_top_block.get(), 
			    Q_top_block.lda());
	matview_type R_view (R.numRows(), R.numCols(), R.values(), R.stride());
	distTsqr_->factorExplicit (R_view, Q_top, forceNonnegativeDiagonal);
      }
      nodeTsqr_->apply (ApplyType::NoTranspose, 
			A_numRows, A_numCols, A_ptr.getRawPtr(), A_stride,
			nodeResults, Q_numCols, Q_ptr.getRawPtr(), Q_stride,
			contiguousCacheBlocks);

      // If necessary, force the R factor to have a nonnegative diagonal.
      if (forceNonnegativeDiagonal && 
//.........这里部分代码省略.........