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>
77 const ordinal_type prec_iter_,
78 const value_type tol_,
79 const ordinal_type PrecNum_,
80 const ordinal_type max_it_,
81 const ordinal_type linear_,
82 const ordinal_type diag_,
83 const ordinal_type equil_);
91 const value_type& alpha,
94 const value_type& beta);
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,
110 ordinal_type max_iter,
111 value_type tolerance,
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;
154template <
typename ordinal_type,
typename value_type,
typename node_type>
159 const ordinal_type prec_iter_,
160 const value_type tol_,
161 const ordinal_type PrecNum_,
162 const ordinal_type max_it_,
163 const ordinal_type linear_,
164 const ordinal_type diag_,
165 const ordinal_type equil_):
168 prec_iter(prec_iter_),
176 ordinal_type sz =
basis->size();
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>(
189template <
typename ordinal_type,
typename value_type,
typename node_type>
193 const value_type& alpha,
196 const value_type& beta)
198#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
199 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::GMRESDivisionStrategy::divide()");
202 ordinal_type sz = basis->size();
203 ordinal_type pa = a.
size();
204 ordinal_type pb = b.
size();
213 const value_type* ca = a.
coeff();
214 const value_type* cb = b.
coeff();
215 value_type* cc = c.
coeff();
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];
242 for (ordinal_type i=0; i<pa; i++)
243 (*B)(i,0) = ca[i]*basis->norm_squared(i);
245 Teuchos::SerialDenseMatrix<ordinal_type,value_type> D(sz, 1);
249 for (ordinal_type i=0; i<sz; i++){
250 Teuchos::SerialDenseMatrix<ordinal_type, value_type> r(Teuchos::View, *A, 1, sz, i, 0);
251 D(i,0)=sqrt(r.normOne());
254 for (ordinal_type i=0; i<sz; i++){
255 for (ordinal_type
j=0;
j<sz;
j++){
256 (*A)(i,
j)=(*A)(i,
j)/(D(i,0)*D(
j,0));
261 for (ordinal_type i=0; i<sz; i++){
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];
290 for (ordinal_type i=0; i<sz; i++){
291 for (ordinal_type
j=0;
j<sz;
j++){
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);
308 for (ordinal_type i=0; i<sz; i++){
309 (*X)(i,0)=(*X)(i,0)/D(i,0);
314 for (ordinal_type i=0; i<pc; i++)
315 cc[i] = alpha*(*X)(i,0) + beta*cc[i];
318 for (ordinal_type i=0; i<pc; i++)
319 cc[i] = alpha*ca[i]/cb[0] + beta*cc[i];
324template <
typename ordinal_type,
typename value_type,
typename node_type>
327GMRES(
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,
330 ordinal_type max_iter,
331 value_type tolerance,
332 ordinal_type prec_iter,
335 ordinal_type PrecNum,
336 const Teuchos::SerialDenseMatrix<ordinal_type, value_type> & M,
339 ordinal_type n = A.numRows();
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);
355 for (ordinal_type i=0; i<n; i++){
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);
392 for (ordinal_type i=0; i<k; i++){
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));
405 for (ordinal_type i=0; i<n; i++){
410 for (ordinal_type i=0; i<k-1; i++){
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);
420 value_type l = sqrt(h(k-1,k-1)*h(k-1,k-1)+h(k,k-1)*h(k,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){
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.
Strategy interface for computing PCE of a/b using only b[0].
ordinal_type prec_iter
Tolerance for GMRES.
GMRESDivisionExpansionStrategy(const GMRESDivisionExpansionStrategy &)
Stokhos::Sparse3Tensor< ordinal_type, value_type > Cijk_type
Short-hand for Cijk.
Teuchos::RCP< const Cijk_type > Cijk
Triple product.
Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > M
Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > X
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.
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)
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis
Basis.
virtual ~GMRESDivisionExpansionStrategy()
Destructor.
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 > > B
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 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.
Specialization for Sacado::UQ::PCE< Storage<...> >
Top-level namespace for Stokhos classes and functions.
Bi-directional iterator for traversing a sparse array.