44 #ifndef STOKHOS_GMRES_DIVISION_EXPANSION_STRATEGY_HPP 45 #define STOKHOS_GMRES_DIVISION_EXPANSION_STRATEGY_HPP 56 #include "Teuchos_TimeMonitor.hpp" 57 #include "Teuchos_RCP.hpp" 58 #include "Teuchos_SerialDenseMatrix.hpp" 59 #include "Teuchos_BLAS.hpp" 60 #include "Teuchos_LAPACK.hpp" 68 template <
typename ordinal_type,
typename value_type,
typename node_type>
107 const Teuchos::SerialDenseMatrix<ordinal_type, value_type> &
A,
108 Teuchos::SerialDenseMatrix<ordinal_type,value_type> &
X,
109 const Teuchos::SerialDenseMatrix<ordinal_type,value_type> &
B,
116 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>&
M,
122 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type, value_type> >
basis;
128 Teuchos::RCP<const Cijk_type>
Cijk;
131 Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type,value_type> >
A,
X,
B,
M;
154 template <
typename ordinal_type,
typename value_type,
typename node_type>
168 prec_iter(prec_iter_),
177 A = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(
179 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(
181 X = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(
183 M = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(
189 template <
typename ordinal_type,
typename value_type,
typename node_type>
198 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 199 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::GMRESDivisionStrategy::divide()");
222 if (pb < Cijk->num_k())
223 k_end = Cijk->find_k(pb);
229 j_it != Cijk->j_end(k_it); ++j_it) {
232 i_it != Cijk->i_end(j_it); ++i_it) {
235 (*A)(i,
j) +=
cijk*cb[k];
243 (*B)(i,0) = ca[i]*basis->norm_squared(i);
245 Teuchos::SerialDenseMatrix<ordinal_type,value_type> D(sz, 1);
250 Teuchos::SerialDenseMatrix<ordinal_type, value_type> r(Teuchos::View, *A, 1, sz, i, 0);
251 D(i,0)=
sqrt(r.normOne());
256 (*A)(i,
j)=(*A)(i,
j)/(D(i,0)*D(
j,0));
262 (*B)(i,0)=(*B)(i,0)/D(i,0);
269 pb = basis->dimension()+1;
271 if (pb < Cijk->num_k())
272 k_end = Cijk->find_k(pb);
276 j_it != Cijk->j_end(k_it); ++j_it) {
279 i_it != Cijk->i_end(j_it); ++i_it) {
282 (*M)(i,
j) +=
cijk*cb[k];
292 (*M)(i,
j)=(*M)(i,
j)/(D(i,0)*D(
j,0));
299 GMRES(*A,*X,*B, max_it, tol, prec_iter, basis->order(), basis->dimension(), PrecNum, *M,
diag);
303 GMRES(*A,*X,*B, max_it, tol, prec_iter, basis->order(), basis->dimension(), PrecNum, *A,
diag);
309 (*X)(i,0)=(*X)(i,0)/D(i,0);
315 cc[i] = alpha*(*X)(i,0) + beta*cc[i];
319 cc[i] = alpha*ca[i]/cb[0] + beta*cc[i];
324 template <
typename ordinal_type,
typename value_type,
typename node_type>
327 GMRES(
const Teuchos::SerialDenseMatrix<ordinal_type, value_type> & A,
328 Teuchos::SerialDenseMatrix<ordinal_type,value_type> & X,
329 const Teuchos::SerialDenseMatrix<ordinal_type,value_type> & B,
336 const Teuchos::SerialDenseMatrix<ordinal_type, value_type> & M,
342 Teuchos::SerialDenseMatrix<ordinal_type, value_type> P(n,n);
343 Teuchos::SerialDenseMatrix<ordinal_type, value_type>
Ax(n,1);
344 Ax.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,1.0, A, X, 0.0);
345 Teuchos::SerialDenseMatrix<ordinal_type, value_type> r0(Teuchos::Copy,B);
347 resid=r0.normFrobenius();
349 Teuchos::SerialDenseMatrix<ordinal_type, value_type> v(n,1);
351 Teuchos::SerialDenseMatrix<ordinal_type, value_type> h(1,1);
353 Teuchos::SerialDenseMatrix<ordinal_type, value_type> V(n,1);
359 Teuchos::SerialDenseMatrix<ordinal_type, value_type> bb(1,1);
361 Teuchos::SerialDenseMatrix<ordinal_type, value_type> w(n,1);
362 Teuchos::SerialDenseMatrix<ordinal_type, value_type> c;
363 Teuchos::SerialDenseMatrix<ordinal_type, value_type> s;
364 while (resid > tolerance && k < max_iter){
369 Teuchos::SerialDenseMatrix<ordinal_type, value_type> vk(Teuchos::Copy, V, n,1,0,k-1);
371 Teuchos::SerialDenseMatrix<ordinal_type, value_type> z(vk);
376 else if (PrecNum == 2){
380 else if (PrecNum == 3){
384 else if (PrecNum == 4){
389 w.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1, A, z, 0.0);
390 Teuchos::SerialDenseMatrix<ordinal_type, value_type> vi(n,1);
391 Teuchos::SerialDenseMatrix<ordinal_type, value_type> ip(1,1);
394 Teuchos::SerialDenseMatrix<ordinal_type, value_type> vi(Teuchos::Copy, V, n,1,0,i);
396 ip.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, vi, w, 0.0);
402 h(k,k-1)=w.normFrobenius();
403 w.scale(1.0/h(k,k-1));
411 value_type q=c(i,0)*h(i,k-1)+s(i,0)*h(i+1,k-1);
412 h(i+1,k-1)=-1*s(i,0)*h(i,k-1)+c(i,0)*h(i+1,k-1);
421 c(k-1,0)=h(k-1,k-1)/l;
428 bb(k,0)=-s(k-1,0)*bb(k-1,0);
429 bb(k-1,0)=c(k-1,0)*bb(k-1,0);
432 resid =
fabs(bb(k,0));
436 bb.reshape(h.numRows()-1 ,1);
440 lapack.TRTRS(
'U',
'N',
'N', h.numRows()-1, 1, h.values(), h.stride(), bb.values(), bb.stride(),&info);
441 Teuchos::SerialDenseMatrix<ordinal_type, value_type> ans(X);
443 ans.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, V, bb, 0.0);
448 else if (PrecNum == 2){
452 else if (PrecNum == 3){
456 else if (PrecNum == 4){
466 #endif // STOKHOS_DIVISION_EXPANSION_STRATEGY_HPP KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > fabs(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)
Teuchos::RCP< const Cijk_type > Cijk
Triple product.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
ordinal_type prec_iter
Tolerance for GMRES.
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...
Stokhos::Sparse3Tensor< ordinal_type, value_type > Cijk_type
Short-hand for Cijk.
ordinal_type GMRES(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)
pointer coeff()
Return coefficient array.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis
Basis.
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.
Abstract base class for multivariate orthogonal polynomials.
Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > B
Top-level namespace for Stokhos classes and functions.
Strategy interface for computing PCE of a/b.
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)
Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > M
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 > > X
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)
Strategy interface for computing PCE of a/b using only b[0].
Class to store coefficients of a projection onto an orthogonal polynomial basis.
virtual ~GMRESDivisionExpansionStrategy()
Destructor.
GMRESDivisionExpansionStrategy & operator=(const GMRESDivisionExpansionStrategy &b)
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...
GMRESDivisionExpansionStrategy(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.
Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > A
Dense matrices for linear system.