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

      /// \brief Verify the result of the "thin" QR factorization \f$A = QR\f$.
      ///
      /// This method returns a list of three magnitudes: 
      /// - \f$\| A - QR \|_F\f$
      /// - \f$\|I - Q^* Q\|_F\f$
      /// - \f$\|A\|_F\f$
      ///
      /// The notation $\f\| X \|\f$ denotes the Frobenius norm
      /// (square root of sum of squares) of a matrix \f$X\f$.
      /// Returning the Frobenius norm of \f$A\f$ allows you to scale
      /// or not scale the residual \f$\|A - QR\|\f$ as you prefer.
      virtual std::vector< magnitude_type >
      verify (const multivector_type& A,
	      const multivector_type& Q,
	      const Teuchos::SerialDenseMatrix< local_ordinal_type, scalar_type >& R)
      {
	using Teuchos::ArrayRCP;

	local_ordinal_type nrowsLocal_A, ncols_A, LDA;
	local_ordinal_type nrowsLocal_Q, ncols_Q, LDQ;
	fetchDims (A, nrowsLocal_A, ncols_A, LDA);
	fetchDims (Q, nrowsLocal_Q, ncols_Q, LDQ);
	if (nrowsLocal_A != nrowsLocal_Q)
	  throw std::runtime_error ("A and Q must have same number of rows");
	else if (ncols_A != ncols_Q)
	  throw std::runtime_error ("A and Q must have same number of columns");
	else if (ncols_A != R.numCols())
	  throw std::runtime_error ("A and R must have same number of columns");
	else if (R.numRows() < R.numCols())
	  throw std::runtime_error ("R must have no fewer rows than columns");

	// Const views suffice for verification
	ArrayRCP< const scalar_type > A_ptr = fetchConstView (A);
	ArrayRCP< const scalar_type > Q_ptr = fetchConstView (Q);
	return global_verify (nrowsLocal_A, ncols_A, A_ptr.get(), LDA,
			      Q_ptr.get(), LDQ, R.values(), R.stride(), 
			      pScalarMessenger_.get());
      }

void
Stokhos::SmolyakPseudoSpectralOperator<ordinal_type,value_type,point_compare_type>::
apply_direct(
  const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& A,
  const value_type& alpha, 
  const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input,
  Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result, 
  const value_type& beta,
  bool trans) const {
  if (trans) {
    TEUCHOS_ASSERT(input.numCols() <= A.numCols());
    TEUCHOS_ASSERT(result.numCols() == A.numRows());
    TEUCHOS_ASSERT(result.numRows() == input.numRows());
    blas.GEMM(Teuchos::NO_TRANS, Teuchos::TRANS, input.numRows(), 
	      A.numRows(), input.numCols(), alpha, input.values(), 
	      input.stride(), A.values(), A.stride(), beta, 
	      result.values(), result.stride());
  }
  else {
    TEUCHOS_ASSERT(input.numRows() <= A.numCols());
    TEUCHOS_ASSERT(result.numRows() == A.numRows());
    TEUCHOS_ASSERT(result.numCols() == input.numCols());
    blas.GEMM(Teuchos::NO_TRANS, Teuchos::NO_TRANS, A.numRows(), 
	      input.numCols(), input.numRows(), alpha, A.values(), 
	      A.stride(), input.values(), input.stride(), beta, 
	      result.values(), result.stride());
  }
}

 void EpetraOpMultiVec::MvTimesMatAddMv ( double alpha, const MultiVec<double>& A, 
     const Teuchos::SerialDenseMatrix<int,double>& B, double beta ) 
 {
   Epetra_LocalMap LocalMap(B.numRows(), 0, Epetra_MV->Map().Comm());
   Epetra_MultiVector B_Pvec(Epetra_DataAccess::View, LocalMap, B.values(), B.stride(), B.numCols());
   
   EpetraOpMultiVec *A_vec = dynamic_cast<EpetraOpMultiVec *>(&const_cast<MultiVec<double> &>(A)); 
   TEUCHOS_TEST_FOR_EXCEPTION( A_vec==NULL,  std::invalid_argument, "Anasazi::EpetraOpMultiVec::SetBlocks() cast of MultiVec<double> to EpetraOpMultiVec failed.");
   
   TEUCHOS_TEST_FOR_EXCEPTION( 
       Epetra_MV->Multiply( 'N', 'N', alpha, *(A_vec->GetEpetraMultiVector()), B_Pvec, beta ) != 0,
       EpetraSpecializedMultiVecFailure, "Anasazi::EpetraOpMultiVec::MvTimesMatAddMv() call to Epetra_MultiVec::Multiply() returned a nonzero value.");
 }

