C++ teuchos::SerialDenseMatrix类代码示例

      /// \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 assembleIRKState(
  const int stageIndex,
  const Teuchos::SerialDenseMatrix<int,Scalar> &A_in,
  const Scalar dt,
  const Thyra::VectorBase<Scalar> &x_base,
  const Thyra::ProductVectorBase<Scalar> &x_stage_bar,
  Teuchos::Ptr<Thyra::VectorBase<Scalar> > x_out_ptr
  )
{

  typedef ScalarTraits<Scalar> ST;

  const int numStages_in = A_in.numRows();
  TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( stageIndex, 0, numStages_in );
  TEUCHOS_ASSERT_EQUALITY( A_in.numRows(), numStages_in );
  TEUCHOS_ASSERT_EQUALITY( A_in.numCols(), numStages_in );
  TEUCHOS_ASSERT_EQUALITY( x_stage_bar.productSpace()->numBlocks(), numStages_in );
  Thyra::VectorBase<Scalar>& x_out = *x_out_ptr;

  V_V( outArg(x_out), x_base );
  for ( int j = 0; j < numStages_in; ++j ) {
    Vp_StV( outArg(x_out), dt * A_in(stageIndex,j), *x_stage_bar.getVectorBlock(j) );
  }

}

void
Stokhos::SmolyakPseudoSpectralOperator<ordinal_type,value_type,point_compare_type>::
transformPCE2QP_smolyak(
  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 {
  Teuchos::SerialDenseMatrix<ordinal_type,value_type> op_input, op_result;
  result.scale(beta);

  for (ordinal_type i=0; i<operators.size(); i++) {
    Teuchos::RCP<operator_type> op = operators[i];
    if (trans) {
      op_input.reshape(input.numRows(), op->coeff_size());
      op_result.reshape(result.numRows(), op->point_size());
    }
    else {
      op_input.reshape(op->coeff_size(), input.numCols());
      op_result.reshape(op->point_size(), result.numCols());
    }
    
    gather(scatter_maps[i], input, trans, op_input);
    op->transformPCE2QP(smolyak_coeffs[i], op_input, op_result, 0.0, trans);
    scatter(gather_maps[i], op_result, trans, result);
  }
}

    virtual ordinal_type ApplyInverse(
      const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Input, 
      Teuchos::SerialDenseMatrix<ordinal_type, value_type>& Result, 
      ordinal_type m) const {
      ordinal_type n=Input.numRows();
      Teuchos::SerialDenseMatrix<ordinal_type, value_type> G(A);
      Teuchos::SerialDenseMatrix<ordinal_type, value_type> z(n,1);
      for (ordinal_type j=0; j<m; j++){
	if (j==0){  // Compute z=D-1r
	  for (ordinal_type i=0; i<n; i++)
	    z(i,0)=Input(i,0)/A(i,i);
	}
	else {
	  //Compute G=invD(-L-U)=I-inv(D)A 
	  for (ordinal_type i=0; i<n; i++){
	    for (ordinal_type j=0; j<n; j++){
	      if (j==i)
		G(i,j)=0;
	      else 
		G(i,j)=-A(i,j)/A(i,i);
	    }
	  }
	  
	  Result.assign(z);
	  //z=Gz+inv(D)r
	  Result.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,1.0, G, z, 1.0);
	  
	}
      }

      return 0;
    }

 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.");
 }

  // Update *this with alpha * A * B + beta * (*this). 
  void MvTimesMatAddMv (ScalarType alpha, const Anasazi::MultiVec<ScalarType> &A, 
                        const Teuchos::SerialDenseMatrix<int, ScalarType> &B, 
                        ScalarType beta)
  {
    
    assert (Length_ == A.GetVecLength());
    assert (B.numRows() == A.GetNumberVecs());
    assert (B.numCols() <= NumberVecs_);

    MyMultiVec* MyA;
    MyA = dynamic_cast<MyMultiVec*>(&const_cast<Anasazi::MultiVec<ScalarType> &>(A)); 
    assert(MyA!=NULL);

    if ((*this)[0] == (*MyA)[0]) {
      // If this == A, then need additional storage ...
      // This situation is a bit peculiar but it may be required by
      // certain algorithms.
      
      std::vector<ScalarType> tmp(NumberVecs_);

      for (int i = 0 ; i < Length_ ; ++i) {
        for (int v = 0; v < A.GetNumberVecs() ; ++v) {
          tmp[v] = (*MyA)(i, v);
        }

        for (int v = 0 ; v < B.numCols() ; ++v) {
          (*this)(i, v) *= beta; 
          ScalarType res = Teuchos::ScalarTraits<ScalarType>::zero();

          for (int j = 0 ; j < A.GetNumberVecs() ; ++j) {
            res +=  tmp[j] * B(j, v);
          }

          (*this)(i, v) += alpha * res;
        }
      }
    }
    else {
      for (int i = 0 ; i < Length_ ; ++i) {
        for (int v = 0 ; v < B.numCols() ; ++v) {
          (*this)(i, v) *= beta; 
          ScalarType res = 0.0;
          for (int j = 0 ; j < A.GetNumberVecs() ; ++j) {
            res +=  (*MyA)(i, j) * B(j, v);
          }

          (*this)(i, v) += alpha * res;
        }
      }
    }
  }

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.");
  }
}

        // Generic BLAS level 3 matrix multiplication
        // \f$\text{this}\leftarrow \alpha A B+\beta\text{this}\f$   
        void gemm(const Real alpha,
                  const MV& A,
                  const Teuchos::SerialDenseMatrix<int,Real> &B,
                  const Real beta) {

           // Scale this by beta
            this->scale(beta);

            for(int i=0;i<B.numRows();++i) {
                for(int j=0;j<B.numCols();++j) {
                    mvec_[j]->axpy(alpha*B(i,j),*A.getVector(i));  
                }
            }
        } 

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
Stokhos::SmolyakPseudoSpectralOperator<ordinal_type,value_type,point_compare_type>::
scatter(
  const Teuchos::Array<ordinal_type>& map, 
  const Teuchos::SerialDenseMatrix<ordinal_type,value_type>& input, 
  bool trans, 
  Teuchos::SerialDenseMatrix<ordinal_type,value_type>& result) const {
  if (trans) {
    for (ordinal_type j=0; j<map.size(); j++)
      for (ordinal_type i=0; i<input.numRows(); i++)
	result(i,map[j]) += input(i,j);
  }
  else {
    for (ordinal_type j=0; j<input.numCols(); j++)
      for (ordinal_type i=0; i<map.size(); i++)
	result(map[i],j) += input(i,j);
  }
}

  // Compute a dense matrix B through the matrix-matrix multiply alpha * A^H * (*this). 
  void MvTransMv (ScalarType alpha, const Anasazi::MultiVec<ScalarType>& A, 
                  Teuchos::SerialDenseMatrix< int, ScalarType >& B
#ifdef HAVE_ANASAZI_EXPERIMENTAL
                  , Anasazi::ConjType conj
#endif
                 ) const
  {
    MyMultiVec* MyA;
    MyA = dynamic_cast<MyMultiVec*>(&const_cast<Anasazi::MultiVec<ScalarType> &>(A)); 
    assert (MyA != 0);
    
    assert (A.GetVecLength() == Length_);
    assert (NumberVecs_ <= B.numCols());
    assert (A.GetNumberVecs() <= B.numRows());
    
#ifdef HAVE_ANASAZI_EXPERIMENTAL
    if (conj == Anasazi::CONJ) {
#endif
      for (int v = 0 ; v < A.GetNumberVecs() ; ++v) {
        for (int w = 0 ; w < NumberVecs_ ; ++w) {
          ScalarType value = 0.0;
          for (int i = 0 ; i < Length_ ; ++i) {
            value += Teuchos::ScalarTraits<ScalarType>::conjugate((*MyA)(i, v)) * (*this)(i, w);
          }
          B(v, w) = alpha * value;
        }
      }
#ifdef HAVE_ANASAZI_EXPERIMENTAL
    } else {
      for (int v = 0 ; v < A.GetNumberVecs() ; ++v) {
        for (int w = 0 ; w < NumberVecs_ ; ++w) {
          ScalarType value = 0.0;
          for (int i = 0 ; i < Length_ ; ++i) {
            value += (*MyA)(i, v) * (*this)(i, w);
          }
          B(v, w) = alpha * value;
        }
      }
    }
#endif
  }

  // Compute a dense matrix B through the matrix-matrix multiply alpha * A^H * (*this).
  void MvTransMv (const ScalarType alpha, const Belos::MultiVec<ScalarType>& A,
                  Teuchos::SerialDenseMatrix< int, ScalarType >& B) const
  {
    MyMultiVec* MyA;
    MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A));
    TEUCHOS_ASSERT(MyA != NULL);

    assert (A.GetGlobalLength() == Length_);
    assert (NumberVecs_ <= B.numCols());
    assert (A.GetNumberVecs() <= B.numRows());

      for (int v = 0 ; v < A.GetNumberVecs() ; ++v) {
        for (int w = 0 ; w < NumberVecs_ ; ++w) {
          ScalarType value = 0.0;
          for (int i = 0 ; i < Length_ ; ++i) {
            value += Teuchos::ScalarTraits<ScalarType>::conjugate((*MyA)(i, v)) * (*this)(i, w);
          }
          B(v, w) = alpha * 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.");
    }
  }