This file is indexed.

/usr/include/trilinos/Stokhos_TotalOrderBasis.hpp is in libtrilinos-stokhos-dev 12.12.1-5.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
// @HEADER
// ***********************************************************************
// 
//                           Stokhos Package
//                 Copyright (2009) Sandia Corporation
// 
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
// 
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
// 
// ***********************************************************************
// @HEADER

#ifndef STOKHOS_TOTAL_ORDER_BASIS_HPP
#define STOKHOS_TOTAL_ORDER_BASIS_HPP

#include "Teuchos_RCP.hpp"

#include "Stokhos_ProductBasis.hpp"
#include "Stokhos_OneDOrthogPolyBasis.hpp"
#include "Stokhos_ProductBasisUtils.hpp"

namespace Stokhos {

  /*!
   * \brief Multivariate orthogonal polynomial basis generated from a
   * total order tensor product of univariate polynomials.
   */
  /*!
   * The multivariate polynomials are given by 
   * \f[
   *     \Psi_i(x) = \psi_{i_1}(x_1)\dots\psi_{i_d}(x_d)
   * \f]
   * where \f$d\f$ is the dimension of the basis and \f$i_j\leq p_j\f$,
   * where \f$p_j\f$ is the order of the $j$th basis.  
   */
  template <typename ordinal_type, typename value_type,
	    typename coeff_compare_type = 
	    TotalOrderLess<MultiIndex<ordinal_type> > >
  class TotalOrderBasis : 
    public ProductBasis<ordinal_type,value_type> {
  public:

    //! Constructor
    /*!
     * \param bases array of 1-D coordinate bases
     * \param sparse_tol tolerance used to drop terms in sparse triple-product
     *                   tensors
     */
    TotalOrderBasis(
      const Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type,
 value_type> > >& bases,
      const value_type& sparse_tol = 1.0e-12,
      const coeff_compare_type& coeff_compare = coeff_compare_type());

    //! Destructor
    virtual ~TotalOrderBasis();

    //! \name Implementation of Stokhos::OrthogPolyBasis methods
    //@{

    //! Return order of basis
    ordinal_type order() const;

    //! Return dimension of basis
    ordinal_type dimension() const;

    //! Return total size of basis
    virtual ordinal_type size() const;

    //! Return array storing norm-squared of each basis polynomial
    /*!
     * Entry \f$l\f$ of returned array is given by \f$\langle\Psi_l^2\rangle\f$
     * for \f$l=0,\dots,P\f$ where \f$P\f$ is size()-1.
     */
    virtual const Teuchos::Array<value_type>& norm_squared() const;

    //! Return norm squared of basis polynomial \c i.
    virtual const value_type& norm_squared(ordinal_type i) const;

    //! Compute triple product tensor
    /*!
     * The \f$(i,j,k)\f$ entry of the tensor \f$C_{ijk}\f$ is given by
     * \f$C_{ijk} = \langle\Psi_i\Psi_j\Psi_k\rangle\f$ where \f$\Psi_l\f$
     * represents basis polynomial \f$l\f$ and \f$i,j,k=0,\dots,P\f$ where
     * \f$P\f$ is size()-1.
     */
    virtual 
    Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > 
    computeTripleProductTensor() const;

    //! Compute linear triple product tensor where k = 0,1,..,d
    virtual 
    Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > 
    computeLinearTripleProductTensor() const;

    //! Evaluate basis polynomial \c i at zero
    virtual value_type evaluateZero(ordinal_type i) const;

    //! Evaluate basis polynomials at given point \c point
    /*!
     * Size of returned array is given by size(), and coefficients are
     * ordered from order 0 up to size size()-1.
     */
    virtual void evaluateBases(
      const Teuchos::ArrayView<const value_type>& point,
      Teuchos::Array<value_type>& basis_vals) const;

    //! Print basis to stream \c os
    virtual void print(std::ostream& os) const;

    //! Return string name of basis
    virtual const std::string& getName() const;

    //@}

    //! \name Implementation of Stokhos::ProductBasis methods
    //@{

    //! Get orders of each coordinate polynomial given an index \c i
    /*!
     * The returned array is of size \f$d\f$, where \f$d\f$ is the dimension of
     * the basis, and entry \f$l\f$ is given by \f$i_l\f$ where
     * \f$\Psi_i(x) = \psi_{i_1}(x_1)\dots\psi_{i_d}(x_d)\f$.
     */
    virtual const MultiIndex<ordinal_type>& term(ordinal_type i) const;

    //! Get index of the multivariate polynomial given orders of each coordinate
    /*!
     * Given the array \c term storing \f$i_1,\dots,\i_d\f$, returns the index
     * \f$i\f$ such that \f$\Psi_i(x) = \psi_{i_1}(x_1)\dots\psi_{i_d}(x_d)\f$.
     */
    virtual ordinal_type index(const MultiIndex<ordinal_type>& term) const;

    //! Return coordinate bases
    /*!
     * Array is of size dimension().
     */
    Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type, 
							   value_type> > > 
    getCoordinateBases() const;

    //! Return maximum order allowable for each coordinate basis
    virtual MultiIndex<ordinal_type> getMaxOrders() const;

    //@}

  private:

    // Prohibit copying
    TotalOrderBasis(const TotalOrderBasis&);

    // Prohibit Assignment
    TotalOrderBasis& operator=(const TotalOrderBasis& b);
    
  protected:

    typedef MultiIndex<ordinal_type> coeff_type;
    typedef std::map<coeff_type,ordinal_type,coeff_compare_type> coeff_set_type;
    typedef Teuchos::Array<coeff_type> coeff_map_type;

    //! Name of basis
    std::string name;

    //! Total order of basis
    ordinal_type p;

    //! Total dimension of basis
    ordinal_type d;

    //! Total size of basis
    ordinal_type sz;

    //! Array of bases
    Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type, value_type> > > bases;

    //! Tolerance for computing sparse Cijk
    value_type sparse_tol;

    //! Maximum orders for each dimension
    coeff_type max_orders;

    //! Basis set
    coeff_set_type basis_set;

    //! Basis map
    coeff_map_type basis_map;

    //! Norms
    Teuchos::Array<value_type> norms;

    //! Temporary array used in basis evaluation
    mutable Teuchos::Array< Teuchos::Array<value_type> > basis_eval_tmp;

  }; // class TotalOrderBasis

