Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Stokhos_DerivOrthogPolyExpansionImp.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 #include "Teuchos_Assert.hpp"
46 
47 template <typename ordinal_type, typename value_type>
54  basis(basis_),
55  Bij(Bij_),
56  Cijk(Cijk_),
57  Dijk(Dijk_),
58  sz(basis->size()),
59  A(2*sz,2*sz),
60  B(2*sz,2),
61  piv(2*sz),
62  lapack()
63 {
64 }
65 
66 extern "C" {
67  double dlange_(char*, int*, int*, double*, int*, double*);
68 }
69 
70 template <typename ordinal_type, typename value_type>
74 {
75  if (s == 0 || nrhs == 0)
76  return 0;
77 
78  ordinal_type info;
79 // lapack.GESV(s, nrhs, A.values(), A.numRows(), &(piv[0]), b.values(),
80 // b.numRows(), &info);
81  lapack.GETRF(s, s, A.values(), A.numRows(), &(piv[0]), &info);
82  value_type norm, rcond;
85  norm = 1.0;
86  ordinal_type n = A.numRows();
87  char t = '1';
88  norm = dlange_(&t, &s, &s, A.values(), &n, &work[0]);
89  lapack.GECON('1', s, A.values(), A.numRows(), norm, &rcond, &work[0],
90  &iwork[0], &info);
91  //std::cout << "condition number = " << 1.0/rcond << std::endl;
92  lapack.GETRS('N', s, nrhs, A.values(), A.numRows(), &(piv[0]), B.values(),
93  B.numRows(), &info);
94  return info;
95 }
96 
97 template <typename ordinal_type, typename value_type>
98 void
103 {
104  ordinal_type pc = a.size();
105  if (c.size() != pc)
106  c.resize(pc);
107 
108  value_type* cc = c.coeff();
109  const value_type* ca = a.coeff();
110 
111  for (ordinal_type i=0; i<pc; i++)
112  cc[i] = -ca[i];
113 }
114 
115 template <typename ordinal_type, typename value_type>
116 void
119  const value_type& val)
120 {
121  c[0] += val;
122 }
123 
124 template <typename ordinal_type, typename value_type>
125 void
128  const value_type& val)
129 {
130  c[0] -= val;
131 }
132 
133 template <typename ordinal_type, typename value_type>
134 void
137  const value_type& val)
138 {
139  ordinal_type pc = c.size();
140  value_type* cc = c.coeff();
141  for (ordinal_type i=0; i<pc; i++)
142  cc[i] *= val;
143 }
144 
145 template <typename ordinal_type, typename value_type>
146 void
149  const value_type& val)
150 {
151  ordinal_type pc = c.size();
152  value_type* cc = c.coeff();
153  for (ordinal_type i=0; i<pc; i++)
154  cc[i] /= val;
155 }
156 
157 template <typename ordinal_type, typename value_type>
158 void
163 {
164  ordinal_type xp = x.size();
165  if (c.size() < xp)
166  c.resize(xp);
167 
168  value_type* cc = c.coeff();
169  const value_type* xc = x.coeff();
170  for (ordinal_type i=0; i<xp; i++)
171  cc[i] += xc[i];
172 }
173 
174 template <typename ordinal_type, typename value_type>
175 void
180 {
181  ordinal_type xp = x.size();
182  if (c.size() < xp)
183  c.resize(xp);
184 
185  value_type* cc = c.coeff();
186  const value_type* xc = x.coeff();
187  for (ordinal_type i=0; i<xp; i++)
188  cc[i] -= xc[i];
189 }
190 
191 template <typename ordinal_type, typename value_type>
192 void
197 {
198  ordinal_type p = c.size();
199  ordinal_type xp = x.size();
200  ordinal_type pc;
201  if (p > 1 && xp > 1)
202  pc = sz;
203  else
204  pc = p*xp;
205  TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
206  "Stokhos::DerivOrthogPolyExpansion::timesEqual()" <<
207  ": Expansion size (" << sz <<
208  ") is too small for computation.");
209  if (c.size() != pc)
210  c.resize(pc);
211 
212  value_type* cc = c.coeff();
213  const value_type* xc = x.coeff();
214 
215  if (p > 1 && xp > 1) {
216  // Copy c coefficients into temporary array
218  value_type tmp, cijk;
219  ordinal_type i,j;
220  for (ordinal_type k=0; k<pc; k++) {
221  tmp = value_type(0.0);
222  ordinal_type n = Cijk->num_values(k);
223  for (ordinal_type l=0; l<n; l++) {
224  Cijk->value(k,l,i,j,cijk);
225  if (i < p && j < xp)
226  tmp += cijk*tc[i]*xc[j];
227  }
228  cc[k] = tmp / basis->norm_squared(k);
229  }
230  }
231  else if (p > 1) {
232  for (ordinal_type i=0; i<p; i++)
233  cc[i] *= xc[0];
234  }
235  else if (xp > 1) {
236  for (ordinal_type i=1; i<xp; i++)
237  cc[i] = cc[0]*xc[i];
238  cc[0] *= xc[0];
239  }
240  else {
241  cc[0] *= xc[0];
242  }
243 }
244 
245 template <typename ordinal_type, typename value_type>
246 void
250 {
251  const char* func = "Stokhos::DerivOrthogPolyExpansion::divide()";
252  ordinal_type p = c.size();
253  ordinal_type xp = x.size();
254  ordinal_type pc;
255  if (xp > 1)
256  pc = sz;
257  else
258  pc = p;
259  TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
260  "Stokhos::DerivOrthogPolyExpansion::divideEqual()" <<
261  ": Expansion size (" << sz <<
262  ") is too small for computation.");
263  if (c.size() != pc)
264  c.resize(pc);
265 
266  value_type* cc = c.coeff();
267  const value_type* xc = x.coeff();
268 
269  if (xp > 1) {
270 
271  // Fill A
273  ordinal_type i,j;
274  for (ordinal_type i=0; i<pc; i++)
275  for (ordinal_type j=0; j<pc; j++)
276  A(i,j) = 0.0;
277  for (ordinal_type k=0; k<xp; k++) {
278  ordinal_type n = Cijk->num_values(k);
279  for (ordinal_type l=0; l<n; l++) {
280  Cijk->value(k,l,i,j,cijk);
281  A(i,j) += cijk*xc[k]/basis->norm_squared(i);
282  }
283  }
284 
285  // Fill B
286  for (ordinal_type i=0; i<p; i++)
287  B(i,0) = cc[i];
288  for (ordinal_type i=p; i<pc; i++)
289  B(i,0) = value_type(0.0);
290 
291  // Solve system
292  int info = solve(pc, 1);
293 
294  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
295  func << ": Argument " << info
296  << " for solve had illegal value");
297  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
298  func << ": Diagonal entry " << info
299  << " in LU factorization is exactly zero");
300 
301  // Get coefficients
302  for (ordinal_type i=0; i<pc; i++)
303  cc[i] = B(i,0);
304  }
305  else {
306  for (ordinal_type i=0; i<p; i++)
307  cc[i] /= xc[0];
308  }
309 }
310 
311 template <typename ordinal_type, typename value_type>
312 void
317 {
318  ordinal_type pa = a.size();
319  ordinal_type pb = b.size();
320  ordinal_type pc = pa > pb ? pa : pb;
321  if (c.size() != pc)
322  c.resize(pc);
323 
324  const value_type* ca = a.coeff();
325  const value_type* cb = b.coeff();
326  value_type* cc = c.coeff();
327 
328  if (pa > pb) {
329  for (ordinal_type i=0; i<pb; i++)
330  cc[i] = ca[i] + cb[i];
331  for (ordinal_type i=pb; i<pc; i++)
332  cc[i] = ca[i];
333  }
334  else {
335  for (ordinal_type i=0; i<pa; i++)
336  cc[i] = ca[i] + cb[i];
337  for (ordinal_type i=pa; i<pc; i++)
338  cc[i] = cb[i];
339  }
340 }
341 
342 template <typename ordinal_type, typename value_type>
343 void
346  const value_type& a,
348 {
349  ordinal_type pc = b.size();
350  if (c.size() != pc)
351  c.resize(pc);
352 
353  const value_type* cb = b.coeff();
354  value_type* cc = c.coeff();
355 
356  cc[0] = a + cb[0];
357  for (ordinal_type i=1; i<pc; i++)
358  cc[i] = cb[i];
359 }
360 
361 template <typename ordinal_type, typename value_type>
362 void
366  const value_type& b)
367 {
368  ordinal_type pc = a.size();
369  if (c.size() != pc)
370  c.resize(pc);
371 
372  const value_type* ca = a.coeff();
373  value_type* cc = c.coeff();
374 
375  cc[0] = ca[0] + b;
376  for (ordinal_type i=1; i<pc; i++)
377  cc[i] = ca[i];
378 }
379 
380 template <typename ordinal_type, typename value_type>
381 void
386 {
387  ordinal_type pa = a.size();
388  ordinal_type pb = b.size();
389  ordinal_type pc = pa > pb ? pa : pb;
390  if (c.size() != pc)
391  c.resize(pc);
392 
393  const value_type* ca = a.coeff();
394  const value_type* cb = b.coeff();
395  value_type* cc = c.coeff();
396 
397  if (pa > pb) {
398  for (ordinal_type i=0; i<pb; i++)
399  cc[i] = ca[i] - cb[i];
400  for (ordinal_type i=pb; i<pc; i++)
401  cc[i] = ca[i];
402  }
403  else {
404  for (ordinal_type i=0; i<pa; i++)
405  cc[i] = ca[i] - cb[i];
406  for (ordinal_type i=pa; i<pc; i++)
407  cc[i] = -cb[i];
408  }
409 }
410 
411 template <typename ordinal_type, typename value_type>
412 void
415  const value_type& a,
417 {
418  ordinal_type pc = b.size();
419  if (c.size() != pc)
420  c.resize(pc);
421 
422  const value_type* cb = b.coeff();
423  value_type* cc = c.coeff();
424 
425  cc[0] = a - cb[0];
426  for (ordinal_type i=1; i<pc; i++)
427  cc[i] = -cb[i];
428 }
429 
430 template <typename ordinal_type, typename value_type>
431 void
435  const value_type& b)
436 {
437  ordinal_type pc = a.size();
438  if (c.size() != pc)
439  c.resize(pc);
440 
441  const value_type* ca = a.coeff();
442  value_type* cc = c.coeff();
443 
444  cc[0] = ca[0] - b;
445  for (ordinal_type i=1; i<pc; i++)
446  cc[i] = ca[i];
447 }
448 
449 template <typename ordinal_type, typename value_type>
450 void
455 {
456  ordinal_type pa = a.size();
457  ordinal_type pb = b.size();
458  ordinal_type pc;
459  if (pa > 1 && pb > 1)
460  pc = sz;
461  else
462  pc = pa*pb;
463  TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
464  "Stokhos::DerivOrthogPolyExpansion::times()" <<
465  ": Expansion size (" << sz <<
466  ") is too small for computation.");
467  if (c.size() != pc)
468  c.resize(pc);
469 
470  const value_type* ca = a.coeff();
471  const value_type* cb = b.coeff();
472  value_type* cc = c.coeff();
473 
474  if (pa > 1 && pb > 1) {
475  value_type tmp, cijk;
476  ordinal_type i,j;
477  for (ordinal_type k=0; k<pc; k++) {
478  tmp = value_type(0.0);
479  ordinal_type n = Cijk->num_values(k);
480  for (ordinal_type l=0; l<n; l++) {
481  Cijk->value(k,l,i,j,cijk);
482  if (i < pa && j < pb)
483  tmp += cijk*ca[i]*cb[j];
484  }
485  cc[k] = tmp / basis->norm_squared(k);
486  }
487  }
488  else if (pa > 1) {
489  for (ordinal_type i=0; i<pc; i++)
490  cc[i] = ca[i]*cb[0];
491  }
492  else if (pb > 1) {
493  for (ordinal_type i=0; i<pc; i++)
494  cc[i] = ca[0]*cb[i];
495  }
496  else {
497  cc[0] = ca[0]*cb[0];
498  }
499 }
500 
501 template <typename ordinal_type, typename value_type>
502 void
505  const value_type& a,
507 {
508  ordinal_type pc = b.size();
509  if (c.size() != pc)
510  c.resize(pc);
511 
512  const value_type* cb = b.coeff();
513  value_type* cc = c.coeff();
514 
515  for (ordinal_type i=0; i<pc; i++)
516  cc[i] = a*cb[i];
517 }
518 
519 template <typename ordinal_type, typename value_type>
520 void
524  const value_type& b)
525 {
526  ordinal_type pc = a.size();
527  if (c.size() != pc)
528  c.resize(pc);
529 
530  const value_type* ca = a.coeff();
531  value_type* cc = c.coeff();
532 
533  for (ordinal_type i=0; i<pc; i++)
534  cc[i] = ca[i]*b;
535 }
536 
537 template <typename ordinal_type, typename value_type>
538 void
543 {
544  const char* func = "Stokhos::DerivOrthogPolyExpansion::divide()";
545  ordinal_type pa = a.size();
546  ordinal_type pb = b.size();
547  ordinal_type pc;
548  if (pb > 1)
549  pc = sz;
550  else
551  pc = pa;
552  TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
553  "Stokhos::DerivOrthogPolyExpansion::divide()" <<
554  ": Expansion size (" << sz <<
555  ") is too small for computation.");
556  if (c.size() != pc)
557  c.resize(pc);
558 
559  const value_type* ca = a.coeff();
560  const value_type* cb = b.coeff();
561  value_type* cc = c.coeff();
562 
563  if (pb > 1) {
564 
565  // Fill A
567  ordinal_type i,j;
568  for (ordinal_type i=0; i<pc; i++)
569  for (ordinal_type j=0; j<pc; j++)
570  A(i,j) = 0.0;
571 
572  for (ordinal_type k=0; k<pb; k++) {
573  ordinal_type n = Cijk->num_values(k);
574  for (ordinal_type l=0; l<n; l++) {
575  Cijk->value(k,l,i,j,cijk);
576  A(i,j) += cijk*cb[k] / basis->norm_squared(i);
577  }
578  }
579 
580  // Fill B
581  for (ordinal_type i=0; i<pa; i++)
582  B(i,0) = ca[i];
583  for (ordinal_type i=pa; i<pc; i++)
584  B(i,0) = value_type(0.0);
585 
586  // Solve system
587  int info = solve(pc, 1);
588 
589  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
590  func << ": Argument " << info
591  << " for solve had illegal value");
592  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
593  func << ": Diagonal entry " << info
594  << " in LU factorization is exactly zero");
595 
596  // Get coefficients
597  for (ordinal_type i=0; i<pc; i++)
598  cc[i] = B(i,0);
599  }
600  else {
601  for (ordinal_type i=0; i<pa; i++)
602  cc[i] = ca[i]/cb[0];
603  }
604 }
605 
606 template <typename ordinal_type, typename value_type>
607 void
610  const value_type& a,
612 {
613  const char* func = "Stokhos::DerivOrthogPolyExpansion::divide()";
614  ordinal_type pb = b.size();
615  ordinal_type pc;
616  if (pb > 1)
617  pc = sz;
618  else
619  pc = 1;
620  if (c.size() != pc)
621  c.resize(pc);
622 
623  const value_type* cb = b.coeff();
624  value_type* cc = c.coeff();
625 
626  if (pb > 1) {
627 
628  // Fill A
630  ordinal_type i,j;
631  for (ordinal_type i=0; i<pc; i++)
632  for (ordinal_type j=0; j<pc; j++)
633  A(i,j) = 0.0;
634 
635  for (ordinal_type k=0; k<pb; k++) {
636  ordinal_type n = Cijk->num_values(k);
637  for (ordinal_type l=0; l<n; l++) {
638  Cijk->value(k,l,i,j,cijk);
639  A(i,j) += cijk*cb[k] / basis->norm_squared(i);
640  }
641  }
642 
643  // Fill B
644  B(0,0) = a;
645  for (ordinal_type i=1; i<pc; i++)
646  B(i,0) = value_type(0.0);
647 
648  // Solve system
649  int info = solve(pc, 1);
650 
651  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
652  func << ": Argument " << info
653  << " for solve had illegal value");
654  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
655  func << ": Diagonal entry " << info
656  << " in LU factorization is exactly zero");
657 
658  // Get coefficients
659  for (ordinal_type i=0; i<pc; i++)
660  cc[i] = B(i,0);
661  }
662  else
663  cc[0] = a / cb[0];
664 }
665 
666 template <typename ordinal_type, typename value_type>
667 void
671  const value_type& b)
672 {
673  ordinal_type pc = a.size();
674  if (c.size() != pc)
675  c.resize(pc);
676 
677  const value_type* ca = a.coeff();
678  value_type* cc = c.coeff();
679 
680  for (ordinal_type i=0; i<pc; i++)
681  cc[i] = ca[i]/b;
682 }
683 
684 template <typename ordinal_type, typename value_type>
685 void
689 {
690  const char* func = "Stokhos::DerivOrthogPolyExpansion::exp()";
691  ordinal_type pa = a.size();
692  ordinal_type pc;
693  if (pa > 1)
694  pc = sz;
695  else
696  pc = 1;
697  if (c.size() != pc)
698  c.resize(pc);
699 
700  const value_type* ca = a.coeff();
701  value_type* cc = c.coeff();
702 
703  if (pa > 1) {
704  value_type psi_0 = basis->evaluateZero(0);
705 
706  // Fill A and B
707  for (ordinal_type i=1; i<pc; i++) {
708  B(i-1,0) = 0.0;
709  for (ordinal_type j=1; j<pc; j++) {
710  A(i-1,j-1) = (*Bij)(i-1,j);
711  for (ordinal_type k=1; k<pa; k++)
712  A(i-1,j-1) -= ca[k]*(*Dijk)(i-1,j,k);
713  B(i-1,0) += ca[j]*(*Bij)(i-1,j);
714  }
715  B(i-1,0) *= psi_0;
716  }
717 
718  // Solve system
719  int info = solve(pc-1, 1);
720 
721  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
722  func << ": Argument " << info
723  << " for solve had illegal value");
724  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
725  func << ": Diagonal entry " << info
726  << " in LU factorization is exactly zero");
727 
728  // Compute order-0 coefficient
729  value_type s = psi_0 * ca[0];
730  value_type t = psi_0;
731  for (ordinal_type i=1; i<pc; i++) {
732  s += basis->evaluateZero(i) * ca[i];
733  t += basis->evaluateZero(i) * B(i-1,0);
734  }
735  s = std::exp(s);
736  cc[0] = (s/t);
737 
738  // Compute remaining coefficients
739  for (ordinal_type i=1; i<pc; i++)
740  cc[i] = B(i-1,0) * cc[0];
741  }
742  else
743  cc[0] = std::exp(ca[0]);
744 }
745 
746 template <typename ordinal_type, typename value_type>
747 void
751 {
752  const char* func = "Stokhos::DerivOrthogPolyExpansion::log()";
753  ordinal_type pa = a.size();
754  ordinal_type pc;
755  if (pa > 1)
756  pc = sz;
757  else
758  pc = 1;
759  if (c.size() != pc)
760  c.resize(pc);
761 
762  const value_type* ca = a.coeff();
763  value_type* cc = c.coeff();
764 
765  if (pa > 1) {
766  value_type psi_0 = basis->evaluateZero(0);
767 
768  // Fill A and B
769  for (ordinal_type i=1; i<pc; i++) {
770  B(i-1,0) = 0.0;
771  for (ordinal_type j=1; j<pc; j++) {
772  A(i-1,j-1) = 0.0;
773  for (ordinal_type k=0; k<pa; k++)
774  A(i-1,j-1) += ca[k]*(*Dijk)(i-1,k,j);
775  B(i-1,0) += ca[j]*(*Bij)(i-1,j);
776  }
777  }
778 
779  // Solve system
780  int info = solve(pc-1, 1);
781 
782  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
783  func << ": Argument " << info
784  << " for solve had illegal value");
785  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
786  func << ": Diagonal entry " << info
787  << " in LU factorization is exactly zero");
788 
789  // Compute order-0 coefficient
790  value_type s = psi_0 * ca[0];
791  value_type t = value_type(0.0);
792  for (ordinal_type i=1; i<pc; i++) {
793  s += basis->evaluateZero(i) * ca[i];
794  t += basis->evaluateZero(i) * B(i-1,0);
795  }
796  cc[0] = (std::log(s) - t) / psi_0;
797 
798  // Compute remaining coefficients
799  for (ordinal_type i=1; i<pc; i++)
800  cc[i] = B(i-1,0);
801  }
802  else
803  cc[0] = std::log(ca[0]);
804 }
805 
806 template <typename ordinal_type, typename value_type>
807 void
811 {
812  if (a.size() > 1) {
813  log(c,a);
814  divide(c,c,std::log(10.0));
815  }
816  else {
817  if (c.size() != 1)
818  c.resize(1);
819  c[0] = std::log10(a[0]);
820  }
821 }
822 
823 template <typename ordinal_type, typename value_type>
824 void
828 {
829  if (a.size() > 1) {
830  log(c,a);
831  timesEqual(c,value_type(0.5));
832  exp(c,c);
833  }
834  else {
835  if (c.size() != 1)
836  c.resize(1);
837  c[0] = std::sqrt(a[0]);
838  }
839 }
840 
841 template <typename ordinal_type, typename value_type>
842 void
846 {
847  if (a.size() > 1) {
848  log(c,a);
849  timesEqual(c,value_type(1.0/3.0));
850  exp(c,c);
851  }
852  else {
853  if (c.size() != 1)
854  c.resize(1);
855  c[0] = std::cbrt(a[0]);
856  }
857 }
858 
859 template <typename ordinal_type, typename value_type>
860 void
865 {
866  if (a.size() > 1 || b.size() > 1) {
867  log(c,a);
868  timesEqual(c,b);
869  exp(c,c);
870  }
871  else {
872  if (c.size() != 1)
873  c.resize(1);
874  c[0] = std::pow(a[0], b[0]);
875  }
876 }
877 
878 template <typename ordinal_type, typename value_type>
879 void
882  const value_type& a,
884 {
885  if (b.size() > 1) {
886  times(c,std::log(a),b);
887  exp(c,c);
888  }
889  else {
890  if (c.size() != 1)
891  c.resize(1);
892  c[0] = std::pow(a, b[0]);
893  }
894 }
895 
896 template <typename ordinal_type, typename value_type>
897 void
901  const value_type& b)
902 {
903  if (a.size() > 1) {
904  log(c,a);
905  timesEqual(c,b);
906  exp(c,c);
907  }
908  else {
909  if (c.size() != 1)
910  c.resize(1);
911  c[0] = std::pow(a[0], b);
912  }
913 }
914 
915 template <typename ordinal_type, typename value_type>
916 void
921 {
922  const char* func = "Stokhos::DerivOrthogPolyExpansion::sincos()";
923  ordinal_type pa = a.size();
924  ordinal_type pc;
925  if (pa > 1)
926  pc = sz;
927  else
928  pc = 1;
929  if (c.size() != pc)
930  c.resize(pc);
931  if (s.size() != pc)
932  s.resize(pc);
933 
934  const value_type* ca = a.coeff();
935  value_type* cs = s.coeff();
936  value_type* cc = c.coeff();
937 
938  if (pa > 1) {
939  value_type psi_0 = basis->evaluateZero(0);
940  ordinal_type offset = pc-1;
941 
942  // Fill A and b
943  B.putScalar(value_type(0.0));
944  value_type tmp, tmp2;
945  for (ordinal_type i=1; i<pc; i++) {
946  tmp2 = value_type(0.0);
947  for (ordinal_type j=1; j<pc; j++) {
948  tmp = (*Bij)(i-1,j);
949  A(i-1,j-1) = tmp;
950  A(i-1+offset,j-1+offset) = tmp;
951  tmp = value_type(0.0);
952  for (ordinal_type k=1; k<pa; k++)
953  tmp += ca[k]*(*Dijk)(i-1,j,k);
954  A(i-1+offset,j-1) = tmp;
955  A(i-1,j-1+offset) = -tmp;
956  tmp2 += ca[j]*(*Bij)(i-1,j);
957  }
958  B(i-1,0) = tmp2*psi_0;
959  B(i-1+offset,1) = tmp2*psi_0;
960  }
961 
962  // Solve system
963  int info = solve(2*pc-2, 2);
964 
965  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
966  func << ": Argument " << info
967  << " for solve had illegal value");
968  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
969  func << ": Diagonal entry " << info
970  << " in LU factorization is exactly zero");
971 
972  // Compute order-0 coefficients
973  value_type t = psi_0 * ca[0];
974  value_type a00 = psi_0;
975  value_type a01 = value_type(0.0);
976  value_type a10 = value_type(0.0);
977  value_type a11 = psi_0;
978  value_type b0 = B(0,0);
979  value_type b1 = B(1,0);
980  for (ordinal_type i=1; i<pc; i++) {
981  t += basis->evaluateZero(i) * ca[i];
982  a00 -= basis->evaluateZero(i) * B(i-1,1);
983  a01 += basis->evaluateZero(i) * B(i-1,0);
984  a10 -= basis->evaluateZero(i) * B(i-1+offset,1);
985  a11 += basis->evaluateZero(i) * B(i-1+offset,0);
986  }
987  A(0,0) = a00;
988  A(0,1) = a01;
989  A(1,0) = a10;
990  A(1,1) = a11;
991  B(0,0) = std::sin(t);
992  B(1,0) = std::cos(t);
993 
994  info = solve(2, 1);
995 
996  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
997  func << ": Argument " << info
998  << " for (2x2) solve had illegal value");
999  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
1000  func << ": Diagonal entry " << info
1001  << " in (2x2) LU factorization is exactly zero");
1002  cs[0] = B(0,0);
1003  cc[0] = B(1,0);
1004 
1005  // Compute remaining coefficients
1006  B(0,0) = b0;
1007  B(1,0) = b1;
1008  for (ordinal_type i=1; i<pc; i++) {
1009  cs[i] = cc[0]*B(i-1,0) - cs[0]*B(i-1,1);
1010  cc[i] = cc[0]*B(i-1+offset,0) - cs[0]*B(i-1+offset,1);
1011  }
1012  }
1013  else {
1014  cs[0] = std::sin(ca[0]);
1015  cc[0] = std::cos(ca[0]);
1016  }
1017 }
1018 
1019 template <typename ordinal_type, typename value_type>
1020 void
1024 {
1025  if (a.size() > 1) {
1027  sincos(s, c, a);
1028  }
1029  else {
1030  if (s.size() != 1)
1031  s.resize(1);
1032  s[0] = std::sin(a[0]);
1033  }
1034 }
1035 
1036 template <typename ordinal_type, typename value_type>
1037 void
1041 {
1042  if (a.size() > 1) {
1044  sincos(s, c, a);
1045  }
1046  else {
1047  if (c.size() != 1)
1048  c.resize(1);
1049  c[0] = std::cos(a[0]);
1050  }
1051 }
1052 
1053 template <typename ordinal_type, typename value_type>
1054 void
1058 {
1059  if (a.size() > 1) {
1061  sincos(t, c, a);
1062  divideEqual(t,c);
1063  }
1064  else {
1065  if (t.size() != 1)
1066  t.resize(1);
1067  t[0] = std::tan(a[0]);
1068  }
1069 }
1070 
1071 template <typename ordinal_type, typename value_type>
1072 void
1077 {
1078  const char* func = "Stokhos::DerivOrthogPolyExpansion::sinhcosh()";
1079  ordinal_type pa = a.size();
1080  ordinal_type pc;
1081  if (pa > 1)
1082  pc = sz;
1083  else
1084  pc = 1;
1085  if (c.size() != pc)
1086  c.resize(pc);
1087  if (s.size() != pc)
1088  s.resize(pc);
1089 
1090  const value_type* ca = a.coeff();
1091  value_type* cs = s.coeff();
1092  value_type* cc = c.coeff();
1093 
1094  if (pa > 1) {
1095  value_type psi_0 = basis->evaluateZero(0);
1096  ordinal_type offset = pc-1;
1097 
1098  // Fill A and b
1099  B.putScalar(value_type(0.0));
1100  value_type tmp, tmp2;
1101  for (ordinal_type i=1; i<pc; i++) {
1102  tmp2 = value_type(0.0);
1103  for (ordinal_type j=1; j<pc; j++) {
1104  tmp = (*Bij)(i-1,j);
1105  A(i-1,j-1) = tmp;
1106  A(i-1+offset,j-1+offset) = tmp;
1107  tmp = value_type(0.0);
1108  for (ordinal_type k=1; k<pa; k++)
1109  tmp += ca[k]*(*Dijk)(i-1,j,k);
1110  A(i-1+offset,j-1) = -tmp;
1111  A(i-1,j-1+offset) = -tmp;
1112  tmp2 += ca[j]*(*Bij)(i-1,j);
1113  }
1114  B(i-1,0) = tmp2*psi_0;
1115  B(i-1+offset,1) = tmp2*psi_0;
1116  }
1117 
1118  // Solve system
1119  int info = solve(2*pc-2, 2);
1120 
1121  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
1122  func << ": Argument " << info
1123  << " for solve had illegal value");
1124  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
1125  func << ": Diagonal entry " << info
1126  << " in LU factorization is exactly zero");
1127 
1128  // Compute order-0 coefficients
1129  value_type t = psi_0 * ca[0];
1130  value_type a00 = psi_0;
1131  value_type a01 = value_type(0.0);
1132  value_type a10 = value_type(0.0);
1133  value_type a11 = psi_0;
1134  value_type b0 = B(0,0);
1135  value_type b1 = B(1,0);
1136  for (ordinal_type i=1; i<pc; i++) {
1137  t += basis->evaluateZero(i) * ca[i];
1138  a00 += basis->evaluateZero(i) * B(i-1,1);
1139  a01 += basis->evaluateZero(i) * B(i-1,0);
1140  a10 += basis->evaluateZero(i) * B(i-1+offset,1);
1141  a11 += basis->evaluateZero(i) * B(i-1+offset,0);
1142  }
1143  A(0,0) = a00;
1144  A(0,1) = a01;
1145  A(1,0) = a10;
1146  A(1,1) = a11;
1147  B(0,0) = std::sinh(t);
1148  B(1,0) = std::cosh(t);
1149  info = solve(2, 1);
1150  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
1151  func << ": Argument " << info
1152  << " for (2x2) solve had illegal value");
1153  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
1154  func << ": Diagonal entry " << info
1155  << " in (2x2) LU factorization is exactly zero");
1156  cs[0] = B(0,0);
1157  cc[0] = B(1,0);
1158 
1159  // Compute remaining coefficients
1160  B(0,0) = b0;
1161  B(1,0) = b1;
1162  for (ordinal_type i=1; i<pc; i++) {
1163  cs[i] = cc[0]*B(i-1,0) + cs[0]*B(i-1,1);
1164  cc[i] = cc[0]*B(i-1+offset,0) + cs[0]*B(i-1+offset,1);
1165  }
1166  }
1167  else {
1168  cs[0] = std::sinh(ca[0]);
1169  cc[0] = std::cosh(ca[0]);
1170  }
1171 }
1172 
1173 template <typename ordinal_type, typename value_type>
1174 void
1178 {
1179  if (a.size() > 1) {
1181  sinhcosh(s, c, a);
1182  }
1183  else {
1184  if (s.size() != 1)
1185  s.resize(1);
1186  s[0] = std::sinh(a[0]);
1187  }
1188 }
1189 
1190 template <typename ordinal_type, typename value_type>
1191 void
1195 {
1196  if (a.size() > 1) {
1198  sinhcosh(s, c, a);
1199  }
1200  else {
1201  if (c.size() != 1)
1202  c.resize(1);
1203  c[0] = std::cosh(a[0]);
1204  }
1205 }
1206 
1207 template <typename ordinal_type, typename value_type>
1208 void
1212 {
1213  if (a.size() > 1) {
1215  sinhcosh(t, c, a);
1216  divideEqual(t,c);
1217  }
1218  else {
1219  if (t.size() != 1)
1220  t.resize(1);
1221  t[0] = std::tanh(a[0]);
1222  }
1223 }
1224 
1225 template <typename ordinal_type, typename value_type>
1226 template <typename OpT>
1227 void
1229 quad(const OpT& quad_func,
1233 {
1234  const char* func = "Stokhos::DerivOrthogPolyExpansion::quad()";
1235  ordinal_type pa = a.size();
1236  ordinal_type pb = b.size();
1237  ordinal_type pc = sz;
1238  if (c.size() != pc)
1239  c.resize(pc);
1240 
1241  const value_type* ca = a.coeff();
1242  const value_type* cb = b.coeff();
1243  value_type* cc = c.coeff();
1244 
1245  value_type psi_0 = basis->evaluateZero(0);
1246 
1247  // Fill A and B
1248  for (ordinal_type i=1; i<pc; i++) {
1249  B(i-1,0) = 0.0;
1250  for (ordinal_type j=1; j<pc; j++) {
1251  A(i-1,j-1) = 0.0;
1252  for (ordinal_type k=0; k<pb; k++)
1253  A(i-1,j-1) += cb[k]*(*Dijk)(i-1,k,j);
1254  }
1255  for (ordinal_type j=1; j<pa; j++) {
1256  B(i-1,0) += ca[j]*(*Bij)(i-1,j);
1257  }
1258  }
1259 
1260  // Solve system
1261  int info = solve(pc-1, 1);
1262 
1263  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
1264  func << ": Argument " << info
1265  << " for solve had illegal value");
1266  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
1267  func << ": Diagonal entry " << info
1268  << " in LU factorization is exactly zero");
1269 
1270  // Compute order-0 coefficient
1271  value_type s = psi_0 * ca[0];
1272  value_type t = value_type(0.0);
1273  for (ordinal_type i=1; i<pc; i++) {
1274  s += basis->evaluateZero(i) * ca[i];
1275  t += basis->evaluateZero(i) * B(i-1,0);
1276  }
1277  cc[0] = (quad_func(s) - t) / psi_0;
1278 
1279  // Get coefficients
1280  for (ordinal_type i=1; i<pc; i++)
1281  cc[i] = B(i-1,0);
1282 }
1283 
1284 template <typename ordinal_type, typename value_type>
1285 void
1289 {
1290  if (a.size() > 1) {
1291  times(c,a,a);
1292  minus(c,value_type(1.0),c);
1293  sqrt(c,c);
1294  timesEqual(c,value_type(-1.0));
1295  quad(acos_quad_func(), c, a, c);
1296  }
1297  else {
1298  if (c.size() != 1)
1299  c.resize(1);
1300  c[0] = std::acos(a[0]);
1301  }
1302 }
1303 
1304 template <typename ordinal_type, typename value_type>
1305 void
1309 {
1310  if (a.size() > 1) {
1311  times(c,a,a);
1312  minus(c,value_type(1.0),c);
1313  sqrt(c,c);
1314  quad(asin_quad_func(), c, a, c);
1315  }
1316  else {
1317  if (c.size() != 1)
1318  c.resize(1);
1319  c[0] = std::asin(a[0]);
1320  }
1321 }
1322 
1323 template <typename ordinal_type, typename value_type>
1324 void
1328 {
1329  if (a.size() > 1) {
1330  times(c,a,a);
1331  plusEqual(c,value_type(1.0));
1332  quad(atan_quad_func(), c, a, c);
1333  }
1334  else {
1335  if (c.size() != 1)
1336  c.resize(1);
1337  c[0] = std::atan(a[0]);
1338  }
1339 }
1340 
1341 // template <typename ordinal_type, typename value_type>
1342 // Hermite<ordinal_type, value_type>
1343 // atan2(const Hermite<ordinal_type, value_type>& a,
1344 // const Hermite<ordinal_type, value_type>& b)
1345 // {
1346 // Hermite<ordinal_type, value_type> c = atan(a/b);
1347 // c.fastAccessCoeff(0) = atan2(a.coeff(0),b.coeff(0));
1348 // }
1349 
1350 // template <typename ordinal_type, typename value_type>
1351 // Hermite<ordinal_type, value_type>
1352 // atan2(const T& a,
1353 // const Hermite<ordinal_type, value_type>& b)
1354 // {
1355 // Hermite<ordinal_type, value_type> c = atan(a/b);
1356 // c.fastAccessCoeff(0) = atan2(a,b.coeff(0));
1357 // }
1358 
1359 // template <typename ordinal_type, typename value_type>
1360 // Hermite<ordinal_type, value_type>
1361 // atan2(const Hermite<ordinal_type, value_type>& a,
1362 // const T& b)
1363 // {
1364 // Hermite<ordinal_type, value_type> c = atan(a/b);
1365 // c.fastAccessCoeff(0) = atan2(a.coeff(0),b);
1366 // }
1367 
1368 template <typename ordinal_type, typename value_type>
1369 void
1373 {
1374  if (a.size() > 1) {
1375  times(c,a,a);
1376  minusEqual(c,value_type(1.0));
1377  sqrt(c,c);
1378  quad(acosh_quad_func(), c, a, c);
1379  }
1380  else {
1381  if (c.size() != 1)
1382  c.resize(1);
1383  c[0] = std::log(a[0]+std::sqrt(a[0]*a[0]-value_type(1.0)));
1384  }
1385 }
1386 
1387 template <typename ordinal_type, typename value_type>
1388 void
1392 {
1393  if (a.size() > 1) {
1394  times(c,a,a);
1395  plusEqual(c,value_type(1.0));
1396  sqrt(c,c);
1397  quad(asinh_quad_func(), c, a, c);
1398  }
1399  else {
1400  if (c.size() != 1)
1401  c.resize(1);
1402  c[0] = std::log(a[0]+std::sqrt(a[0]*a[0]+value_type(1.0)));
1403  }
1404 }
1405 
1406 template <typename ordinal_type, typename value_type>
1407 void
1411 {
1412  if (a.size() > 1) {
1413  times(c,a,a);
1414  minus(c,value_type(1.0),c);
1415  quad(atanh_quad_func(), c, a, c);
1416  }
1417  else {
1418  if (c.size() != 1)
1419  c.resize(1);
1420  c[0] = 0.5*std::log((value_type(1.0)+a[0])/(value_type(1.0)-a[0]));
1421  }
1422 }
1423 
1424 template <typename ordinal_type, typename value_type>
1425 void
1429 {
1430  if (a[0] >= 0)
1431  c = a;
1432  else
1433  unaryMinus(c,a);
1434 }
1435 
1436 template <typename ordinal_type, typename value_type>
1437 void
1441 {
1442  if (a[0] >= 0)
1443  c = a;
1444  else
1445  unaryMinus(c,a);
1446 }
1447 
1448 template <typename ordinal_type, typename value_type>
1449 void
1454 {
1455  if (a[0] >= b[0])
1456  c = a;
1457  else
1458  c = b;
1459 }
1460 
1461 template <typename ordinal_type, typename value_type>
1462 void
1465  const value_type& a,
1467 {
1468  if (a >= b[0]) {
1470  c[0] = a;
1471  }
1472  else
1473  c = b;
1474 }
1475 
1476 template <typename ordinal_type, typename value_type>
1477 void
1481  const value_type& b)
1482 {
1483  if (a[0] >= b)
1484  c = a;
1485  else {
1487  c[0] = b;
1488  }
1489 }
1490 
1491 template <typename ordinal_type, typename value_type>
1492 void
1497 {
1498  if (a[0] <= b[0])
1499  c = a;
1500  else
1501  c = b;
1502 }
1503 
1504 template <typename ordinal_type, typename value_type>
1505 void
1508  const value_type& a,
1510 {
1511  if (a <= b[0]) {
1513  c[0] = a;
1514  }
1515  else
1516  c = b;
1517 }
1518 
1519 template <typename ordinal_type, typename value_type>
1520 void
1524  const value_type& b)
1525 {
1526  if (a[0] <= b)
1527  c = a;
1528  else {
1530  c[0] = b;
1531  }
1532 }
1533 
1534 template <typename ordinal_type, typename value_type>
1535 void
1539 {
1540  ordinal_type pc = a.size();
1541 
1542  const value_type* ca = a.coeff();
1543  value_type* cc = c.coeff();
1544 
1545  for (ordinal_type i=0; i<pc; i++) {
1546  cc[i] = 0.0;
1547  for (ordinal_type j=0; j<pc; j++)
1548  cc[i] += ca[j]*(*Bij)(i,j);
1549  cc[i] /= basis->norm_squared(i);
1550  }
1551 }
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
void cbrt(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > tan(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > sinh(const PCE< Storage > &a)
DerivOrthogPolyExpansion(const Teuchos::RCP< const DerivBasis< ordinal_type, value_type > > &basis, const Teuchos::RCP< const Teuchos::SerialDenseMatrix< ordinal_type, value_type > > &Bij, const Teuchos::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > &Cijk, const Teuchos::RCP< const Stokhos::Dense3Tensor< ordinal_type, value_type > > &Dijk)
Constructor.
void acos(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void resize(ordinal_type sz)
Resize coefficient array (coefficients are preserved)
Abstract base class for multivariate orthogonal polynomials that support computing double and triple ...
void minusEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
double dlange_(char *, int *, int *, double *, int *, double *)
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
void quad(const OpT &quad_func, OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
KOKKOS_INLINE_FUNCTION PCE< Storage > pow(const PCE< Storage > &a, const PCE< Storage > &b)
void cos(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void plus(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void timesEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
void asin(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
pointer coeff()
Return coefficient array.
KOKKOS_INLINE_FUNCTION PCE< Storage > tanh(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > cbrt(const PCE< Storage > &a)
void cosh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
static T * get_and_fill(int sz)
Get memory for new array of length sz and fill with zeros.
KOKKOS_INLINE_FUNCTION PCE< Storage > acos(const PCE< Storage > &a)
void unaryMinus(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void acosh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void derivative(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void divideEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
void tan(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void log10(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void divide(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
KOKKOS_INLINE_FUNCTION PCE< Storage > cosh(const PCE< Storage > &a)
void min(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
void atan(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void sinhcosh(OrthogPolyApprox< ordinal_type, value_type, node_type > &s, OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void sqrt(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void minus(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
KOKKOS_INLINE_FUNCTION PCE< Storage > atan(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > exp(const PCE< Storage > &a)
void tanh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
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)
void asinh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void abs(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Class to store coefficients of a projection onto an orthogonal polynomial basis.
expr val()
void sincos(OrthogPolyApprox< ordinal_type, value_type, node_type > &s, OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void fabs(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void times(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
void max(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
void log(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void pow(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
Data structure storing a dense 3-tensor C(i,j,k).
KOKKOS_INLINE_FUNCTION PCE< Storage > sin(const PCE< Storage > &a)
void sin(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
ordinal_type size() const
Return size.
void atanh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > log(const PCE< Storage > &a)
void sinh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > log10(const PCE< Storage > &a)
void plusEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
int n
KOKKOS_INLINE_FUNCTION PCE< Storage > asin(const PCE< Storage > &a)
ordinal_type solve(ordinal_type s, ordinal_type nrhs)
Solve linear system.
void exp(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > cos(const PCE< Storage > &a)