Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Stokhos_DenseDirectDivisionExpansionStrategy.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_DENSE_DIRECT_DIVISION_EXPANSION_STRATEGY_HPP
45 #define STOKHOS_DENSE_DIRECT_DIVISION_EXPANSION_STRATEGY_HPP
46 
50 
51 #include "Teuchos_TimeMonitor.hpp"
52 #include "Teuchos_RCP.hpp"
55 
56 namespace 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 
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 
97 
100 
103 
106 
109 
110  }; // class DenseDirectDivisionExpansionStrategy
111 
112 } // namespace Stokhos
113 
114 template <typename ordinal_type, typename value_type, typename node_type>
119  basis(basis_),
120  Cijk(Cijk_),
121  solver()
122 {
123  ordinal_type sz = basis->size();
125  sz, sz));
127  sz, 1));
129  sz, 1));
130 }
131 
132 template <typename ordinal_type, typename value_type, typename node_type>
133 void
136  const value_type& alpha,
139  const value_type& beta)
140 {
141 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
142  TEUCHOS_FUNC_TIME_MONITOR("Stokhos::DenseDirectDivisionStrategy::divide()");
143 #endif
144 
145  ordinal_type sz = basis->size();
146  ordinal_type pa = a.size();
147  ordinal_type pb = b.size();
148  ordinal_type pc;
149  if (pb > 1)
150  pc = sz;
151  else
152  pc = pa;
153  if (c.size() != pc)
154  c.resize(pc);
155 
156  const value_type* ca = a.coeff();
157  const value_type* cb = b.coeff();
158  value_type* cc = c.coeff();
159 
160  if (pb > 1) {
161  // Compute A
162  A->putScalar(0.0);
163  typename Cijk_type::k_iterator k_begin = Cijk->k_begin();
164  typename Cijk_type::k_iterator k_end = Cijk->k_end();
165  if (pb < Cijk->num_k())
166  k_end = Cijk->find_k(pb);
168  ordinal_type i,j,k;
169  for (typename Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
170  k = index(k_it);
171  for (typename Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
172  j_it != Cijk->j_end(k_it); ++j_it) {
173  j = index(j_it);
174  for (typename Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
175  i_it != Cijk->i_end(j_it); ++i_it) {
176  i = index(i_it);
177  cijk = value(i_it);
178  (*A)(i,j) += cijk*cb[k];
179  }
180  }
181  }
182 
183  // Compute B
184  B->putScalar(0.0);
185  for (ordinal_type i=0; i<pa; i++)
186  (*B)(i,0) = ca[i]*basis->norm_squared(i);
187 
188  // Setup solver
189  solver.setMatrix(A);
190  solver.setVectors(X, B);
191  if (solver.shouldEquilibrate()) {
192  solver.factorWithEquilibration(true);
193  solver.equilibrateMatrix();
194  }
195 
196  // Compute X = A^{-1}*B
197  solver.solve();
198 
199  // value_type kappa;
200  // solver.reciprocalConditionEstimate(kappa);
201  // std::cout << "kappa = " << 1.0/kappa << std::endl;
202 
203  // Compute c
204  for (ordinal_type i=0; i<pc; i++)
205  cc[i] = alpha*(*X)(i,0) + beta*cc[i];
206  }
207  else {
208  for (ordinal_type i=0; i<pc; i++)
209  cc[i] = alpha*ca[i]/cb[0] + beta*cc[i];
210  }
211 }
212 
213 #endif // STOKHOS_DIVISION_EXPANSION_STRATEGY_HPP
#define TEUCHOS_FUNC_TIME_MONITOR(FUNCNAME)
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)
void resize(ordinal_type sz)
Resize coefficient array (coefficients are preserved)
Teuchos::SerialDenseSolver< ordinal_type, value_type > solver
Serial dense solver.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
pointer coeff()
Return coefficient array.
Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > X
Bi-directional iterator for traversing a sparse array.
Stokhos::Sparse3Tensor< ordinal_type, value_type > Cijk_type
Short-hand for Cijk.
Abstract base class for multivariate orthogonal polynomials.
DenseDirectDivisionExpansionStrategy & operator=(const DenseDirectDivisionExpansionStrategy &b)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
DenseDirectDivisionExpansionStrategy(const Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > &basis_, const Teuchos::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > &Cijk_)
Constructor.
Strategy interface for computing PCE of a/b.
Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > B
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< view_type >::value, typename CijkType< view_type >::type >::type cijk(const view_type &view)
Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > A
Dense matrices for linear system.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis
Basis.
ordinal_type size() const
Return size.
Strategy interface for computing PCE of a/b using only b[0].