Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_GMRESDivisionExpansionStrategy.hpp
Go to the documentation of this file.
1// $Id$
2// $Source$
3// @HEADER
4// ***********************************************************************
5//
6// Stokhos Package
7// Copyright (2009) Sandia Corporation
8//
9// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10// license for use of this work by or on behalf of the U.S. Government.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
40//
41// ***********************************************************************
42// @HEADER
43
44#ifndef STOKHOS_GMRES_DIVISION_EXPANSION_STRATEGY_HPP
45#define STOKHOS_GMRES_DIVISION_EXPANSION_STRATEGY_HPP
46
55
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"
61
62namespace Stokhos {
63
65
68 template <typename ordinal_type, typename value_type, typename node_type>
70 public DivisionExpansionStrategy<ordinal_type,value_type,node_type> {
71 public:
72
75 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type, value_type> >& basis_,
76 const Teuchos::RCP<const Stokhos::Sparse3Tensor<ordinal_type, value_type> >& Cijk_,
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_);
84
87
88 // Division operation: c = \alpha*(a/b) + beta*c
89 virtual void divide(
91 const value_type& alpha,
94 const value_type& beta);
95
96 private:
97
98 // Prohibit copying
101
102 // Prohibit Assignment
105
106 ordinal_type GMRES(
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,
112 ordinal_type prec_iter,
113 ordinal_type order,
114 ordinal_type dim,
115 ordinal_type PrecNum,
116 const Teuchos::SerialDenseMatrix<ordinal_type, value_type>& M,
117 ordinal_type diag);
118
119 protected:
120
122 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type, value_type> > basis;
123
126
128 Teuchos::RCP<const Cijk_type> Cijk;
129
131 Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type,value_type> > A, X, B, M;
132
134
135 ordinal_type prec_iter;
136
137 value_type tol;
138
139 ordinal_type PrecNum;
140
141 ordinal_type max_it;
142
143 ordinal_type linear;
144
145 ordinal_type diag;
146
147 ordinal_type equil;
148
149
150 }; // class GMRESDivisionExpansionStrategy
151
152} // namespace Stokhos
153
154template <typename ordinal_type, typename value_type, typename node_type>
157 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type, value_type> >& basis_,
158 const Teuchos::RCP<const Stokhos::Sparse3Tensor<ordinal_type, value_type> >& Cijk_,
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_):
166 basis(basis_),
167 Cijk(Cijk_),
168 prec_iter(prec_iter_),
169 tol(tol_),
170 PrecNum(PrecNum_),
171 max_it(max_it_),
172 linear(linear_),
173 diag(diag_),
174 equil(equil_)
175{
176 ordinal_type sz = basis->size();
177 A = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(
178 sz, sz));
179 B = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(
180 sz, 1));
181 X = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(
182 sz, 1));
183 M = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(
184 sz, sz));
185
186}
187
188
189template <typename ordinal_type, typename value_type, typename node_type>
190void
193 const value_type& alpha,
196 const value_type& beta)
197{
198#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
199 TEUCHOS_FUNC_TIME_MONITOR("Stokhos::GMRESDivisionStrategy::divide()");
200#endif
201
202 ordinal_type sz = basis->size();
203 ordinal_type pa = a.size();
204 ordinal_type pb = b.size();
205 ordinal_type pc;
206 if (pb > 1)
207 pc = sz;
208 else
209 pc = pa;
210 if (c.size() != pc)
211 c.resize(pc);
212
213 const value_type* ca = a.coeff();
214 const value_type* cb = b.coeff();
215 value_type* cc = c.coeff();
216
217 if (pb > 1) {
218 // Compute A
219 A->putScalar(0.0);
220 typename Cijk_type::k_iterator k_begin = Cijk->k_begin();
221 typename Cijk_type::k_iterator k_end = Cijk->k_end();
222 if (pb < Cijk->num_k())
223 k_end = Cijk->find_k(pb);
224 value_type cijk;
225 ordinal_type i,j,k;
226 for (typename Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
227 k = index(k_it);
228 for (typename Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
229 j_it != Cijk->j_end(k_it); ++j_it) {
230 j = index(j_it);
231 for (typename Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
232 i_it != Cijk->i_end(j_it); ++i_it) {
233 i = index(i_it);
234 cijk = value(i_it);
235 (*A)(i,j) += cijk*cb[k];
236 }
237 }
238 }
239
240 // Compute B
241 B->putScalar(0.0);
242 for (ordinal_type i=0; i<pa; i++)
243 (*B)(i,0) = ca[i]*basis->norm_squared(i);
244
245 Teuchos::SerialDenseMatrix<ordinal_type,value_type> D(sz, 1);
246 //Equilibrate the linear system
247 if (equil == 1){
248 //Create diag mtx of max row entries
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());
252 }
253 //Compute inv(D)*A*inv(D)
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));
257 }
258 }
259
260 //Scale b by inv(D)
261 for (ordinal_type i=0; i<sz; i++){
262 (*B)(i,0)=(*B)(i,0)/D(i,0);
263 }
264
265 }
266
267 if (linear == 1){
268 //Compute M, the linear matrix to be used in the preconditioner
269 pb = basis->dimension()+1;
270 M->putScalar(0.0);
271 if (pb < Cijk->num_k())
272 k_end = Cijk->find_k(pb);
273 for (typename Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
274 k = index(k_it);
275 for ( typename Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
276 j_it != Cijk->j_end(k_it); ++j_it) {
277 j = index(j_it);
278 for ( typename Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
279 i_it != Cijk->i_end(j_it); ++i_it) {
280 i = index(i_it);
281 cijk = value(i_it);
282 (*M)(i,j) += cijk*cb[k];
283 }
284 }
285 }
286
287 //Scale M
288 if (equil == 1){
289 //Compute inv(D)*M*inv(D)
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));
293 }
294 }
295 }
296
297
298 // Compute X = A^{-1}*B
299 GMRES(*A,*X,*B, max_it, tol, prec_iter, basis->order(), basis->dimension(), PrecNum, *M, diag);
300 }
301
302 else{
303 GMRES(*A,*X,*B, max_it, tol, prec_iter, basis->order(), basis->dimension(), PrecNum, *A, diag);
304 }
305
306 if (equil == 1 ) {
307 //Rescale X
308 for (ordinal_type i=0; i<sz; i++){
309 (*X)(i,0)=(*X)(i,0)/D(i,0);
310 }
311 }
312
313 // Compute c
314 for (ordinal_type i=0; i<pc; i++)
315 cc[i] = alpha*(*X)(i,0) + beta*cc[i];
316 }
317 else {
318 for (ordinal_type i=0; i<pc; i++)
319 cc[i] = alpha*ca[i]/cb[0] + beta*cc[i];
320 }
321}
322
323
324template <typename ordinal_type, typename value_type, typename node_type>
325ordinal_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,
333 ordinal_type order,
334 ordinal_type dim,
335 ordinal_type PrecNum,
336 const Teuchos::SerialDenseMatrix<ordinal_type, value_type> & M,
337 ordinal_type diag)
338{
339 ordinal_type n = A.numRows();
340 ordinal_type k = 1;
341 value_type resid;
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);
346 r0-=Ax;
347 resid=r0.normFrobenius();
348 //define vector v=r/norm(r) where r=b-Ax
349 Teuchos::SerialDenseMatrix<ordinal_type, value_type> v(n,1);
350 r0.scale(1/resid);
351 Teuchos::SerialDenseMatrix<ordinal_type, value_type> h(1,1);
352 //Matrix of orthog basis vectors V
353 Teuchos::SerialDenseMatrix<ordinal_type, value_type> V(n,1);
354 //Set v=r0/norm(r0) to be 1st col of V
355 for (ordinal_type i=0; i<n; i++){
356 V(i,0)=r0(i,0);
357 }
358 //right hand side
359 Teuchos::SerialDenseMatrix<ordinal_type, value_type> bb(1,1);
360 bb(0,0)=resid;
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){
365 h.reshape(k+1,k);
366 //Arnoldi iteration(Gram-Schmidt )
367 V.reshape(n,k+1);
368 //set vk to be kth col of V
369 Teuchos::SerialDenseMatrix<ordinal_type, value_type> vk(Teuchos::Copy, V, n,1,0,k-1);
370 //Preconditioning step: solve Mz=vk
371 Teuchos::SerialDenseMatrix<ordinal_type, value_type> z(vk);
372 if (PrecNum == 1){
374 precond.ApplyInverse(vk,z,prec_iter);
375 }
376 else if (PrecNum == 2){
378 precond.ApplyInverse(vk,z,2);
379 }
380 else if (PrecNum == 3){
382 precond.ApplyInverse(vk,z,1);
383 }
384 else if (PrecNum == 4){
386 precond.ApplyInverse(vk,z,prec_iter);
387 }
388
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++){
393 //set vi to be ith col of V
394 Teuchos::SerialDenseMatrix<ordinal_type, value_type> vi(Teuchos::Copy, V, n,1,0,i);
395 //Calculate inner product
396 ip.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, vi, w, 0.0);
397 h(i,k-1)= ip(0,0);
398 //scale vi by h(i,k-1)
399 vi.scale(ip(0,0));
400 w-=vi;
401 }
402 h(k,k-1)=w.normFrobenius();
403 w.scale(1.0/h(k,k-1));
404 //add column vk+1=w to V
405 for (ordinal_type i=0; i<n; i++){
406 V(i,k)=w(i,0);
407 }
408 //Solve upper hessenberg least squares problem via Givens rotations
409 //Compute previous Givens rotations
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);
413 h(i,k-1)=q;
414
415 }
416 //Compute next Givens rotations
417 c.reshape(k,1);
418 s.reshape(k,1);
419 bb.reshape(k+1,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;
422 s(k-1,0)=h(k,k-1)/l;
423
424 // Givens rotation on h and bb
425 h(k-1,k-1)=l;
426 h(k,k-1)=0;
427
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);
430
431 //Determine residual
432 resid = fabs(bb(k,0));
433 k++;
434 }
435 //Extract upper triangular square matrix
436 bb.reshape(h.numRows()-1 ,1);
437 //Solve linear system
438 ordinal_type info;
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);
442 V.reshape(n,k-1);
443 ans.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, V, bb, 0.0);
444 if (PrecNum == 1){
446 precond.ApplyInverse(ans,ans,prec_iter);
447 }
448 else if (PrecNum == 2){
450 precond.ApplyInverse(ans,ans,2);
451 }
452 else if (PrecNum == 3){
454 precond.ApplyInverse(ans,ans,1);
455 }
456 else if (PrecNum == 4){
458 precond.ApplyInverse(ans,ans,prec_iter);}
459 X+=ans;
460
461 //std::cout << "iteration count= " << k-1 << std::endl;
462
463 return 0;
464}
465
466#endif // STOKHOS_DIVISION_EXPANSION_STRATEGY_HPP
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].
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 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.