  // An approach for building a sparse 3-tensor only for lexicographically
  // ordered total order basis
  // To-do:
  //   * Remove the n_choose_k() calls
  //   * Remove the loops in the Cijk_1D_Iterator::increment() functions
  //   * Store the 1-D Cijk tensors in a compressed format and eliminate
  //     the implicit searches with getValue()
  //     * Instead of looping over (i,j,k) multi-indices we could just store
  //       the 1-D Cijk tensors as an array of (i,j,k,c) tuples.
  template <typename ordinal_type, 
	    typename value_type>
  Teuchos::RCP< Sparse3Tensor<ordinal_type, value_type> >
  computeTripleProductTensorLTO(
    const TotalOrderBasis<ordinal_type, value_type,LexographicLess<MultiIndex<ordinal_type> > >& product_basis,
    bool symmetric = false) {
#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
    TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total Triple-Product Tensor Time");
#endif
    using Teuchos::RCP;
    using Teuchos::rcp;
    using Teuchos::Array;

    typedef MultiIndex<ordinal_type> coeff_type;
    const Array< RCP<const OneDOrthogPolyBasis<ordinal_type, value_type> > >& bases = product_basis.getCoordinateBases();
    ordinal_type d = bases.size();
    //ordinal_type p = product_basis.order();
    Array<ordinal_type> basis_orders(d);
    for (int i=0; i<d; ++i)
      basis_orders[i] = bases[i]->order();

    // Create 1-D triple products
    Array< RCP<Sparse3Tensor<ordinal_type,value_type> > > Cijk_1d(d);
    for (ordinal_type i=0; i<d; i++) {
      Cijk_1d[i] = 
	bases[i]->computeSparseTripleProductTensor(bases[i]->order()+1);
    }

    RCP< Sparse3Tensor<ordinal_type, value_type> > Cijk = 
      rcp(new Sparse3Tensor<ordinal_type, value_type>);
      
    // Create i, j, k iterators for each dimension
    typedef ProductBasisUtils::Cijk_1D_Iterator<ordinal_type> Cijk_Iterator;
    Array<Cijk_Iterator> Cijk_1d_iterators(d);
    coeff_type terms_i(d,0), terms_j(d,0), terms_k(d,0);
    Array<ordinal_type> sum_i(d,0), sum_j(d,0), sum_k(d,0);
    for (ordinal_type dim=0; dim<d; dim++) {
      Cijk_1d_iterators[dim] = Cijk_Iterator(bases[dim]->order(), symmetric);
    }

    ordinal_type I = 0;
    ordinal_type J = 0;
    ordinal_type K = 0;
    ordinal_type cnt = 0;
    bool stop = false;
    while (!stop) {

      // Fill out terms from 1-D iterators
      for (ordinal_type dim=0; dim<d; ++dim) {
	terms_i[dim] = Cijk_1d_iterators[dim].i;
	terms_j[dim] = Cijk_1d_iterators[dim].j;
	terms_k[dim] = Cijk_1d_iterators[dim].k;
      }

      // Compute global I,J,K
      /*
      ordinal_type II = lexicographicMapping(terms_i, p);
      ordinal_type JJ = lexicographicMapping(terms_j, p);
      ordinal_type KK = lexicographicMapping(terms_k, p);
      if (I != II || J != JJ || K != KK) {
	std::cout << "DIFF!!!" << std::endl;
	std::cout << terms_i << ":  I = " << I << ", II = " << II << std::endl;
	std::cout << terms_j << ":  J = " << J << ", JJ = " << JJ << std::endl;
	std::cout << terms_k << ":  K = " << K << ", KK = " << KK << std::endl;
      }
      */

      // Compute triple-product value
      value_type c = value_type(1.0);
      for (ordinal_type dim=0; dim<d; dim++) {
	c *= Cijk_1d[dim]->getValue(Cijk_1d_iterators[dim].i, 
				    Cijk_1d_iterators[dim].j, 
				    Cijk_1d_iterators[dim].k);
      }

      TEUCHOS_TEST_FOR_EXCEPTION(
	std::abs(c) <= 1.0e-12,
	std::logic_error,
	"Got 0 triple product value " << c 
	<< ", I = " << I << " = " << terms_i
	<< ", J = " << J << " = " << terms_j
	<< ", K = " << K << " = " << terms_k
	<< std::endl);
      
      // Add term to global Cijk
      Cijk->add_term(I,J,K,c);
      // Cijk->add_term(I,K,J,c);
      // Cijk->add_term(J,I,K,c);
      // Cijk->add_term(J,K,I,c);
      // Cijk->add_term(K,I,J,c);
      // Cijk->add_term(K,J,I,c);
      
      // Increment iterators to the next valid term
      ordinal_type cdim = d-1;
      bool inc = true;
      while (inc && cdim >= 0) {
	ordinal_type delta_i, delta_j, delta_k;
	bool more = 
	  Cijk_1d_iterators[cdim].increment(delta_i, delta_j, delta_k);

	// Update number of terms used for computing global index
	if (cdim == d-1) {
	  I += delta_i;
	  J += delta_j;
	  K += delta_k;
	}
	else {
	  if (delta_i > 0) {
	    for (ordinal_type ii=0; ii<delta_i; ++ii)
	      I += 
		Stokhos::n_choose_k(
		  basis_orders[cdim+1]-sum_i[cdim] -
		  (Cijk_1d_iterators[cdim].i-ii)+d-cdim,
		  d-cdim-1);
	  }
	  else {
	    for (ordinal_type ii=0; ii<-delta_i; ++ii)
	      I -= 
		Stokhos::n_choose_k(
		  basis_orders[cdim+1]-sum_i[cdim] -
		  (Cijk_1d_iterators[cdim].i+ii)+d-cdim-1,
		  d-cdim-1);
	  }

	  if (delta_j > 0) {
	    for (ordinal_type jj=0; jj<delta_j; ++jj)
	      J +=  
		Stokhos::n_choose_k(
		  basis_orders[cdim+1]-sum_j[cdim] -
		  (Cijk_1d_iterators[cdim].j-jj)+d-cdim,
		  d-cdim-1);
	  }
	  else {
	    for (ordinal_type jj=0; jj<-delta_j; ++jj)
	      J -=  
		Stokhos::n_choose_k(
		  basis_orders[cdim+1]-sum_j[cdim] -
		  (Cijk_1d_iterators[cdim].j+jj)+d-cdim-1,
		  d-cdim-1);
	  }
	  
	  if (delta_k > 0) {
	    for (ordinal_type kk=0; kk<delta_k; ++kk)
	      K +=  
		Stokhos::n_choose_k(
		  basis_orders[cdim+1]-sum_k[cdim] -
		  (Cijk_1d_iterators[cdim].k-kk)+d-cdim,
		  d-cdim-1);
	  }
	  else {
	    for (ordinal_type kk=0; kk<-delta_k; ++kk)
	      K -=  
		Stokhos::n_choose_k(
		  basis_orders[cdim+1]-sum_k[cdim] -
		  (Cijk_1d_iterators[cdim].k+kk)+d-cdim-1,
		  d-cdim-1);
	  }
	}
	
	if (!more) {
	  // If no more terms in this dimension, go to previous one
	  --cdim;
	}
	else {
	  // cdim has more terms, so reset iterators for all dimensions > cdim
	  // adjusting max order based on sum of i,j,k for previous dims
	  inc = false;

	  for (ordinal_type dim=cdim+1; dim<d; ++dim) {

	    // Update sums of orders for previous dimension
	    sum_i[dim] = sum_i[dim-1] + Cijk_1d_iterators[dim-1].i;
	    sum_j[dim] = sum_j[dim-1] + Cijk_1d_iterators[dim-1].j;
	    sum_k[dim] = sum_k[dim-1] + Cijk_1d_iterators[dim-1].k;

	    // Reset iterator for this dimension
	    Cijk_1d_iterators[dim] = 
	      Cijk_Iterator(basis_orders[dim]-sum_i[dim],
			    basis_orders[dim]-sum_j[dim],
			    basis_orders[dim]-sum_k[dim],
			    symmetric);
	  }
	}
      }
      
      if (cdim < 0)
	stop = true;
      
      cnt++;
    }
    
    Cijk->fillComplete();
    
    return Cijk;
  }

} // Namespace Stokhos

// Include template definitions
#include "Stokhos_TotalOrderBasisImp.hpp"

#endif