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>
80 const ordinal_type prec_iter_,
81 const value_type tol_,
82 const ordinal_type PrecNum_,
83 const ordinal_type max_it_,
84 const ordinal_type linear_,
85 const ordinal_type diag_,
86 const ordinal_type equil_);
94 const value_type& alpha,
97 const value_type& beta);
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,
113 ordinal_type max_iter,
114 value_type tolerance,
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;
155template <
typename ordinal_type,
typename value_type,
typename node_type>
160 const ordinal_type prec_iter_,
161 const value_type tol_,
162 const ordinal_type PrecNum_,
163 const ordinal_type max_it_,
164 const ordinal_type linear_,
165 const ordinal_type diag_,
166 const ordinal_type equil_):
169 prec_iter(prec_iter_),
178 ordinal_type sz =
basis->size();
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>(
191template <
typename ordinal_type,
typename value_type,
typename node_type>
195 const value_type& alpha,
198 const value_type& beta)
200#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
201 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::CGDivisionStrategy::divide()");
205 ordinal_type sz = basis->size();
206 ordinal_type pa = a.
size();
207 ordinal_type pb = b.
size();
217 const value_type* ca = a.
coeff();
218 const value_type* cb = b.
coeff();
221 value_type* cc = c.
coeff();
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];
249 for (ordinal_type i=0; i<pa; i++)
250 (*B)(i,0) = ca[i]*basis->norm_squared(i);
252 Teuchos::SerialDenseMatrix<ordinal_type,value_type> D(sz, 1);
256 for (ordinal_type i=0; i<sz; i++){
257 Teuchos::SerialDenseMatrix<ordinal_type, value_type> r(Teuchos::View, *A, 1, sz, i, 0);
258 D(i,0)=sqrt(r.normOne());
263 for (ordinal_type i=0; i<sz; i++){
264 for (ordinal_type
j=0;
j<sz;
j++){
265 (*A)(i,
j)=(*A)(i,
j)/(D(i,0)*D(
j,0));
270 for (ordinal_type i=0; i<sz; i++){
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];
301 for (ordinal_type i=0; i<sz; i++){
302 for (ordinal_type
j=0;
j<sz;
j++){
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);
317 for (ordinal_type i=0; i<sz; i++){
318 (*X)(i,0)=(*X)(i,0)/D(i,0);
323 for (ordinal_type i=0; i<pc; i++)
324 cc[i] = alpha*(*X)(i,0) + beta*cc[i];
327 for (ordinal_type i=0; i<pc; i++)
328 cc[i] = alpha*ca[i]/cb[0] + beta*cc[i];
333template <
typename ordinal_type,
typename value_type,
typename node_type>
336CG(
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,
339 ordinal_type max_iter,
340 value_type tolerance,
341 ordinal_type prec_iter,
344 ordinal_type PrecNum,
345 const Teuchos::SerialDenseMatrix<ordinal_type, value_type> & M,
349 ordinal_type n = A.numRows();
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();
Strategy interface for computing PCE of a/b using only b[0].
Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > M
Teuchos::RCP< const Cijk_type > Cijk
Triple product.
Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > B
Stokhos::Sparse3Tensor< ordinal_type, value_type > Cijk_type
Short-hand for Cijk.
virtual ~CGDivisionExpansionStrategy()
Destructor.
ordinal_type prec_iter
Tolerance for CG.
CGDivisionExpansionStrategy(const CGDivisionExpansionStrategy &)
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.
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)
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis
Basis.
CGDivisionExpansionStrategy & operator=(const CGDivisionExpansionStrategy &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 > > X
Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > A
Dense matrices for linear system.
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.
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.
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.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
void resize(ordinal_type sz)
Resize coefficient array (coefficients are preserved)
pointer coeff()
Return coefficient array.
ordinal_type size() const
Return size.
Abstract base class for multivariate orthogonal polynomials.
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.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
kji_sparse_array::const_iterator k_iterator
Iterator for looping over k entries.
Top-level namespace for Stokhos classes and functions.
Bi-directional iterator for traversing a sparse array.