44 #ifndef STOKHOS_CG_DIVISION_EXPANSION_STRATEGY_HPP 45 #define STOKHOS_CG_DIVISION_EXPANSION_STRATEGY_HPP 57 #include "Teuchos_TimeMonitor.hpp" 58 #include "Teuchos_RCP.hpp" 59 #include "Teuchos_SerialDenseMatrix.hpp" 60 #include "Teuchos_BLAS.hpp" 61 #include "Teuchos_LAPACK.hpp" 71 template <
typename ordinal_type,
typename value_type,
typename node_type>
110 const Teuchos::SerialDenseMatrix<ordinal_type, value_type> &
A,
111 Teuchos::SerialDenseMatrix<ordinal_type,value_type> &
X,
112 const Teuchos::SerialDenseMatrix<ordinal_type,value_type> &
B,
119 const Teuchos::SerialDenseMatrix<ordinal_type, value_type> &
M,
125 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type, value_type> >
basis;
131 Teuchos::RCP<const Cijk_type>
Cijk;
134 Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type,value_type> >
A,
X,
B,
M;
155 template <
typename ordinal_type,
typename value_type,
typename node_type>
169 prec_iter(prec_iter_),
179 A = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(
181 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(
183 X = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(
185 M = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(
191 template <
typename ordinal_type,
typename value_type,
typename node_type>
200 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 201 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::CGDivisionStrategy::divide()");
229 if (pb < Cijk->num_k())
230 k_end = Cijk->find_k(pb);
236 j_it != Cijk->j_end(k_it); ++j_it) {
239 i_it != Cijk->i_end(j_it); ++i_it) {
242 (*A)(i,
j) +=
cijk*cb[k];
250 (*B)(i,0) = ca[i]*basis->norm_squared(i);
252 Teuchos::SerialDenseMatrix<ordinal_type,value_type> D(sz, 1);
257 Teuchos::SerialDenseMatrix<ordinal_type, value_type> r(Teuchos::View, *A, 1, sz, i, 0);
258 D(i,0)=
sqrt(r.normOne());
265 (*A)(i,
j)=(*A)(i,
j)/(D(i,0)*D(
j,0));
271 (*B)(i,0)=(*B)(i,0)/D(i,0);
279 pb = basis->dimension()+1;
282 if (pb < Cijk->num_k())
283 k_end = Cijk->find_k(pb);
287 j_it != Cijk->j_end(k_it); ++j_it) {
290 i_it != Cijk->i_end(j_it); ++i_it) {
293 (*M)(i,
j) +=
cijk*cb[k];
303 (*M)(i,
j)=(*M)(i,
j)/(D(i,0)*D(
j,0));
307 CG(*A,*X,*B, max_it, tol, prec_iter, basis->order(), basis->dimension(), PrecNum, *M,
diag);
312 CG(*A,*X,*B, max_it, tol, prec_iter, basis->order(), basis->dimension(), PrecNum, *A,
diag);
318 (*X)(i,0)=(*X)(i,0)/D(i,0);
324 cc[i] = alpha*(*X)(i,0) + beta*cc[i];
328 cc[i] = alpha*ca[i]/cb[0] + beta*cc[i];
333 template <
typename ordinal_type,
typename value_type,
typename node_type>
336 CG(
const Teuchos::SerialDenseMatrix<ordinal_type, value_type> & A,
337 Teuchos::SerialDenseMatrix<ordinal_type,value_type> & X,
338 const Teuchos::SerialDenseMatrix<ordinal_type,value_type> & B,
345 const Teuchos::SerialDenseMatrix<ordinal_type, value_type> & M,
352 Teuchos::SerialDenseMatrix<ordinal_type, value_type>
Ax(n,1);
353 Ax.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,1.0, A, X, 0.0);
354 Teuchos::SerialDenseMatrix<ordinal_type, value_type> r(Teuchos::Copy,B);
356 resid=r.normFrobenius();
357 Teuchos::SerialDenseMatrix<ordinal_type, value_type> p(r);
358 Teuchos::SerialDenseMatrix<ordinal_type, value_type> rho(1,1);
359 Teuchos::SerialDenseMatrix<ordinal_type, value_type> oldrho(1,1);
360 Teuchos::SerialDenseMatrix<ordinal_type, value_type> pAp(1,1);
361 Teuchos::SerialDenseMatrix<ordinal_type, value_type>
Ap(n,1);
364 while (resid > tolerance && k < max_iter){
365 Teuchos::SerialDenseMatrix<ordinal_type, value_type> z(r);
372 else if (PrecNum == 2){
376 else if (PrecNum == 3){
380 else if (PrecNum == 4){
385 rho.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,1.0, r, z, 0.0);
390 rho.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, r, z, 0.0);
393 b=rho(0,0)/oldrho(0,0);
397 Ap.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,1.0, A, p, 0.0);
398 pAp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,1.0, p,
Ap, 0.0);
400 Teuchos::SerialDenseMatrix<ordinal_type, value_type> scalep(p);
406 resid=r.normFrobenius();
414 #endif // STOKHOS_DIVISION_EXPANSION_STRATEGY_HPP KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
virtual ordinal_type ApplyInverse(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Input, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Result, ordinal_type m) const
Returns the result of a Operator inverse applied to a Teuchos::SerialDenseMatrix Input in Result...
void resize(ordinal_type sz)
Resize coefficient array (coefficients are preserved)
Strategy interface for computing PCE of a/b using only b[0].
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
virtual ordinal_type ApplyInverse(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Input, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Result, ordinal_type m) const
Returns the result of a Operator inverse applied to a Teuchos::SerialDenseMatrix Input in Result...
Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > M
Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > X
pointer coeff()
Return coefficient array.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis
Basis.
Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > A
Dense matrices for linear system.
ordinal_type size() const
Return size.
const IndexType const IndexType const IndexType const IndexType const ValueType * Ax
Bi-directional iterator for traversing a sparse array.
ordinal_type prec_iter
Tolerance for CG.
Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > B
Abstract base class for multivariate orthogonal polynomials.
virtual ~CGDivisionExpansionStrategy()
Destructor.
Top-level namespace for Stokhos classes and functions.
Strategy interface for computing PCE of a/b.
virtual ordinal_type ApplyInverse(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Input, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Result, ordinal_type m) const
Returns the result of a Operator inverse applied to a Teuchos::SerialDenseMatrix Input in Result...
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< view_type >::value, typename CijkType< view_type >::type >::type cijk(const view_type &view)
const IndexType const IndexType const IndexType * Ap
Class to store coefficients of a projection onto an orthogonal polynomial basis.
CGDivisionExpansionStrategy & operator=(const CGDivisionExpansionStrategy &b)
Teuchos::RCP< const Cijk_type > Cijk
Triple product.
CGDivisionExpansionStrategy(const Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > &basis_, const Teuchos::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > &Cijk_, const ordinal_type prec_iter_, const value_type tol_, const ordinal_type PrecNum_, const ordinal_type max_it_, const ordinal_type linear_, const ordinal_type diag_, const ordinal_type equil_)
Constructor.
Stokhos::Sparse3Tensor< ordinal_type, value_type > Cijk_type
Short-hand for Cijk.
ordinal_type CG(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &A, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &X, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &B, ordinal_type max_iter, value_type tolerance, ordinal_type prec_iter, ordinal_type order, ordinal_type dim, ordinal_type PrecNum, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &M, ordinal_type diag)
virtual void divide(Stokhos::OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &alpha, const Stokhos::OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const Stokhos::OrthogPolyApprox< ordinal_type, value_type, node_type > &b, const value_type &beta)
virtual ordinal_type ApplyInverse(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Input, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Result, ordinal_type prec_iters) const
Returns the result of a Operator inverse applied to a Teuchos::SerialDenseMatrix Input in Result...