1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43
44 #ifndef ROL_QUADRATUREDEF_HPP
45 #define ROL_QUADRATUREDEF_HPP
46
47 /** \file ROL_QuadratureDef.hpp
48 \brief Definition file for the ROL::Quadrature class.
49 \author Created by D. Kouri and D. Ridzal.
50 */
51
52 #include "ROL_QuadratureHelpers.hpp"
53
54 namespace ROL {
55
56 template <class Real>
Quadrature(int dimension)57 Quadrature<Real>::Quadrature(int dimension) : dimension_(dimension) {
58 points_.clear(); weights_.clear(); accuracy_.clear();
59 accuracy_.resize(dimension);
60 }
61
62 template <class Real>
addPointAndWeight(std::vector<Real> point,Real weight,int loc)63 void Quadrature<Real>::addPointAndWeight(std::vector<Real> point, Real weight, int loc) {
64 this->points_.insert(std::pair<std::vector<Real>,int>(point,loc));
65 this->weights_.push_back(weight);
66 }
67
68 template <class Real>
getNumPoints() const69 int Quadrature<Real>::getNumPoints() const {
70 return this->weights_.size();
71 }
72
73 template <class Real>
getAccuracy(std::vector<int> & accuracy) const74 void Quadrature<Real>::getAccuracy(std::vector<int> & accuracy) const {
75 accuracy = this->accuracy_;
76 }
77
78 template <class Real>
getCubature(std::vector<std::vector<Real>> & points,std::vector<Real> & weights) const79 void Quadrature<Real>::getCubature(std::vector<std::vector<Real> >& points,
80 std::vector<Real>& weights) const {
81 points.resize(this->weights_.size());
82 weights.resize(this->weights_.size());
83 typename std::map<std::vector<Real>,int>::const_iterator it;
84 for (it=(this->points_).begin(); it!=(this->points_).end();it++) {
85 points[it->second] = it->first;
86 weights[it->second] = (this->weights_)[it->second];
87 }
88 }
89
90 template <class Real>
getDimension() const91 int Quadrature<Real>::getDimension() const {
92 return this->dimension_;
93 }
94
95 template <class Real>
begin()96 typename std::map<std::vector<Real>,int>::iterator Quadrature<Real>::begin() {
97 return (this->points_).begin();
98 }
99
100 template <class Real>
end()101 typename std::map<std::vector<Real>,int>::iterator Quadrature<Real>::end() {
102 return (this->points_).end();
103 }
104
105 template <class Real>
insert(typename std::map<std::vector<Real>,int>::iterator it,std::vector<Real> point,Real weight)106 void Quadrature<Real>::insert(typename std::map<std::vector<Real>,int>::iterator it,
107 std::vector<Real> point, Real weight) {
108 (this->points_).insert(it,std::pair<std::vector<Real>,int>(point,(int)(this->points_).size()));
109 (this->weights_).push_back(weight);
110 }
111
112 template <class Real>
getNode(typename std::map<std::vector<Real>,int>::iterator it)113 std::vector<Real> Quadrature<Real>::getNode(typename std::map<std::vector<Real>,int>::iterator it) {
114 return it->first;
115 }
116
117 template <class Real>
getWeight(int node)118 Real Quadrature<Real>::getWeight(int node) {
119 return (this->weights_)[node];
120 }
121
122 template <class Real>
getWeight(std::vector<Real> point)123 Real Quadrature<Real>::getWeight(std::vector<Real> point) {
124 return (this->weights_)[(this->points_)[point]];
125 }
126
127 template <class Real>
update(Real alpha,Quadrature<Real> & rule)128 void Quadrature<Real>::update(Real alpha, Quadrature<Real> &rule) {
129 if ( rule.begin() != rule.end() ) {
130 // Initialize an iterator on std::map<std::vector<Real>,Real>
131 typename std::map<std::vector<Real>,int>::iterator it;
132 typename std::map<std::vector<Real>,int>::iterator it2;
133
134 // Get location of last point and weight
135 int loc = points_.size();
136
137 // Intersection of rule1 and rule2 and set difference rule2 \ rule1
138 for ( it2=rule.begin(); it2!=rule.end(); ++it2 ) {
139 it = points_.find(it2->first);
140 if ( it != points_.end() ) {
141 weights_[it->second] += alpha*rule.getWeight(it2->second);
142 }
143 else {
144 points_.insert(std::pair<std::vector<Real>, int>(it2->first,loc));
145 weights_.push_back(alpha*rule.getWeight(it2->second));
146 loc++;
147 }
148 }
149 }
150 }
151
152 template <class Real>
normalize()153 void Quadrature<Real>::normalize(){
154 Real sum = 0.0;
155 typename std::vector<Real>::iterator it;
156 for (it=weights_.begin(); it!=weights_.end(); it++) {
157 sum += *it;
158 }
159 for (it=weights_.begin(); it!=weights_.end(); it++) {
160 *it /= sum;
161 }
162 }
163
164 // PRIVATE MEMBER FUNCTION DEFINITIONS
165 template<class Real>
buildInitial(const int dimension,const int maxNumPoints,const std::vector<EQuadrature> & rule1D,const std::vector<EGrowth> & growth1D,const bool isNormalized)166 void Quadrature<Real>::buildInitial(const int dimension,
167 const int maxNumPoints,
168 const std::vector<EQuadrature> &rule1D,
169 const std::vector<EGrowth> &growth1D,
170 const bool isNormalized) {
171 accuracy_.clear();
172 accuracy_.resize(dimension);
173 std::vector<int> degree(1);
174 Quadrature<Real> newRule(1);
175 for (int i=0; i<dimension; i++) {
176 // Compute 1D rules
177 int numPoints = growthRule1D(maxNumPoints,growth1D[i],rule1D[i]);
178 Quadrature<Real> rule1(rule1D[i],numPoints,isNormalized);
179 rule1.getAccuracy(degree);
180 accuracy_.push_back(degree[0]);
181 // Build Tensor Rule
182 newRule = kron_prod<Real>(newRule,rule1);
183 }
184 typename std::map<std::vector<Real>,int>::iterator it;
185 int loc = 0;
186 std::vector<Real> node;
187 for (it=newRule.begin(); it!=newRule.end(); it++) {
188 node = it->first;
189 addPointAndWeight(node,newRule.getWeight(node),loc);
190 loc++;
191 }
192 }
193
194 template<class Real>
buildSGMGA(const int dimension,const int maxLevel,const std::vector<EQuadrature> & rule1D,const std::vector<EGrowth> & growth1D)195 void Quadrature<Real>::buildSGMGA(const int dimension,
196 const int maxLevel,
197 const std::vector<EQuadrature> &rule1D,
198 const std::vector<EGrowth> &growth1D) {
199 std::vector<Real> sparse_point;
200 std::vector<Real> sparse_weight;
201 int point_num = 0;
202
203 Real tol = std::sqrt(ROL_EPSILON<Real>());
204
205 GWPointer *gw_compute_points;
206 gw_compute_points = new GWPointer[dimension];
207 GWPointer *gw_compute_weights;
208 gw_compute_weights = new GWPointer[dimension];
209 GWPointer2 *gw_compute_order;
210 gw_compute_order = new GWPointer2[dimension];
211
212 int growth = 0;
213
214 for (int i=0; i<dimension; ++i) {
215 if(rule1D[i] == QUAD_CLENSHAWCURTIS) {
216 gw_compute_points[i] = &Quadrature<Real>::ClenshawCurtisPoints;
217 gw_compute_weights[i] = &Quadrature<Real>::ClenshawCurtisWeights;
218 gw_compute_order[i] = &webbur::level_to_order_exp_cc;
219 }
220 else if (rule1D[i] == QUAD_FEJER2) {
221 gw_compute_points[i] = &Quadrature<Real>::Fejer2Points;
222 gw_compute_weights[i] = &Quadrature<Real>::Fejer2Weights;
223 gw_compute_order[i] = &webbur::level_to_order_exp_f2;
224 }
225 else if (rule1D[i] == QUAD_LEGENDRE) {
226 gw_compute_points[i] = &Quadrature<Real>::LegendrePoints;
227 gw_compute_weights[i] = &Quadrature<Real>::LegendreWeights;
228 gw_compute_order[i] = &webbur::level_to_order_exp_gauss;
229 }
230 else if (rule1D[i] == QUAD_PATTERSON) {
231 gw_compute_points[i] = &Quadrature<Real>::PattersonPoints;
232 gw_compute_weights[i] = &Quadrature<Real>::PattersonWeights;
233 gw_compute_order[i] = &webbur::level_to_order_exp_gp;
234 }
235 else if (rule1D[i] == QUAD_HERMITE) {
236 gw_compute_points[i] = &Quadrature<Real>::HermitePoints;
237 gw_compute_weights[i] = &Quadrature<Real>::HermiteWeights;
238 gw_compute_order[i] = &webbur::level_to_order_exp_gauss;
239 }
240 else if (rule1D[i] == QUAD_GENZKEISTER) {
241 gw_compute_points[i] = &Quadrature<Real>::GenzKeisterPoints;
242 gw_compute_weights[i] = &Quadrature<Real>::GenzKeisterWeights;
243 gw_compute_order[i] = &webbur::level_to_order_exp_hgk;
244 }
245 else if (rule1D[i] == QUAD_LAGUERRE) {
246 gw_compute_points[i] = &Quadrature<Real>::LaguerrePoints;
247 gw_compute_weights[i] = &Quadrature<Real>::LaguerreWeights;
248 gw_compute_order[i] = &webbur::level_to_order_exp_gauss;
249 }
250
251 if( growth1D[i] == GROWTH_DEFAULT ||
252 growth1D[i] == GROWTH_FULLEXP ) {
253 growth = 2;
254 }
255 else if ( growth1D[i] == GROWTH_SLOWLIN ||
256 growth1D[i] == GROWTH_SLOWLINODD ||
257 growth1D[i] == GROWTH_SLOWEXP ) {
258 growth = 0;
259 }
260 else if ( growth1D[i] == GROWTH_MODLIN ||
261 growth1D[i] == GROWTH_MODEXP ) {
262 growth = 1;
263 }
264 }
265
266 std::vector<Real> level_weight(dimension,1);
267 std::vector<int> np(dimension,0);
268 int np_sum = webbur::i4vec_sum(dimension,&np[0]);
269 std::vector<Real> p(np_sum,0);
270
271 // Compute necessary data.
272 int point_total_num
273 = webbur::sandia_sgmga_size_total(dimension, &level_weight[0], maxLevel,
274 growth, gw_compute_order);
275
276 point_num
277 = webbur::sandia_sgmga_size(dimension, &level_weight[0], maxLevel,
278 gw_compute_points, tol, growth,
279 gw_compute_order);
280
281 std::vector<int> sparse_unique_index(point_total_num);
282 webbur::sandia_sgmga_unique_index(dimension, &level_weight[0], maxLevel,
283 gw_compute_points, tol, point_num,
284 point_total_num, growth, gw_compute_order,
285 &sparse_unique_index[0]);
286
287 std::vector<int> sparse_order(dimension*point_num,0);
288 std::vector<int> sparse_index(dimension*point_num,0);
289 webbur::sandia_sgmga_index(dimension, &level_weight[0], maxLevel, point_num,
290 point_total_num, &sparse_unique_index[0], growth,
291 gw_compute_order, &sparse_order[0],
292 &sparse_index[0]);
293
294 // Compute points and weights.
295 sparse_point.resize(dimension*point_num,0.0);
296 webbur::sandia_sgmga_point(dimension, &level_weight[0], maxLevel,
297 gw_compute_points, point_num, &sparse_order[0],
298 &sparse_index[0], growth, gw_compute_order,
299 &sparse_point[0]);
300
301 sparse_weight.resize(point_num,0.0);
302 webbur::sandia_sgmga_weight(dimension, &level_weight[0], maxLevel,
303 gw_compute_weights, point_num, point_total_num,
304 &sparse_unique_index[0], growth,
305 gw_compute_order, &sparse_weight[0]);
306
307 delete [] gw_compute_points;
308 delete [] gw_compute_weights;
309 delete [] gw_compute_order;
310 /***********************************************************************/
311 /******** END BUILD SPARSE GRID USING SANDIA_SGMGA *********************/
312 /***********************************************************************/
313 typename std::map<std::vector<Real>,int>::iterator it;
314 Real weight(0);
315 std::vector<Real> node(dimension);
316 for (int i = 0; i < point_num; ++i) {
317 weight = sparse_weight[i];
318 for (int j = 0; j < dimension; ++j) {
319 node[j] = sparse_point[j+i*dimension];
320 }
321 addPointAndWeight(node,weight,i);
322 }
323 }
324
325 template<class Real>
ClenshawCurtisPoints(int n,int dim,double x[])326 void Quadrature<Real>::ClenshawCurtisPoints(int n, int dim, double x[]) {
327 webbur::clenshaw_curtis_compute_points(n,x);
328 }
329
330 template<class Real>
ClenshawCurtisWeights(int n,int dim,double w[])331 void Quadrature<Real>::ClenshawCurtisWeights(int n, int dim, double w[]) {
332 webbur::clenshaw_curtis_compute_weights(n,w);
333 }
334
335 template<class Real>
Fejer2Points(int n,int dim,double x[])336 void Quadrature<Real>::Fejer2Points(int n, int dim, double x[]) {
337 webbur::fejer2_compute_points(n,x);
338 }
339
340 template<class Real>
Fejer2Weights(int n,int dim,double w[])341 void Quadrature<Real>::Fejer2Weights(int n, int dim, double w[]) {
342 webbur::fejer2_compute_weights(n,w);
343 }
344
345 template<class Real>
LegendrePoints(int n,int dim,double x[])346 void Quadrature<Real>::LegendrePoints(int n, int dim, double x[]) {
347 webbur::legendre_compute_points(n,x);
348 }
349
350 template<class Real>
LegendreWeights(int n,int dim,double w[])351 void Quadrature<Real>::LegendreWeights(int n, int dim, double w[]) {
352 webbur::legendre_compute_weights(n,w);
353 }
354
355 template<class Real>
PattersonPoints(int n,int dim,double x[])356 void Quadrature<Real>::PattersonPoints(int n, int dim, double x[]) {
357 webbur::patterson_lookup_points(n,x);
358 }
359
360 template<class Real>
PattersonWeights(int n,int dim,double w[])361 void Quadrature<Real>::PattersonWeights(int n, int dim, double w[]) {
362 webbur::patterson_lookup_weights(n,w);
363 }
364
365 template<class Real>
HermitePoints(int n,int dim,double x[])366 void Quadrature<Real>::HermitePoints(int n, int dim, double x[]) {
367 webbur::hermite_compute_points(n,x);
368 }
369
370 template<class Real>
HermiteWeights(int n,int dim,double w[])371 void Quadrature<Real>::HermiteWeights(int n, int dim, double w[]) {
372 webbur::hermite_compute_weights(n,w);
373 }
374
375 template<class Real>
GenzKeisterPoints(int n,int dim,double x[])376 void Quadrature<Real>::GenzKeisterPoints(int n, int dim, double x[]) {
377 webbur::hermite_genz_keister_lookup_points(n,x);
378 }
379
380 template<class Real>
GenzKeisterWeights(int n,int dim,double w[])381 void Quadrature<Real>::GenzKeisterWeights(int n, int dim, double w[]) {
382 webbur::hermite_genz_keister_lookup_weights(n,w);
383 }
384
385 template<class Real>
LaguerrePoints(int n,int dim,double x[])386 void Quadrature<Real>::LaguerrePoints(int n, int dim, double x[]) {
387 webbur::laguerre_compute_points(n,x);
388 }
389
390 template<class Real>
LaguerreWeights(int n,int dim,double w[])391 void Quadrature<Real>::LaguerreWeights(int n, int dim, double w[]) {
392 webbur::laguerre_compute_weights(n,w);
393 }
394
395 } // end ROL namespace
396
397 #endif
398