void updateGuess(Teuchos::SerialDenseVector<int, std::complex<double> >& myCurrentGuess,
		Teuchos::SerialDenseVector<int, std::complex<double> >& myTargetsCalculated,
		Teuchos::SerialDenseMatrix<int, std::complex<double> >& myJacobian, 
		Teuchos::LAPACK<int, std::complex<double> >& myLAPACK
		 )
{
	//v = J(inverse) * (-F(x))
	//new guess = v + old guess
	myTargetsCalculated *= -1.0;

	//Perform an LU factorization of this matrix. 
	int ipiv[NUMDIMENSIONS], info;
	char TRANS = 'N';
	myLAPACK.GETRF( NUMDIMENSIONS, NUMDIMENSIONS, myJacobian.values(), myJacobian.stride(), ipiv, &info ); 

	// Solve the linear system.
	myLAPACK.GETRS( TRANS, NUMDIMENSIONS, 1, myJacobian.values(), myJacobian.stride(),
		       	ipiv, myTargetsCalculated.values(), myTargetsCalculated.stride(), &info );  

	//We have overwritten myTargetsCalculated with guess update values
	//myBLAS.AXPY(NUMDIMENSIONS, 1.0, myGuessAdjustment.values(), 1, myCurrentGuess.values(), 1);
	myCurrentGuess += myTargetsCalculated;
}

void EpetraMultiVec::MvTransMv ( const double alpha, const MultiVec<double>& A,
				 Teuchos::SerialDenseMatrix<int,double>& B) const
{    
  EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
  
  if (A_vec) {
    Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
    Epetra_MultiVector B_Pvec(View, LocalMap, B.values(), B.stride(), B.numCols());
    
    int info = B_Pvec.Multiply( 'T', 'N', alpha, *A_vec, *this, 0.0 );
    TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
		       "Belos::EpetraMultiVec::MvTransMv call to Multiply() returned a nonzero value.");
  }
}

