Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_CompletePolynomialBasisImp.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) 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 Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41#include "Teuchos_TimeMonitor.hpp"
42
43template <typename ordinal_type, typename value_type>
46 const Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type, value_type> > >& bases_,
47 const value_type& sparse_tol_,
48 bool use_old_cijk_alg_,
49 const Teuchos::RCP< Teuchos::Array<value_type> >& deriv_coeffs_) :
50 p(0),
51 d(bases_.size()),
52 sz(0),
53 bases(bases_),
54 basis_orders(d),
55 sparse_tol(sparse_tol_),
56 use_old_cijk_alg(use_old_cijk_alg_),
57 deriv_coeffs(deriv_coeffs_),
58 norms(),
59 terms()
60{
61
62 // Compute total order
63 for (ordinal_type i=0; i<d; i++) {
64 basis_orders[i] = bases[i]->order();
65 if (basis_orders[i] > p)
66 p = basis_orders[i];
67 }
68
69 // Compute basis terms
71
72 // Compute norms
73 norms.resize(sz);
74 value_type nrm;
75 for (ordinal_type k=0; k<sz; k++) {
76 nrm = value_type(1.0);
77 for (ordinal_type i=0; i<d; i++)
78 nrm = nrm * bases[i]->norm_squared(terms[k][i]);
79 norms[k] = nrm;
80 }
81
82 // Create name
83 name = "Complete polynomial basis (";
84 for (ordinal_type i=0; i<d-1; i++)
85 name += bases[i]->getName() + ", ";
86 name += bases[d-1]->getName() + ")";
87
88 // Allocate array for basis evaluation
89 basis_eval_tmp.resize(d);
90 for (ordinal_type j=0; j<d; j++)
91 basis_eval_tmp[j].resize(basis_orders[j]+1);
92
93 // Set up deriv_coeffs
94 if (deriv_coeffs == Teuchos::null) {
95 deriv_coeffs = Teuchos::rcp(new Teuchos::Array<value_type>(d));
96 for (ordinal_type j=0; j<d; j++)
97 (*deriv_coeffs)[j] = value_type(1.0);
98 }
99}
100
101template <typename ordinal_type, typename value_type>
106
107template <typename ordinal_type, typename value_type>
108ordinal_type
114
115template <typename ordinal_type, typename value_type>
116ordinal_type
122
123template <typename ordinal_type, typename value_type>
124ordinal_type
130
131template <typename ordinal_type, typename value_type>
132const Teuchos::Array<value_type>&
138
139template <typename ordinal_type, typename value_type>
140const value_type&
142norm_squared(ordinal_type i) const
143{
144 return norms[i];
145}
146
147template <typename ordinal_type, typename value_type>
148Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
151{
152#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
153 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total Triple-Product Tensor Fill Time");
154#endif
155 if (use_old_cijk_alg)
156 return computeTripleProductTensorOld(sz);
157 else
158 return computeTripleProductTensorNew(sz);
159}
160
161template <typename ordinal_type, typename value_type>
162Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
165{
166#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
167 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total Triple-Product Tensor Fill Time");
168#endif
169 if (use_old_cijk_alg)
170 return computeTripleProductTensorOld(d+1);
171 else
172 return computeTripleProductTensorNew(d+1);
173}
174
175template <typename ordinal_type, typename value_type>
176Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
178computeTripleProductTensorOld(ordinal_type order) const
179{
180 // Compute Cijk = < \Psi_i \Psi_j \Psi_k >
181 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
182 Teuchos::rcp(new Sparse3Tensor<ordinal_type, value_type>);
183
184 // Create 1-D triple products
185 Teuchos::Array< Teuchos::RCP<Dense3Tensor<ordinal_type,value_type> > > Cijk_1d(d);
186 for (ordinal_type i=0; i<d; i++)
187 Cijk_1d[i] = bases[i]->computeTripleProductTensor();
188
189 for (ordinal_type j=0; j<sz; j++) {
190 for (ordinal_type i=0; i<sz; i++) {
191 for (ordinal_type k=0; k<order; k++) {
192 value_type c = value_type(1.0);
193 for (ordinal_type l=0; l<d; l++)
194 c *= (*Cijk_1d[l])(terms[i][l],terms[j][l],terms[k][l]);
195 if (std::abs(c/norms[i]) > sparse_tol)
196 Cijk->add_term(i,j,k,c);
197 }
198 }
199 }
200
201 Cijk->fillComplete();
202
203 return Cijk;
204}
205
206template <typename ordinal_type, typename value_type>
207Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
209computeTripleProductTensorNew(ordinal_type order) const
210{
211 // The algorithm for computing Cijk = < \Psi_i \Psi_j \Psi_k > here works
212 // by factoring
213 // < \Psi_i \Psi_j \Psi_k > =
214 // < \psi^1_{i_1}\psi^1_{j_1}\psi^1_{k_1} >_1 x ... x
215 // < \psi^d_{i_d}\psi^d_{j_d}\psi^d_{k_d} >_d
216 // We compute the sparse triple product < \psi^l_i\psi^l_j\psi^l_k >_l
217 // for each dimension l, and then compute all non-zero products of these
218 // terms. The complexity arises from iterating through all possible
219 // combinations, throwing out terms that aren't in the basis and are beyond
220 // the k-order limit provided
221 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
222 Teuchos::rcp(new Sparse3Tensor<ordinal_type, value_type>);
223
224 // Map the specified order limit to a limit on each dimension
225 // Subtract 1 to get the term for the last order we want to include,
226 // add up the orders for each term to get the total order, then add 1
227 MultiIndex<ordinal_type> term = this->term(order-1);
228 ordinal_type k_lim = 0;
229 for (ordinal_type i=0; i<d; i++)
230 k_lim = k_lim + term[i];
231 k_lim++;
232
233 // Create 1-D triple products
234 Teuchos::Array< Teuchos::RCP<Sparse3Tensor<ordinal_type,value_type> > > Cijk_1d(d);
235 for (ordinal_type i=0; i<d; i++) {
236 if (k_lim <= basis_orders[i]+1)
237 Cijk_1d[i] = bases[i]->computeSparseTripleProductTensor(k_lim);
238 else
239 Cijk_1d[i] = bases[i]->computeSparseTripleProductTensor(basis_orders[i]+1);
240 }
241
242 // Create i, j, k iterators for each dimension
243 // Note: we have to supply an initializer in the arrays of iterators to
244 // avoid checked-stl errors about singular iterators
245 typedef Sparse3Tensor<ordinal_type,value_type> Cijk_type;
246 typedef typename Cijk_type::k_iterator k_iterator;
247 typedef typename Cijk_type::kj_iterator kj_iterator;
248 typedef typename Cijk_type::kji_iterator kji_iterator;
249 Teuchos::Array<k_iterator> k_iterators(d, Cijk_1d[0]->k_begin());
250 Teuchos::Array<kj_iterator > j_iterators(d, Cijk_1d[0]->j_begin(k_iterators[0]));
251 Teuchos::Array<kji_iterator > i_iterators(d, Cijk_1d[0]->i_begin(j_iterators[0]));
252 MultiIndex<ordinal_type> terms_i(d), terms_j(d), terms_k(d);
253 ordinal_type sum_i = 0;
254 ordinal_type sum_j = 0;
255 ordinal_type sum_k = 0;
256 for (ordinal_type dim=0; dim<d; dim++) {
257 k_iterators[dim] = Cijk_1d[dim]->k_begin();
258 j_iterators[dim] = Cijk_1d[dim]->j_begin(k_iterators[dim]);
259 i_iterators[dim] = Cijk_1d[dim]->i_begin(j_iterators[dim]);
260 terms_i[dim] = Stokhos::index(i_iterators[dim]);
261 terms_j[dim] = Stokhos::index(j_iterators[dim]);
262 terms_k[dim] = Stokhos::index(k_iterators[dim]);
263 sum_i += terms_i[dim];
264 sum_j += terms_j[dim];
265 sum_k += terms_k[dim];
266 }
267
268 ordinal_type I = 0;
269 ordinal_type J = 0;
270 ordinal_type K = 0;
271 bool inc_i = true;
272 bool inc_j = true;
273 bool inc_k = true;
274 bool stop = false;
275 ordinal_type cnt = 0;
276 while (!stop) {
277
278 // Add term if it is in the basis
279 if (sum_i <= p && sum_j <= p && sum_k <= p) {
280 if (inc_k)
281 K = CPBUtils::compute_index(terms_k, terms, num_terms, p);
282 if (K < order) {
283 if (inc_i)
284 I = CPBUtils::compute_index(terms_i, terms, num_terms, p);
285 if (inc_j)
286 J = CPBUtils::compute_index(terms_j, terms, num_terms, p);
287 value_type c = value_type(1.0);
288 for (ordinal_type dim=0; dim<d; dim++)
289 c *= value(i_iterators[dim]);
290 if (std::abs(c/norms[I]) > sparse_tol)
291 Cijk->add_term(I,J,K,c);
292 }
293 }
294
295 // Increment iterators to the next valid term
296 ordinal_type cdim = 0;
297 bool inc = true;
298 inc_i = false;
299 inc_j = false;
300 inc_k = false;
301 while (inc && cdim < d) {
302 ordinal_type cur_dim = cdim;
303 ++i_iterators[cdim];
304 inc_i = true;
305 if (i_iterators[cdim] != Cijk_1d[cdim]->i_end(j_iterators[cdim])) {
306 terms_i[cdim] = Stokhos::index(i_iterators[cdim]);
307 sum_i = 0;
308 for (int dim=0; dim<d; dim++)
309 sum_i += terms_i[dim];
310 }
311 if (i_iterators[cdim] == Cijk_1d[cdim]->i_end(j_iterators[cdim]) ||
312 sum_i > p) {
313 ++j_iterators[cdim];
314 inc_j = true;
315 if (j_iterators[cdim] != Cijk_1d[cdim]->j_end(k_iterators[cdim])) {
316 terms_j[cdim] = Stokhos::index(j_iterators[cdim]);
317 sum_j = 0;
318 for (int dim=0; dim<d; dim++)
319 sum_j += terms_j[dim];
320 }
321 if (j_iterators[cdim] == Cijk_1d[cdim]->j_end(k_iterators[cdim]) ||
322 sum_j > p) {
323 ++k_iterators[cdim];
324 inc_k = true;
325 if (k_iterators[cdim] != Cijk_1d[cdim]->k_end()) {
326 terms_k[cdim] = Stokhos::index(k_iterators[cdim]);
327 sum_k = 0;
328 for (int dim=0; dim<d; dim++)
329 sum_k += terms_k[dim];
330 }
331 if (k_iterators[cdim] == Cijk_1d[cdim]->k_end() || sum_k > p) {
332 k_iterators[cdim] = Cijk_1d[cdim]->k_begin();
333 ++cdim;
334 terms_k[cur_dim] = Stokhos::index(k_iterators[cur_dim]);
335 sum_k = 0;
336 for (int dim=0; dim<d; dim++)
337 sum_k += terms_k[dim];
338 }
339 else
340 inc = false;
341 j_iterators[cur_dim] =
342 Cijk_1d[cur_dim]->j_begin(k_iterators[cur_dim]);
343 terms_j[cur_dim] = Stokhos::index(j_iterators[cur_dim]);
344 sum_j = 0;
345 for (int dim=0; dim<d; dim++)
346 sum_j += terms_j[dim];
347 }
348 else
349 inc = false;
350 i_iterators[cur_dim] = Cijk_1d[cur_dim]->i_begin(j_iterators[cur_dim]);
351 terms_i[cur_dim] = Stokhos::index(i_iterators[cur_dim]);
352 sum_i = 0;
353 for (int dim=0; dim<d; dim++)
354 sum_i += terms_i[dim];
355 }
356 else
357 inc = false;
358
359 if (sum_i > p || sum_j > p || sum_k > p)
360 inc = true;
361 }
362
363 if (cdim == d)
364 stop = true;
365
366 cnt++;
367 }
368
369 Cijk->fillComplete();
370
371 return Cijk;
372}
373
374template <typename ordinal_type, typename value_type>
375Teuchos::RCP< Stokhos::Dense3Tensor<ordinal_type, value_type> >
378 const Teuchos::RCP< const Teuchos::SerialDenseMatrix<ordinal_type, value_type> >& Bij,
379 const Teuchos::RCP< const Stokhos::Sparse3Tensor<ordinal_type, value_type> >& Cijk) const
380{
381 // Compute Dijk = < \Psi_i \Psi_j \Psi_k' >
382 Teuchos::RCP< Stokhos::Dense3Tensor<ordinal_type, value_type> > Dijk =
383 Teuchos::rcp(new Dense3Tensor<ordinal_type, value_type>(sz));
384 for (ordinal_type i=0; i<sz; i++)
385 for (ordinal_type j=0; j<sz; j++)
386 for (ordinal_type k=0; k<sz; k++)
387 (*Dijk)(i,j,k) = value_type(0.0);
388
389 ordinal_type i,j,m;
390 value_type c;
391 for (ordinal_type k=0; k<sz; k++) {
392 for (typename Cijk_type::k_iterator m_it=Cijk->k_begin();
393 m_it!=Cijk->k_end(); ++m_it) {
394 m = Stokhos::index(m_it);
395 for (typename Cijk_type::kj_iterator j_it = Cijk->j_begin(m_it);
396 j_it != Cijk->j_end(m_it); ++j_it) {
397 j = Stokhos::index(j_it);
398 for (typename Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
399 i_it != Cijk->i_end(j_it); ++i_it) {
400 i = Stokhos::index(i_it);
401 c = Stokhos::value(i_it);
402 (*Dijk)(i,j,k) += (*Bij)(m,k)*c/norms[m];
403 }
404 }
405 }
406 }
407
408 return Dijk;
409}
410
411template <typename ordinal_type, typename value_type>
412Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type, value_type> >
415{
416 // Compute Bij = < \Psi_i \Psi_j' >
417 Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type, value_type> > Bij =
418 Teuchos::rcp(new Teuchos::SerialDenseMatrix<ordinal_type,value_type>(sz,
419 sz));
420
421 // Create products
422 Teuchos::Array< Teuchos::RCP<const Teuchos::SerialDenseMatrix<ordinal_type,value_type> > > Bij_1d(d);
423 for (ordinal_type i=0; i<d; i++)
424 Bij_1d[i] = bases[i]->computeDerivDoubleProductTensor();
425
426 for (ordinal_type i=0; i<sz; i++) {
427 for (ordinal_type k=0; k<sz; k++) {
428 value_type t = value_type(1.0);
429 value_type c = value_type(0.0);
430 for (ordinal_type j=0; j<d; j++) {
431 bool is_zero = false;
432 for (ordinal_type l=0; l<d; l++) {
433 if (l != j && terms[i][l] != terms[k][l])
434 is_zero = true;
435 if (l != j)
436 t *= bases[l]->norm_squared(terms[k][l]);
437 }
438 if (!is_zero)
439 c += t*(*deriv_coeffs)[j]*(*Bij_1d[j])(terms[k][j],terms[i][j]);
440 }
441 (*Bij)(i,k) = c;
442 }
443 }
444
445 return Bij;
446}
447
448template <typename ordinal_type, typename value_type>
449value_type
451evaluateZero(ordinal_type i) const
452{
453 // z = psi_{i_1}(0) * ... * psi_{i_d}(0) where i_1,...,i_d are the basis
454 // terms for coefficient i
455 value_type z = value_type(1.0);
456 for (ordinal_type j=0; j<d; j++)
457 z = z * bases[j]->evaluate(value_type(0.0), terms[i][j]);
458
459 return z;
460}
461
462template <typename ordinal_type, typename value_type>
463void
465evaluateBases(const Teuchos::ArrayView<const value_type>& point,
466 Teuchos::Array<value_type>& basis_vals) const
467{
468 for (ordinal_type j=0; j<d; j++)
469 bases[j]->evaluateBases(point[j], basis_eval_tmp[j]);
470
471 // Only evaluate basis upto number of terms included in basis_pts
472 for (ordinal_type i=0; i<sz; i++) {
473 value_type t = value_type(1.0);
474 for (ordinal_type j=0; j<d; j++)
475 t *= basis_eval_tmp[j][terms[i][j]];
476 basis_vals[i] = t;
477 }
478}
479
480template <typename ordinal_type, typename value_type>
481void
483print(std::ostream& os) const
484{
485 os << "Complete basis of order " << p << ", dimension " << d
486 << ", and size " << sz << ". Component bases:\n";
487 for (ordinal_type i=0; i<d; i++)
488 os << *bases[i];
489 os << "Basis vector norms (squared):\n\t";
490 for (ordinal_type i=0; i<static_cast<ordinal_type>(norms.size()); i++)
491 os << norms[i] << " ";
492 os << "\n";
493}
494
495template <typename ordinal_type, typename value_type>
498term(ordinal_type i) const
499{
500 return terms[i];
501}
502
503template <typename ordinal_type, typename value_type>
504ordinal_type
507{
508 return CPBUtils::compute_index(term, terms, num_terms, p);
509}
510
511template <typename ordinal_type, typename value_type>
512const std::string&
518
519template <typename ordinal_type, typename value_type>
520Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type, value_type> > >
526
527template <typename ordinal_type, typename value_type>
530getMaxOrders() const
531{
532 MultiIndex<ordinal_type> max_orders(d);
533 for (ordinal_type i=0; i<d; ++i)
534 max_orders[i] = basis_orders[i];
535 return max_orders;
536}
static ordinal_type compute_terms(ordinal_type p, ordinal_type d, ordinal_type &sz, Teuchos::Array< MultiIndex< ordinal_type > > &terms, Teuchos::Array< ordinal_type > &num_terms)
Compute the 2-D array of basis terms which maps a basis index into the orders for each basis dimensio...
virtual value_type evaluateZero(ordinal_type i) const
Evaluate basis polynomial i at zero.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
virtual const MultiIndex< ordinal_type > & term(ordinal_type i) const
Get orders of each coordinate polynomial given an index i.
virtual ordinal_type index(const MultiIndex< ordinal_type > &term) const
Get index of the multivariate polynomial given orders of each coordinate.
virtual void print(std::ostream &os) const
Print basis to stream os.
Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > getCoordinateBases() const
Return coordinate bases.
Teuchos::Array< ordinal_type > num_terms
Number of terms up to each order.
virtual MultiIndex< ordinal_type > getMaxOrders() const
Return maximum order allowable for each coordinate basis.
virtual void evaluateBases(const Teuchos::ArrayView< const value_type > &point, Teuchos::Array< value_type > &basis_vals) const
Evaluate basis polynomials at given point point.
Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > bases
Array of bases.
Teuchos::Array< MultiIndex< ordinal_type > > terms
2-D array of basis terms
CompletePolynomialBasis(const Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > &bases, const value_type &sparse_tol=1.0e-12, bool use_old_cijk_alg=false, const Teuchos::RCP< Teuchos::Array< value_type > > &deriv_coeffs=Teuchos::null)
Constructor.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeLinearTripleProductTensor() const
Compute linear triple product tensor where k = 0,1,..,d.
virtual ordinal_type size() const
Return total size of basis.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorNew(ordinal_type order) const
Compute triple product tensor using new algorithm.
virtual Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > computeDerivDoubleProductTensor() const
Compute double product tensor where represents the derivative of in the direction .
ordinal_type order() const
Return order of basis.
Teuchos::RCP< Teuchos::Array< value_type > > deriv_coeffs
Coefficients for derivative.
virtual const std::string & getName() const
Return string name of basis.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorOld(ordinal_type order) const
Compute triple product tensor using old algorithm.
virtual const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
virtual Teuchos::RCP< Stokhos::Dense3Tensor< ordinal_type, value_type > > computeDerivTripleProductTensor(const Teuchos::RCP< const Teuchos::SerialDenseMatrix< ordinal_type, value_type > > &Bij, const Teuchos::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > &Cijk) const
Compute triple product tensor where represents the derivative of in the direction .
ordinal_type dimension() const
Return dimension of basis.
Teuchos::Array< Teuchos::Array< value_type > > basis_eval_tmp
Temporary array used in basis evaluation.
Teuchos::Array< ordinal_type > basis_orders
Array storing order of each basis.
ordinal_type d
Total dimension of basis.
A multidimensional index.
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.
Bi-directional iterator for traversing a sparse array.