Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_CGDivisionExpansionStrategy.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_CG_DIVISION_EXPANSION_STRATEGY_HPP
45#define STOKHOS_CG_DIVISION_EXPANSION_STRATEGY_HPP
46
56
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"
62
63#include <iostream>
64
65namespace Stokhos {
66
68
71 template <typename ordinal_type, typename value_type, typename node_type>
73 public DivisionExpansionStrategy<ordinal_type,value_type,node_type> {
74 public:
75
78 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type, value_type> >& basis_,
79 const Teuchos::RCP<const Stokhos::Sparse3Tensor<ordinal_type, value_type> >& Cijk_,
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_);
87
90
91 // Division operation: c = \alpha*(a/b) + beta*c
92 virtual void divide(
94 const value_type& alpha,
97 const value_type& beta);
98
99 private:
100
101 // Prohibit copying
104
105 // Prohibit Assignment
108
109 ordinal_type CG(
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,
115 ordinal_type prec_iter,
116 ordinal_type order,
117 ordinal_type dim,
118 ordinal_type PrecNum,
119 const Teuchos::SerialDenseMatrix<ordinal_type, value_type> & M,
120 ordinal_type diag);
121
122 protected:
123
125 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type, value_type> > basis;
126
129
131 Teuchos::RCP<const Cijk_type> Cijk;
132
134 Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type,value_type> > A, X, B, M;
135
137 ordinal_type prec_iter;
138
139 value_type tol;
140
141 ordinal_type PrecNum;
142
143 ordinal_type max_it;
144
145 ordinal_type linear;
146
147 ordinal_type diag;
148
149 ordinal_type equil;
150
151 }; // class CGDivisionExpansionStrategy
152
153} // namespace Stokhos
154
155template <typename ordinal_type, typename value_type, typename node_type>
158 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type, value_type> >& basis_,
159 const Teuchos::RCP<const Stokhos::Sparse3Tensor<ordinal_type, value_type> >& Cijk_,
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_):
167 basis(basis_),
168 Cijk(Cijk_),
169 prec_iter(prec_iter_),
170 tol(tol_),
171 PrecNum(PrecNum_),
172 max_it(max_it_),
173 linear(linear_),
174 diag(diag_),
175 equil(equil_)
176
177{
178 ordinal_type sz = basis->size();
179 A = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(
180 sz, sz));
181 B = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(
182 sz, 1));
183 X = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(
184 sz, 1));
185 M = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(
186 sz, sz));
187
188}
189
190
191template <typename ordinal_type, typename value_type, typename node_type>
192void
195 const value_type& alpha,
198 const value_type& beta)
199{
200#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
201 TEUCHOS_FUNC_TIME_MONITOR("Stokhos::CGDivisionStrategy::divide()");
202#endif
203
204
205 ordinal_type sz = basis->size();
206 ordinal_type pa = a.size();
207 ordinal_type pb = b.size();
208
209 ordinal_type pc;
210 if (pb > 1)
211 pc = sz;
212 else
213 pc = pa;
214 if (c.size() != pc)
215 c.resize(pc);
216
217 const value_type* ca = a.coeff();
218 const value_type* cb = b.coeff();
219
220
221 value_type* cc = c.coeff();
222
223 if (pb > 1) {
224 // Compute A
225 A->putScalar(0.0);
226 typename Cijk_type::k_iterator k_begin = Cijk->k_begin();
227 typename Cijk_type::k_iterator k_end = Cijk->k_end();
228
229 if (pb < Cijk->num_k())
230 k_end = Cijk->find_k(pb);
231 value_type cijk;
232 ordinal_type i,j,k;
233 for (typename Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
234 k = index(k_it);
235 for (typename Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
236 j_it != Cijk->j_end(k_it); ++j_it) {
237 j = index(j_it);
238 for (typename Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
239 i_it != Cijk->i_end(j_it); ++i_it) {
240 i = index(i_it);
241 cijk = value(i_it);
242 (*A)(i,j) += cijk*cb[k];
243 }
244 }
245 }
246
247 // Compute B
248 B->putScalar(0.0);
249 for (ordinal_type i=0; i<pa; i++)
250 (*B)(i,0) = ca[i]*basis->norm_squared(i);
251
252 Teuchos::SerialDenseMatrix<ordinal_type,value_type> D(sz, 1);
253 //Equilibrate the linear system
254 if (equil == 1){
255 //Create diag mtx of max row entries
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());
259 }
260
261
262 //Compute inv(D)*A*inv(D)
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));
266 }
267 }
268
269 //Scale b by inv(D)
270 for (ordinal_type i=0; i<sz; i++){
271 (*B)(i,0)=(*B)(i,0)/D(i,0);
272 }
273
274 }
275
276 if (linear == 1){
277 //Compute M, the linear matrix to be used in the preconditioner
278
279 pb = basis->dimension()+1;
280
281 M->putScalar(0.0);
282 if (pb < Cijk->num_k())
283 k_end = Cijk->find_k(pb);
284 for (typename Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
285 k = index(k_it);
286 for ( typename Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
287 j_it != Cijk->j_end(k_it); ++j_it) {
288 j = index(j_it);
289 for ( typename Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
290 i_it != Cijk->i_end(j_it); ++i_it) {
291 i = index(i_it);
292 cijk = value(i_it);
293 (*M)(i,j) += cijk*cb[k];
294 }
295 }
296 }
297
298 //Scale M
299 if (equil == 1){
300 //Compute inv(D)*M*inv(D)
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));
304 }
305 }
306 }
307 CG(*A,*X,*B, max_it, tol, prec_iter, basis->order(), basis->dimension(), PrecNum, *M, diag);
308 }
309
310 else{
311
312 CG(*A,*X,*B, max_it, tol, prec_iter, basis->order(), basis->dimension(), PrecNum, *A, diag);
313 }
314
315 if (equil == 1 ) {
316 //Rescale X
317 for (ordinal_type i=0; i<sz; i++){
318 (*X)(i,0)=(*X)(i,0)/D(i,0);
319 }
320 }
321
322 // Compute c
323 for (ordinal_type i=0; i<pc; i++)
324 cc[i] = alpha*(*X)(i,0) + beta*cc[i];
325 }
326 else {
327 for (ordinal_type i=0; i<pc; i++)
328 cc[i] = alpha*ca[i]/cb[0] + beta*cc[i];
329 }
330}
331
332
333template <typename ordinal_type, typename value_type, typename node_type>
334ordinal_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,
342 ordinal_type order ,
343 ordinal_type m,
344 ordinal_type PrecNum,
345 const Teuchos::SerialDenseMatrix<ordinal_type, value_type> & M,
346 ordinal_type diag)
347
348{
349 ordinal_type n = A.numRows();
350 ordinal_type k=0;
351 value_type resid;
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);
355 r-=Ax;
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);
362 value_type b;
363 value_type a;
364 while (resid > tolerance && k < max_iter){
365 Teuchos::SerialDenseMatrix<ordinal_type, value_type> z(r);
366 //Solve Mz=r
367 if (PrecNum != 0){
368 if (PrecNum == 1){
370 precond.ApplyInverse(r,z,prec_iter);
371 }
372 else if (PrecNum == 2){
374 precond.ApplyInverse(r,z,2);
375 }
376 else if (PrecNum == 3){
378 precond.ApplyInverse(r,z,1);
379 }
380 else if (PrecNum == 4){
382 precond.ApplyInverse(r,z,prec_iter);
383 }
384 }
385 rho.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,1.0, r, z, 0.0);
386
387
388 if (k==0){
389 p.assign(z);
390 rho.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, r, z, 0.0);
391 }
392 else {
393 b=rho(0,0)/oldrho(0,0);
394 p.scale(b);
395 p+=z;
396 }
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);
399 a=rho(0,0)/pAp(0,0);
400 Teuchos::SerialDenseMatrix<ordinal_type, value_type> scalep(p);
401 scalep.scale(a);
402 X+=scalep;
403 Ap.scale(a);
404 r-=Ap;
405 oldrho.assign(rho);
406 resid=r.normFrobenius();
407 k++;
408 }
409
410 //std::cout << "iteration count " << k << std::endl;
411 return 0;
412}
413
414 #endif // STOKHOS_DIVISION_EXPANSION_STRATEGY_HPP
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.
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.