void EpetraMultiVec::MvTimesMatAddMv ( const double alpha, const MultiVec<double>& A, 
				       const Teuchos::SerialDenseMatrix<int,double>& B, const double beta ) 
{
  Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
  Epetra_MultiVector B_Pvec(View, LocalMap, B.values(), B.stride(), B.numCols());
  
  EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A)); 
  TEST_FOR_EXCEPTION(A_vec==NULL, EpetraMultiVecFailure,
                     "Belos::EpetraMultiVec::MvTimesMatAddMv cast from Belos::MultiVec<> to Belos::EpetraMultiVec failed.");
  
  int info = Multiply( 'N', 'N', alpha, *A_vec, B_Pvec, beta );
  TEST_FOR_EXCEPTION(info!=0, EpetraMultiVecFailure, 
		     "Belos::EpetraMultiVec::MvTimesMatAddMv call to Multiply() returned a nonzero value.");

}

    LocalOrdinal
    revealRank (Kokkos::MultiVector<Scalar, NodeType>& Q,
		Teuchos::SerialDenseMatrix<LocalOrdinal, Scalar>& R,
		const magnitude_type& tol,
		const bool contiguousCacheBlocks = false) const
    {
      typedef Kokkos::MultiVector<Scalar, NodeType> KMV;

      const LocalOrdinal nrows = static_cast<LocalOrdinal> (Q.getNumRows());
      const LocalOrdinal ncols = static_cast<LocalOrdinal> (Q.getNumCols());
      const LocalOrdinal ldq = static_cast<LocalOrdinal> (Q.getStride());
      Teuchos::ArrayRCP<Scalar> Q_ptr = Q.getValuesNonConst();

      // Take the easy exit if available.
      if (ncols == 0)
	return 0;

      //
      // FIXME (mfh 16 Jul 2010) We _should_ compute the SVD of R (as
      // the copy B) on Proc 0 only.  This would ensure that all
      // processors get the same SVD and rank (esp. in a heterogeneous
      // computing environment).  For now, we just do this computation
      // redundantly, and hope that all the returned rank values are
      // the same.
      //
      matrix_type U (ncols, ncols, STS::zero());
      const ordinal_type rank = 
	reveal_R_rank (ncols, R.values(), R.stride(), 
		       U.get(), U.lda(), tol);
      if (rank < ncols)
	{
	  // cerr << ">>> Rank of R: " << rank << " < ncols=" << ncols << endl;
	  // cerr << ">>> Resulting U:" << endl;
	  // print_local_matrix (cerr, ncols, ncols, R, ldr);
	  // cerr << endl;

	  // If R is not full rank: reveal_R_rank() already computed
	  // the SVD \f$R = U \Sigma V^*\f$ of (the input) R, and
	  // overwrote R with \f$\Sigma V^*\f$.  Now, we compute \f$Q
	  // := Q \cdot U\f$, respecting cache blocks of Q.
	  Q_times_B (nrows, ncols, Q_ptr.getRawPtr(), ldq, 
		     U.get(), U.lda(), contiguousCacheBlocks);
	}
      return rank;
    }

  void EpetraMultiVec::MvTransMv ( double alpha, const MultiVec<double>& A,
                                   Teuchos::SerialDenseMatrix<int,double>& B
#ifdef HAVE_ANASAZI_EXPERIMENTAL
                                   , ConjType conj
#endif
                                  ) const
  {    
    EpetraMultiVec *A_vec = dynamic_cast<EpetraMultiVec *>(&const_cast<MultiVec<double> &>(A));
    
    if (A_vec) {
      Epetra_LocalMap LocalMap(B.numRows(), 0, Map().Comm());
      Epetra_MultiVector B_Pvec(View, LocalMap, B.values(), B.stride(), B.numCols());
      
    TEUCHOS_TEST_FOR_EXCEPTION( 
        B_Pvec.Multiply( 'T', 'N', alpha, *A_vec, *this, 0.0 ) != 0,
        EpetraMultiVecFailure, "Anasazi::EpetraMultiVec::MvTransMv() call to Epetra_MultiVec::Multiply() returned a nonzero value.");
    }
  }

  void EpetraOpMultiVec::MvTransMv ( double alpha, const MultiVec<double>& A,
                                   Teuchos::SerialDenseMatrix<int,double>& B
#ifdef HAVE_ANASAZI_EXPERIMENTAL
                                   , ConjType conj
#endif
                                  ) const
  {    
    EpetraOpMultiVec *A_vec = dynamic_cast<EpetraOpMultiVec *>(&const_cast<MultiVec<double> &>(A));
    
    if (A_vec) {
      Epetra_LocalMap LocalMap(B.numRows(), 0, Epetra_MV->Map().Comm());
      Epetra_MultiVector B_Pvec(Epetra_DataAccess::View, LocalMap, B.values(), B.stride(), B.numCols());
     
      int info = Epetra_OP->Apply( *Epetra_MV, *Epetra_MV_Temp );
      TEUCHOS_TEST_FOR_EXCEPTION( info != 0, EpetraSpecializedMultiVecFailure, 
        "Anasazi::EpetraOpMultiVec::MvTransMv(): Error returned from Epetra_Operator::Apply()" );

      TEUCHOS_TEST_FOR_EXCEPTION( 
        B_Pvec.Multiply( 'T', 'N', alpha, *(A_vec->GetEpetraMultiVector()), *Epetra_MV_Temp, 0.0 ) != 0,
        EpetraSpecializedMultiVecFailure, "Anasazi::EpetraOpMultiVec::MvTransMv() call to Epetra_MultiVector::Multiply() returned a nonzero value.");
    }
  }

    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 && 
//.........这里部分代码省略.........