Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_SPDDenseDirectDivisionExpansionStrategy.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_SPD_DENSE_DIRECT_DIVISION_EXPANSION_STRATEGY_HPP
45#define STOKHOS_SPD_DENSE_DIRECT_DIVISION_EXPANSION_STRATEGY_HPP
46
50
51#include "Teuchos_TimeMonitor.hpp"
52#include "Teuchos_RCP.hpp"
53#include "Teuchos_SerialSymDenseMatrix.hpp"
54#include "Teuchos_SerialSpdDenseSolver.hpp"
55
56namespace Stokhos {
57
59
62 template <typename ordinal_type, typename value_type, typename node_type>
64 public DivisionExpansionStrategy<ordinal_type,value_type,node_type> {
65 public:
66
69 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type, value_type> >& basis_,
70 const Teuchos::RCP<const Stokhos::Sparse3Tensor<ordinal_type, value_type> >& Cijk_);
71
74
75 // Division operation: c = \alpha*(a/b) + beta*c
76 virtual void divide(
78 const value_type& alpha,
81 const value_type& beta);
82
83 private:
84
85 // Prohibit copying
88
89 // Prohibit Assignment
92
93 protected:
94
96 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type, value_type> > basis;
97
100
102 Teuchos::RCP<const Cijk_type> Cijk;
103
105 Teuchos::RCP< Teuchos::SerialSymDenseMatrix<ordinal_type,value_type> > A;
106
107 Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type,value_type> > X, B;
108
110 Teuchos::SerialSpdDenseSolver<ordinal_type,value_type> solver;
111
112 }; // class SPDDenseDirectDivisionExpansionStrategy
113
114} // namespace Stokhos
115
116template <typename ordinal_type, typename value_type, typename node_type>
119 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type, value_type> >& basis_,
120 const Teuchos::RCP<const Stokhos::Sparse3Tensor<ordinal_type, value_type> >& Cijk_) :
121 basis(basis_),
122 Cijk(Cijk_),
123 solver()
124{
125 ordinal_type sz = basis->size();
126 A = Teuchos::rcp(new Teuchos::SerialSymDenseMatrix<ordinal_type,value_type>(
127 sz, sz));
128 B = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(
129 sz, 1));
130 X = Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(
131 sz, 1));
132}
133
134template <typename ordinal_type, typename value_type, typename node_type>
135void
138 const value_type& alpha,
141 const value_type& beta)
142{
143#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
144 TEUCHOS_FUNC_TIME_MONITOR("Stokhos::SPDDenseDirectDivisionStrategy::divide()");
145#endif
146
147 ordinal_type sz = basis->size();
148 ordinal_type pa = a.size();
149 ordinal_type pb = b.size();
150 ordinal_type pc;
151 if (pb > 1)
152 pc = sz;
153 else
154 pc = pa;
155 if (c.size() != pc)
156 c.resize(pc);
157
158 const value_type* ca = a.coeff();
159 const value_type* cb = b.coeff();
160 value_type* cc = c.coeff();
161
162 if (pb > 1) {
163 // Compute A
164 A->putScalar(0.0);
165 typename Cijk_type::k_iterator k_begin = Cijk->k_begin();
166 typename Cijk_type::k_iterator k_end = Cijk->k_end();
167 if (pb < Cijk->num_k())
168 k_end = Cijk->find_k(pb);
169 value_type cijk;
170 ordinal_type i,j,k;
171 for (typename Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
172 k = index(k_it);
173 for (typename Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
174 j_it != Cijk->j_end(k_it); ++j_it) {
175 j = index(j_it);
176 for (typename Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
177 i_it != Cijk->i_end(j_it); ++i_it) {
178 i = index(i_it);
179 cijk = value(i_it);
180 (*A)(i,j) += cijk*cb[k];
181 }
182 }
183 }
184 A->setUpper();
185
186 // Compute B
187 B->putScalar(0.0);
188 for (ordinal_type i=0; i<pa; i++)
189 (*B)(i,0) = ca[i]*basis->norm_squared(i);
190
191 // Setup solver
192 solver.setMatrix(A);
193 solver.setVectors(X, B);
194 if (solver.shouldEquilibrate()) {
195 solver.factorWithEquilibration(true);
196 solver.equilibrateMatrix();
197 }
198
199 // Compute X = A^{-1}*B
200 solver.solve();
201
202 // value_type kappa;
203 // solver.reciprocalConditionEstimate(kappa);
204 // std::cout << "kappa = " << 1.0/kappa << std::endl;
205
206 // Compute c
207 for (ordinal_type i=0; i<pc; i++)
208 cc[i] = alpha*(*X)(i,0) + beta*cc[i];
209 }
210 else {
211 for (ordinal_type i=0; i<pc; i++)
212 cc[i] = alpha*ca[i]/cb[0] + beta*cc[i];
213 }
214}
215
216#endif // STOKHOS_DIVISION_EXPANSION_STRATEGY_HPP
Strategy interface for computing PCE of a/b.
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.
Strategy interface for computing PCE of a/b using only b[0].
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis
Basis.
Teuchos::SerialSpdDenseSolver< ordinal_type, value_type > solver
Serial dense solver.
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::SerialSymDenseMatrix< ordinal_type, value_type > > A
Dense matrices for linear system.
Stokhos::Sparse3Tensor< ordinal_type, value_type > Cijk_type
Short-hand for Cijk.
Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > X
SPDDenseDirectDivisionExpansionStrategy(const SPDDenseDirectDivisionExpansionStrategy &)
Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > B
SPDDenseDirectDivisionExpansionStrategy(const Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > &basis_, const Teuchos::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > &Cijk_)
Constructor.
SPDDenseDirectDivisionExpansionStrategy & operator=(const SPDDenseDirectDivisionExpansionStrategy &b)
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.
j_sparse_array::const_iterator kji_iterator
Iterator for looping over i entries given k and j.
Top-level namespace for Stokhos classes and functions.
Bi-directional iterator for traversing a sparse array.