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 // @HEADER
2 // *****************************************************************************
3 // Stokhos Package
4 //
5 // Copyright 2009 NTESS and the Stokhos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #include "Teuchos_Assert.hpp"
12 
13 template <typename ordinal_type, typename value_type>
20  basis(basis_),
21  Bij(Bij_),
22  Cijk(Cijk_),
23  Dijk(Dijk_),
24  sz(basis->size()),
25  A(2*sz,2*sz),
26  B(2*sz,2),
27  piv(2*sz),
28  lapack()
29 {
30 }
31 
32 extern "C" {
33  double dlange_(char*, int*, int*, double*, int*, double*);
34 }
35 
36 template <typename ordinal_type, typename value_type>
40 {
41  if (s == 0 || nrhs == 0)
42  return 0;
43 
44  ordinal_type info;
45 // lapack.GESV(s, nrhs, A.values(), A.numRows(), &(piv[0]), b.values(),
46 // b.numRows(), &info);
47  lapack.GETRF(s, s, A.values(), A.numRows(), &(piv[0]), &info);
48  value_type norm, rcond;
51  norm = 1.0;
52  ordinal_type n = A.numRows();
53  char t = '1';
54  norm = dlange_(&t, &s, &s, A.values(), &n, &work[0]);
55  lapack.GECON('1', s, A.values(), A.numRows(), norm, &rcond, &work[0],
56  &iwork[0], &info);
57  //std::cout << "condition number = " << 1.0/rcond << std::endl;
58  lapack.GETRS('N', s, nrhs, A.values(), A.numRows(), &(piv[0]), B.values(),
59  B.numRows(), &info);
60  return info;
61 }
62 
63 template <typename ordinal_type, typename value_type>
64 void
69 {
70  ordinal_type pc = a.size();
71  if (c.size() != pc)
72  c.resize(pc);
73 
74  value_type* cc = c.coeff();
75  const value_type* ca = a.coeff();
76 
77  for (ordinal_type i=0; i<pc; i++)
78  cc[i] = -ca[i];
79 }
80 
81 template <typename ordinal_type, typename value_type>
82 void
85  const value_type& val)
86 {
87  c[0] += val;
88 }
89 
90 template <typename ordinal_type, typename value_type>
91 void
94  const value_type& val)
95 {
96  c[0] -= val;
97 }
98 
99 template <typename ordinal_type, typename value_type>
100 void
103  const value_type& val)
104 {
105  ordinal_type pc = c.size();
106  value_type* cc = c.coeff();
107  for (ordinal_type i=0; i<pc; i++)
108  cc[i] *= val;
109 }
110 
111 template <typename ordinal_type, typename value_type>
112 void
115  const value_type& val)
116 {
117  ordinal_type pc = c.size();
118  value_type* cc = c.coeff();
119  for (ordinal_type i=0; i<pc; i++)
120  cc[i] /= val;
121 }
122 
123 template <typename ordinal_type, typename value_type>
124 void
129 {
130  ordinal_type xp = x.size();
131  if (c.size() < xp)
132  c.resize(xp);
133 
134  value_type* cc = c.coeff();
135  const value_type* xc = x.coeff();
136  for (ordinal_type i=0; i<xp; i++)
137  cc[i] += xc[i];
138 }
139 
140 template <typename ordinal_type, typename value_type>
141 void
146 {
147  ordinal_type xp = x.size();
148  if (c.size() < xp)
149  c.resize(xp);
150 
151  value_type* cc = c.coeff();
152  const value_type* xc = x.coeff();
153  for (ordinal_type i=0; i<xp; i++)
154  cc[i] -= xc[i];
155 }
156 
157 template <typename ordinal_type, typename value_type>
158 void
163 {
164  ordinal_type p = c.size();
165  ordinal_type xp = x.size();
166  ordinal_type pc;
167  if (p > 1 && xp > 1)
168  pc = sz;
169  else
170  pc = p*xp;
171  TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
172  "Stokhos::DerivOrthogPolyExpansion::timesEqual()" <<
173  ": Expansion size (" << sz <<
174  ") is too small for computation.");
175  if (c.size() != pc)
176  c.resize(pc);
177 
178  value_type* cc = c.coeff();
179  const value_type* xc = x.coeff();
180 
181  if (p > 1 && xp > 1) {
182  // Copy c coefficients into temporary array
184  value_type tmp, cijk;
185  ordinal_type i,j;
186  for (ordinal_type k=0; k<pc; k++) {
187  tmp = value_type(0.0);
188  ordinal_type n = Cijk->num_values(k);
189  for (ordinal_type l=0; l<n; l++) {
190  Cijk->value(k,l,i,j,cijk);
191  if (i < p && j < xp)
192  tmp += cijk*tc[i]*xc[j];
193  }
194  cc[k] = tmp / basis->norm_squared(k);
195  }
196  }
197  else if (p > 1) {
198  for (ordinal_type i=0; i<p; i++)
199  cc[i] *= xc[0];
200  }
201  else if (xp > 1) {
202  for (ordinal_type i=1; i<xp; i++)
203  cc[i] = cc[0]*xc[i];
204  cc[0] *= xc[0];
205  }
206  else {
207  cc[0] *= xc[0];
208  }
209 }
210 
211 template <typename ordinal_type, typename value_type>
212 void
216 {
217  const char* func = "Stokhos::DerivOrthogPolyExpansion::divide()";
218  ordinal_type p = c.size();
219  ordinal_type xp = x.size();
220  ordinal_type pc;
221  if (xp > 1)
222  pc = sz;
223  else
224  pc = p;
225  TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
226  "Stokhos::DerivOrthogPolyExpansion::divideEqual()" <<
227  ": Expansion size (" << sz <<
228  ") is too small for computation.");
229  if (c.size() != pc)
230  c.resize(pc);
231 
232  value_type* cc = c.coeff();
233  const value_type* xc = x.coeff();
234 
235  if (xp > 1) {
236 
237  // Fill A
239  ordinal_type i,j;
240  for (ordinal_type i=0; i<pc; i++)
241  for (ordinal_type j=0; j<pc; j++)
242  A(i,j) = 0.0;
243  for (ordinal_type k=0; k<xp; k++) {
244  ordinal_type n = Cijk->num_values(k);
245  for (ordinal_type l=0; l<n; l++) {
246  Cijk->value(k,l,i,j,cijk);
247  A(i,j) += cijk*xc[k]/basis->norm_squared(i);
248  }
249  }
250 
251  // Fill B
252  for (ordinal_type i=0; i<p; i++)
253  B(i,0) = cc[i];
254  for (ordinal_type i=p; i<pc; i++)
255  B(i,0) = value_type(0.0);
256 
257  // Solve system
258  int info = solve(pc, 1);
259 
260  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
261  func << ": Argument " << info
262  << " for solve had illegal value");
263  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
264  func << ": Diagonal entry " << info
265  << " in LU factorization is exactly zero");
266 
267  // Get coefficients
268  for (ordinal_type i=0; i<pc; i++)
269  cc[i] = B(i,0);
270  }
271  else {
272  for (ordinal_type i=0; i<p; i++)
273  cc[i] /= xc[0];
274  }
275 }
276 
277 template <typename ordinal_type, typename value_type>
278 void
283 {
284  ordinal_type pa = a.size();
285  ordinal_type pb = b.size();
286  ordinal_type pc = pa > pb ? pa : pb;
287  if (c.size() != pc)
288  c.resize(pc);
289 
290  const value_type* ca = a.coeff();
291  const value_type* cb = b.coeff();
292  value_type* cc = c.coeff();
293 
294  if (pa > pb) {
295  for (ordinal_type i=0; i<pb; i++)
296  cc[i] = ca[i] + cb[i];
297  for (ordinal_type i=pb; i<pc; i++)
298  cc[i] = ca[i];
299  }
300  else {
301  for (ordinal_type i=0; i<pa; i++)
302  cc[i] = ca[i] + cb[i];
303  for (ordinal_type i=pa; i<pc; i++)
304  cc[i] = cb[i];
305  }
306 }
307 
308 template <typename ordinal_type, typename value_type>
309 void
312  const value_type& a,
314 {
315  ordinal_type pc = b.size();
316  if (c.size() != pc)
317  c.resize(pc);
318 
319  const value_type* cb = b.coeff();
320  value_type* cc = c.coeff();
321 
322  cc[0] = a + cb[0];
323  for (ordinal_type i=1; i<pc; i++)
324  cc[i] = cb[i];
325 }
326 
327 template <typename ordinal_type, typename value_type>
328 void
332  const value_type& b)
333 {
334  ordinal_type pc = a.size();
335  if (c.size() != pc)
336  c.resize(pc);
337 
338  const value_type* ca = a.coeff();
339  value_type* cc = c.coeff();
340 
341  cc[0] = ca[0] + b;
342  for (ordinal_type i=1; i<pc; i++)
343  cc[i] = ca[i];
344 }
345 
346 template <typename ordinal_type, typename value_type>
347 void
352 {
353  ordinal_type pa = a.size();
354  ordinal_type pb = b.size();
355  ordinal_type pc = pa > pb ? pa : pb;
356  if (c.size() != pc)
357  c.resize(pc);
358 
359  const value_type* ca = a.coeff();
360  const value_type* cb = b.coeff();
361  value_type* cc = c.coeff();
362 
363  if (pa > pb) {
364  for (ordinal_type i=0; i<pb; i++)
365  cc[i] = ca[i] - cb[i];
366  for (ordinal_type i=pb; i<pc; i++)
367  cc[i] = ca[i];
368  }
369  else {
370  for (ordinal_type i=0; i<pa; i++)
371  cc[i] = ca[i] - cb[i];
372  for (ordinal_type i=pa; i<pc; i++)
373  cc[i] = -cb[i];
374  }
375 }
376 
377 template <typename ordinal_type, typename value_type>
378 void
381  const value_type& a,
383 {
384  ordinal_type pc = b.size();
385  if (c.size() != pc)
386  c.resize(pc);
387 
388  const value_type* cb = b.coeff();
389  value_type* cc = c.coeff();
390 
391  cc[0] = a - cb[0];
392  for (ordinal_type i=1; i<pc; i++)
393  cc[i] = -cb[i];
394 }
395 
396 template <typename ordinal_type, typename value_type>
397 void
401  const value_type& b)
402 {
403  ordinal_type pc = a.size();
404  if (c.size() != pc)
405  c.resize(pc);
406 
407  const value_type* ca = a.coeff();
408  value_type* cc = c.coeff();
409 
410  cc[0] = ca[0] - b;
411  for (ordinal_type i=1; i<pc; i++)
412  cc[i] = ca[i];
413 }
414 
415 template <typename ordinal_type, typename value_type>
416 void
421 {
422  ordinal_type pa = a.size();
423  ordinal_type pb = b.size();
424  ordinal_type pc;
425  if (pa > 1 && pb > 1)
426  pc = sz;
427  else
428  pc = pa*pb;
429  TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
430  "Stokhos::DerivOrthogPolyExpansion::times()" <<
431  ": Expansion size (" << sz <<
432  ") is too small for computation.");
433  if (c.size() != pc)
434  c.resize(pc);
435 
436  const value_type* ca = a.coeff();
437  const value_type* cb = b.coeff();
438  value_type* cc = c.coeff();
439 
440  if (pa > 1 && pb > 1) {
441  value_type tmp, cijk;
442  ordinal_type i,j;
443  for (ordinal_type k=0; k<pc; k++) {
444  tmp = value_type(0.0);
445  ordinal_type n = Cijk->num_values(k);
446  for (ordinal_type l=0; l<n; l++) {
447  Cijk->value(k,l,i,j,cijk);
448  if (i < pa && j < pb)
449  tmp += cijk*ca[i]*cb[j];
450  }
451  cc[k] = tmp / basis->norm_squared(k);
452  }
453  }
454  else if (pa > 1) {
455  for (ordinal_type i=0; i<pc; i++)
456  cc[i] = ca[i]*cb[0];
457  }
458  else if (pb > 1) {
459  for (ordinal_type i=0; i<pc; i++)
460  cc[i] = ca[0]*cb[i];
461  }
462  else {
463  cc[0] = ca[0]*cb[0];
464  }
465 }
466 
467 template <typename ordinal_type, typename value_type>
468 void
471  const value_type& a,
473 {
474  ordinal_type pc = b.size();
475  if (c.size() != pc)
476  c.resize(pc);
477 
478  const value_type* cb = b.coeff();
479  value_type* cc = c.coeff();
480 
481  for (ordinal_type i=0; i<pc; i++)
482  cc[i] = a*cb[i];
483 }
484 
485 template <typename ordinal_type, typename value_type>
486 void
490  const value_type& b)
491 {
492  ordinal_type pc = a.size();
493  if (c.size() != pc)
494  c.resize(pc);
495 
496  const value_type* ca = a.coeff();
497  value_type* cc = c.coeff();
498 
499  for (ordinal_type i=0; i<pc; i++)
500  cc[i] = ca[i]*b;
501 }
502 
503 template <typename ordinal_type, typename value_type>
504 void
509 {
510  const char* func = "Stokhos::DerivOrthogPolyExpansion::divide()";
511  ordinal_type pa = a.size();
512  ordinal_type pb = b.size();
513  ordinal_type pc;
514  if (pb > 1)
515  pc = sz;
516  else
517  pc = pa;
518  TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
519  "Stokhos::DerivOrthogPolyExpansion::divide()" <<
520  ": Expansion size (" << sz <<
521  ") is too small for computation.");
522  if (c.size() != pc)
523  c.resize(pc);
524 
525  const value_type* ca = a.coeff();
526  const value_type* cb = b.coeff();
527  value_type* cc = c.coeff();
528 
529  if (pb > 1) {
530 
531  // Fill A
533  ordinal_type i,j;
534  for (ordinal_type i=0; i<pc; i++)
535  for (ordinal_type j=0; j<pc; j++)
536  A(i,j) = 0.0;
537 
538  for (ordinal_type k=0; k<pb; k++) {
539  ordinal_type n = Cijk->num_values(k);
540  for (ordinal_type l=0; l<n; l++) {
541  Cijk->value(k,l,i,j,cijk);
542  A(i,j) += cijk*cb[k] / basis->norm_squared(i);
543  }
544  }
545 
546  // Fill B
547  for (ordinal_type i=0; i<pa; i++)
548  B(i,0) = ca[i];
549  for (ordinal_type i=pa; i<pc; i++)
550  B(i,0) = value_type(0.0);
551 
552  // Solve system
553  int info = solve(pc, 1);
554 
555  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
556  func << ": Argument " << info
557  << " for solve had illegal value");
558  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
559  func << ": Diagonal entry " << info
560  << " in LU factorization is exactly zero");
561 
562  // Get coefficients
563  for (ordinal_type i=0; i<pc; i++)
564  cc[i] = B(i,0);
565  }
566  else {
567  for (ordinal_type i=0; i<pa; i++)
568  cc[i] = ca[i]/cb[0];
569  }
570 }
571 
572 template <typename ordinal_type, typename value_type>
573 void
576  const value_type& a,
578 {
579  const char* func = "Stokhos::DerivOrthogPolyExpansion::divide()";
580  ordinal_type pb = b.size();
581  ordinal_type pc;
582  if (pb > 1)
583  pc = sz;
584  else
585  pc = 1;
586  if (c.size() != pc)
587  c.resize(pc);
588 
589  const value_type* cb = b.coeff();
590  value_type* cc = c.coeff();
591 
592  if (pb > 1) {
593 
594  // Fill A
596  ordinal_type i,j;
597  for (ordinal_type i=0; i<pc; i++)
598  for (ordinal_type j=0; j<pc; j++)
599  A(i,j) = 0.0;
600 
601  for (ordinal_type k=0; k<pb; k++) {
602  ordinal_type n = Cijk->num_values(k);
603  for (ordinal_type l=0; l<n; l++) {
604  Cijk->value(k,l,i,j,cijk);
605  A(i,j) += cijk*cb[k] / basis->norm_squared(i);
606  }
607  }
608 
609  // Fill B
610  B(0,0) = a;
611  for (ordinal_type i=1; i<pc; i++)
612  B(i,0) = value_type(0.0);
613 
614  // Solve system
615  int info = solve(pc, 1);
616 
617  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
618  func << ": Argument " << info
619  << " for solve had illegal value");
620  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
621  func << ": Diagonal entry " << info
622  << " in LU factorization is exactly zero");
623 
624  // Get coefficients
625  for (ordinal_type i=0; i<pc; i++)
626  cc[i] = B(i,0);
627  }
628  else
629  cc[0] = a / cb[0];
630 }
631 
632 template <typename ordinal_type, typename value_type>
633 void
637  const value_type& b)
638 {
639  ordinal_type pc = a.size();
640  if (c.size() != pc)
641  c.resize(pc);
642 
643  const value_type* ca = a.coeff();
644  value_type* cc = c.coeff();
645 
646  for (ordinal_type i=0; i<pc; i++)
647  cc[i] = ca[i]/b;
648 }
649 
650 template <typename ordinal_type, typename value_type>
651 void
655 {
656  const char* func = "Stokhos::DerivOrthogPolyExpansion::exp()";
657  ordinal_type pa = a.size();
658  ordinal_type pc;
659  if (pa > 1)
660  pc = sz;
661  else
662  pc = 1;
663  if (c.size() != pc)
664  c.resize(pc);
665 
666  const value_type* ca = a.coeff();
667  value_type* cc = c.coeff();
668 
669  if (pa > 1) {
670  value_type psi_0 = basis->evaluateZero(0);
671 
672  // Fill A and B
673  for (ordinal_type i=1; i<pc; i++) {
674  B(i-1,0) = 0.0;
675  for (ordinal_type j=1; j<pc; j++) {
676  A(i-1,j-1) = (*Bij)(i-1,j);
677  for (ordinal_type k=1; k<pa; k++)
678  A(i-1,j-1) -= ca[k]*(*Dijk)(i-1,j,k);
679  B(i-1,0) += ca[j]*(*Bij)(i-1,j);
680  }
681  B(i-1,0) *= psi_0;
682  }
683 
684  // Solve system
685  int info = solve(pc-1, 1);
686 
687  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
688  func << ": Argument " << info
689  << " for solve had illegal value");
690  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
691  func << ": Diagonal entry " << info
692  << " in LU factorization is exactly zero");
693 
694  // Compute order-0 coefficient
695  value_type s = psi_0 * ca[0];
696  value_type t = psi_0;
697  for (ordinal_type i=1; i<pc; i++) {
698  s += basis->evaluateZero(i) * ca[i];
699  t += basis->evaluateZero(i) * B(i-1,0);
700  }
701  s = std::exp(s);
702  cc[0] = (s/t);
703 
704  // Compute remaining coefficients
705  for (ordinal_type i=1; i<pc; i++)
706  cc[i] = B(i-1,0) * cc[0];
707  }
708  else
709  cc[0] = std::exp(ca[0]);
710 }
711 
712 template <typename ordinal_type, typename value_type>
713 void
717 {
718  const char* func = "Stokhos::DerivOrthogPolyExpansion::log()";
719  ordinal_type pa = a.size();
720  ordinal_type pc;
721  if (pa > 1)
722  pc = sz;
723  else
724  pc = 1;
725  if (c.size() != pc)
726  c.resize(pc);
727 
728  const value_type* ca = a.coeff();
729  value_type* cc = c.coeff();
730 
731  if (pa > 1) {
732  value_type psi_0 = basis->evaluateZero(0);
733 
734  // Fill A and B
735  for (ordinal_type i=1; i<pc; i++) {
736  B(i-1,0) = 0.0;
737  for (ordinal_type j=1; j<pc; j++) {
738  A(i-1,j-1) = 0.0;
739  for (ordinal_type k=0; k<pa; k++)
740  A(i-1,j-1) += ca[k]*(*Dijk)(i-1,k,j);
741  B(i-1,0) += ca[j]*(*Bij)(i-1,j);
742  }
743  }
744 
745  // Solve system
746  int info = solve(pc-1, 1);
747 
748  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
749  func << ": Argument " << info
750  << " for solve had illegal value");
751  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
752  func << ": Diagonal entry " << info
753  << " in LU factorization is exactly zero");
754 
755  // Compute order-0 coefficient
756  value_type s = psi_0 * ca[0];
757  value_type t = value_type(0.0);
758  for (ordinal_type i=1; i<pc; i++) {
759  s += basis->evaluateZero(i) * ca[i];
760  t += basis->evaluateZero(i) * B(i-1,0);
761  }
762  cc[0] = (std::log(s) - t) / psi_0;
763 
764  // Compute remaining coefficients
765  for (ordinal_type i=1; i<pc; i++)
766  cc[i] = B(i-1,0);
767  }
768  else
769  cc[0] = std::log(ca[0]);
770 }
771 
772 template <typename ordinal_type, typename value_type>
773 void
777 {
778  if (a.size() > 1) {
779  log(c,a);
780  divide(c,c,std::log(10.0));
781  }
782  else {
783  if (c.size() != 1)
784  c.resize(1);
785  c[0] = std::log10(a[0]);
786  }
787 }
788 
789 template <typename ordinal_type, typename value_type>
790 void
794 {
795  if (a.size() > 1) {
796  log(c,a);
797  timesEqual(c,value_type(0.5));
798  exp(c,c);
799  }
800  else {
801  if (c.size() != 1)
802  c.resize(1);
803  c[0] = std::sqrt(a[0]);
804  }
805 }
806 
807 template <typename ordinal_type, typename value_type>
808 void
812 {
813  if (a.size() > 1) {
814  log(c,a);
815  timesEqual(c,value_type(1.0/3.0));
816  exp(c,c);
817  }
818  else {
819  if (c.size() != 1)
820  c.resize(1);
821  c[0] = std::cbrt(a[0]);
822  }
823 }
824 
825 template <typename ordinal_type, typename value_type>
826 void
831 {
832  if (a.size() > 1 || b.size() > 1) {
833  log(c,a);
834  timesEqual(c,b);
835  exp(c,c);
836  }
837  else {
838  if (c.size() != 1)
839  c.resize(1);
840  c[0] = std::pow(a[0], b[0]);
841  }
842 }
843 
844 template <typename ordinal_type, typename value_type>
845 void
848  const value_type& a,
850 {
851  if (b.size() > 1) {
852  times(c,std::log(a),b);
853  exp(c,c);
854  }
855  else {
856  if (c.size() != 1)
857  c.resize(1);
858  c[0] = std::pow(a, b[0]);
859  }
860 }
861 
862 template <typename ordinal_type, typename value_type>
863 void
867  const value_type& b)
868 {
869  if (a.size() > 1) {
870  log(c,a);
871  timesEqual(c,b);
872  exp(c,c);
873  }
874  else {
875  if (c.size() != 1)
876  c.resize(1);
877  c[0] = std::pow(a[0], b);
878  }
879 }
880 
881 template <typename ordinal_type, typename value_type>
882 void
887 {
888  const char* func = "Stokhos::DerivOrthogPolyExpansion::sincos()";
889  ordinal_type pa = a.size();
890  ordinal_type pc;
891  if (pa > 1)
892  pc = sz;
893  else
894  pc = 1;
895  if (c.size() != pc)
896  c.resize(pc);
897  if (s.size() != pc)
898  s.resize(pc);
899 
900  const value_type* ca = a.coeff();
901  value_type* cs = s.coeff();
902  value_type* cc = c.coeff();
903 
904  if (pa > 1) {
905  value_type psi_0 = basis->evaluateZero(0);
906  ordinal_type offset = pc-1;
907 
908  // Fill A and b
909  B.putScalar(value_type(0.0));
910  value_type tmp, tmp2;
911  for (ordinal_type i=1; i<pc; i++) {
912  tmp2 = value_type(0.0);
913  for (ordinal_type j=1; j<pc; j++) {
914  tmp = (*Bij)(i-1,j);
915  A(i-1,j-1) = tmp;
916  A(i-1+offset,j-1+offset) = tmp;
917  tmp = value_type(0.0);
918  for (ordinal_type k=1; k<pa; k++)
919  tmp += ca[k]*(*Dijk)(i-1,j,k);
920  A(i-1+offset,j-1) = tmp;
921  A(i-1,j-1+offset) = -tmp;
922  tmp2 += ca[j]*(*Bij)(i-1,j);
923  }
924  B(i-1,0) = tmp2*psi_0;
925  B(i-1+offset,1) = tmp2*psi_0;
926  }
927 
928  // Solve system
929  int info = solve(2*pc-2, 2);
930 
931  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
932  func << ": Argument " << info
933  << " for solve had illegal value");
934  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
935  func << ": Diagonal entry " << info
936  << " in LU factorization is exactly zero");
937 
938  // Compute order-0 coefficients
939  value_type t = psi_0 * ca[0];
940  value_type a00 = psi_0;
941  value_type a01 = value_type(0.0);
942  value_type a10 = value_type(0.0);
943  value_type a11 = psi_0;
944  value_type b0 = B(0,0);
945  value_type b1 = B(1,0);
946  for (ordinal_type i=1; i<pc; i++) {
947  t += basis->evaluateZero(i) * ca[i];
948  a00 -= basis->evaluateZero(i) * B(i-1,1);
949  a01 += basis->evaluateZero(i) * B(i-1,0);
950  a10 -= basis->evaluateZero(i) * B(i-1+offset,1);
951  a11 += basis->evaluateZero(i) * B(i-1+offset,0);
952  }
953  A(0,0) = a00;
954  A(0,1) = a01;
955  A(1,0) = a10;
956  A(1,1) = a11;
957  B(0,0) = std::sin(t);
958  B(1,0) = std::cos(t);
959 
960  info = solve(2, 1);
961 
962  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
963  func << ": Argument " << info
964  << " for (2x2) solve had illegal value");
965  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
966  func << ": Diagonal entry " << info
967  << " in (2x2) LU factorization is exactly zero");
968  cs[0] = B(0,0);
969  cc[0] = B(1,0);
970 
971  // Compute remaining coefficients
972  B(0,0) = b0;
973  B(1,0) = b1;
974  for (ordinal_type i=1; i<pc; i++) {
975  cs[i] = cc[0]*B(i-1,0) - cs[0]*B(i-1,1);
976  cc[i] = cc[0]*B(i-1+offset,0) - cs[0]*B(i-1+offset,1);
977  }
978  }
979  else {
980  cs[0] = std::sin(ca[0]);
981  cc[0] = std::cos(ca[0]);
982  }
983 }
984 
985 template <typename ordinal_type, typename value_type>
986 void
990 {
991  if (a.size() > 1) {
993  sincos(s, c, a);
994  }
995  else {
996  if (s.size() != 1)
997  s.resize(1);
998  s[0] = std::sin(a[0]);
999  }
1000 }
1001 
1002 template <typename ordinal_type, typename value_type>
1003 void
1007 {
1008  if (a.size() > 1) {
1010  sincos(s, c, a);
1011  }
1012  else {
1013  if (c.size() != 1)
1014  c.resize(1);
1015  c[0] = std::cos(a[0]);
1016  }
1017 }
1018 
1019 template <typename ordinal_type, typename value_type>
1020 void
1024 {
1025  if (a.size() > 1) {
1027  sincos(t, c, a);
1028  divideEqual(t,c);
1029  }
1030  else {
1031  if (t.size() != 1)
1032  t.resize(1);
1033  t[0] = std::tan(a[0]);
1034  }
1035 }
1036 
1037 template <typename ordinal_type, typename value_type>
1038 void
1043 {
1044  const char* func = "Stokhos::DerivOrthogPolyExpansion::sinhcosh()";
1045  ordinal_type pa = a.size();
1046  ordinal_type pc;
1047  if (pa > 1)
1048  pc = sz;
1049  else
1050  pc = 1;
1051  if (c.size() != pc)
1052  c.resize(pc);
1053  if (s.size() != pc)
1054  s.resize(pc);
1055 
1056  const value_type* ca = a.coeff();
1057  value_type* cs = s.coeff();
1058  value_type* cc = c.coeff();
1059 
1060  if (pa > 1) {
1061  value_type psi_0 = basis->evaluateZero(0);
1062  ordinal_type offset = pc-1;
1063 
1064  // Fill A and b
1065  B.putScalar(value_type(0.0));
1066  value_type tmp, tmp2;
1067  for (ordinal_type i=1; i<pc; i++) {
1068  tmp2 = value_type(0.0);
1069  for (ordinal_type j=1; j<pc; j++) {
1070  tmp = (*Bij)(i-1,j);
1071  A(i-1,j-1) = tmp;
1072  A(i-1+offset,j-1+offset) = tmp;
1073  tmp = value_type(0.0);
1074  for (ordinal_type k=1; k<pa; k++)
1075  tmp += ca[k]*(*Dijk)(i-1,j,k);
1076  A(i-1+offset,j-1) = -tmp;
1077  A(i-1,j-1+offset) = -tmp;
1078  tmp2 += ca[j]*(*Bij)(i-1,j);
1079  }
1080  B(i-1,0) = tmp2*psi_0;
1081  B(i-1+offset,1) = tmp2*psi_0;
1082  }
1083 
1084  // Solve system
1085  int info = solve(2*pc-2, 2);
1086 
1087  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
1088  func << ": Argument " << info
1089  << " for solve had illegal value");
1090  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
1091  func << ": Diagonal entry " << info
1092  << " in LU factorization is exactly zero");
1093 
1094  // Compute order-0 coefficients
1095  value_type t = psi_0 * ca[0];
1096  value_type a00 = psi_0;
1097  value_type a01 = value_type(0.0);
1098  value_type a10 = value_type(0.0);
1099  value_type a11 = psi_0;
1100  value_type b0 = B(0,0);
1101  value_type b1 = B(1,0);
1102  for (ordinal_type i=1; i<pc; i++) {
1103  t += basis->evaluateZero(i) * ca[i];
1104  a00 += basis->evaluateZero(i) * B(i-1,1);
1105  a01 += basis->evaluateZero(i) * B(i-1,0);
1106  a10 += basis->evaluateZero(i) * B(i-1+offset,1);
1107  a11 += basis->evaluateZero(i) * B(i-1+offset,0);
1108  }
1109  A(0,0) = a00;
1110  A(0,1) = a01;
1111  A(1,0) = a10;
1112  A(1,1) = a11;
1113  B(0,0) = std::sinh(t);
1114  B(1,0) = std::cosh(t);
1115  info = solve(2, 1);
1116  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
1117  func << ": Argument " << info
1118  << " for (2x2) solve had illegal value");
1119  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
1120  func << ": Diagonal entry " << info
1121  << " in (2x2) LU factorization is exactly zero");
1122  cs[0] = B(0,0);
1123  cc[0] = B(1,0);
1124 
1125  // Compute remaining coefficients
1126  B(0,0) = b0;
1127  B(1,0) = b1;
1128  for (ordinal_type i=1; i<pc; i++) {
1129  cs[i] = cc[0]*B(i-1,0) + cs[0]*B(i-1,1);
1130  cc[i] = cc[0]*B(i-1+offset,0) + cs[0]*B(i-1+offset,1);
1131  }
1132  }
1133  else {
1134  cs[0] = std::sinh(ca[0]);
1135  cc[0] = std::cosh(ca[0]);
1136  }
1137 }
1138 
1139 template <typename ordinal_type, typename value_type>
1140 void
1144 {
1145  if (a.size() > 1) {
1147  sinhcosh(s, c, a);
1148  }
1149  else {
1150  if (s.size() != 1)
1151  s.resize(1);
1152  s[0] = std::sinh(a[0]);
1153  }
1154 }
1155 
1156 template <typename ordinal_type, typename value_type>
1157 void
1161 {
1162  if (a.size() > 1) {
1164  sinhcosh(s, c, a);
1165  }
1166  else {
1167  if (c.size() != 1)
1168  c.resize(1);
1169  c[0] = std::cosh(a[0]);
1170  }
1171 }
1172 
1173 template <typename ordinal_type, typename value_type>
1174 void
1178 {
1179  if (a.size() > 1) {
1181  sinhcosh(t, c, a);
1182  divideEqual(t,c);
1183  }
1184  else {
1185  if (t.size() != 1)
1186  t.resize(1);
1187  t[0] = std::tanh(a[0]);
1188  }
1189 }
1190 
1191 template <typename ordinal_type, typename value_type>
1192 template <typename OpT>
1193 void
1195 quad(const OpT& quad_func,
1199 {
1200  const char* func = "Stokhos::DerivOrthogPolyExpansion::quad()";
1201  ordinal_type pa = a.size();
1202  ordinal_type pb = b.size();
1203  ordinal_type pc = sz;
1204  if (c.size() != pc)
1205  c.resize(pc);
1206 
1207  const value_type* ca = a.coeff();
1208  const value_type* cb = b.coeff();
1209  value_type* cc = c.coeff();
1210 
1211  value_type psi_0 = basis->evaluateZero(0);
1212 
1213  // Fill A and B
1214  for (ordinal_type i=1; i<pc; i++) {
1215  B(i-1,0) = 0.0;
1216  for (ordinal_type j=1; j<pc; j++) {
1217  A(i-1,j-1) = 0.0;
1218  for (ordinal_type k=0; k<pb; k++)
1219  A(i-1,j-1) += cb[k]*(*Dijk)(i-1,k,j);
1220  }
1221  for (ordinal_type j=1; j<pa; j++) {
1222  B(i-1,0) += ca[j]*(*Bij)(i-1,j);
1223  }
1224  }
1225 
1226  // Solve system
1227  int info = solve(pc-1, 1);
1228 
1229  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
1230  func << ": Argument " << info
1231  << " for solve had illegal value");
1232  TEUCHOS_TEST_FOR_EXCEPTION(info > 0, std::logic_error,
1233  func << ": Diagonal entry " << info
1234  << " in LU factorization is exactly zero");
1235 
1236  // Compute order-0 coefficient
1237  value_type s = psi_0 * ca[0];
1238  value_type t = value_type(0.0);
1239  for (ordinal_type i=1; i<pc; i++) {
1240  s += basis->evaluateZero(i) * ca[i];
1241  t += basis->evaluateZero(i) * B(i-1,0);
1242  }
1243  cc[0] = (quad_func(s) - t) / psi_0;
1244 
1245  // Get coefficients
1246  for (ordinal_type i=1; i<pc; i++)
1247  cc[i] = B(i-1,0);
1248 }
1249 
1250 template <typename ordinal_type, typename value_type>
1251 void
1255 {
1256  if (a.size() > 1) {
1257  times(c,a,a);
1258  minus(c,value_type(1.0),c);
1259  sqrt(c,c);
1260  timesEqual(c,value_type(-1.0));
1261  quad(acos_quad_func(), c, a, c);
1262  }
1263  else {
1264  if (c.size() != 1)
1265  c.resize(1);
1266  c[0] = std::acos(a[0]);
1267  }
1268 }
1269 
1270 template <typename ordinal_type, typename value_type>
1271 void
1275 {
1276  if (a.size() > 1) {
1277  times(c,a,a);
1278  minus(c,value_type(1.0),c);
1279  sqrt(c,c);
1280  quad(asin_quad_func(), c, a, c);
1281  }
1282  else {
1283  if (c.size() != 1)
1284  c.resize(1);
1285  c[0] = std::asin(a[0]);
1286  }
1287 }
1288 
1289 template <typename ordinal_type, typename value_type>
1290 void
1294 {
1295  if (a.size() > 1) {
1296  times(c,a,a);
1297  plusEqual(c,value_type(1.0));
1298  quad(atan_quad_func(), c, a, c);
1299  }
1300  else {
1301  if (c.size() != 1)
1302  c.resize(1);
1303  c[0] = std::atan(a[0]);
1304  }
1305 }
1306 
1307 // template <typename ordinal_type, typename value_type>
1308 // Hermite<ordinal_type, value_type>
1309 // atan2(const Hermite<ordinal_type, value_type>& a,
1310 // const Hermite<ordinal_type, value_type>& b)
1311 // {
1312 // Hermite<ordinal_type, value_type> c = atan(a/b);
1313 // c.fastAccessCoeff(0) = atan2(a.coeff(0),b.coeff(0));
1314 // }
1315 
1316 // template <typename ordinal_type, typename value_type>
1317 // Hermite<ordinal_type, value_type>
1318 // atan2(const T& a,
1319 // const Hermite<ordinal_type, value_type>& b)
1320 // {
1321 // Hermite<ordinal_type, value_type> c = atan(a/b);
1322 // c.fastAccessCoeff(0) = atan2(a,b.coeff(0));
1323 // }
1324 
1325 // template <typename ordinal_type, typename value_type>
1326 // Hermite<ordinal_type, value_type>
1327 // atan2(const Hermite<ordinal_type, value_type>& a,
1328 // const T& b)
1329 // {
1330 // Hermite<ordinal_type, value_type> c = atan(a/b);
1331 // c.fastAccessCoeff(0) = atan2(a.coeff(0),b);
1332 // }
1333 
1334 template <typename ordinal_type, typename value_type>
1335 void
1339 {
1340  if (a.size() > 1) {
1341  times(c,a,a);
1342  minusEqual(c,value_type(1.0));
1343  sqrt(c,c);
1344  quad(acosh_quad_func(), c, a, c);
1345  }
1346  else {
1347  if (c.size() != 1)
1348  c.resize(1);
1349  c[0] = std::log(a[0]+std::sqrt(a[0]*a[0]-value_type(1.0)));
1350  }
1351 }
1352 
1353 template <typename ordinal_type, typename value_type>
1354 void
1358 {
1359  if (a.size() > 1) {
1360  times(c,a,a);
1361  plusEqual(c,value_type(1.0));
1362  sqrt(c,c);
1363  quad(asinh_quad_func(), c, a, c);
1364  }
1365  else {
1366  if (c.size() != 1)
1367  c.resize(1);
1368  c[0] = std::log(a[0]+std::sqrt(a[0]*a[0]+value_type(1.0)));
1369  }
1370 }
1371 
1372 template <typename ordinal_type, typename value_type>
1373 void
1377 {
1378  if (a.size() > 1) {
1379  times(c,a,a);
1380  minus(c,value_type(1.0),c);
1381  quad(atanh_quad_func(), c, a, c);
1382  }
1383  else {
1384  if (c.size() != 1)
1385  c.resize(1);
1386  c[0] = 0.5*std::log((value_type(1.0)+a[0])/(value_type(1.0)-a[0]));
1387  }
1388 }
1389 
1390 template <typename ordinal_type, typename value_type>
1391 void
1395 {
1396  if (a[0] >= 0)
1397  c = a;
1398  else
1399  unaryMinus(c,a);
1400 }
1401 
1402 template <typename ordinal_type, typename value_type>
1403 void
1407 {
1408  if (a[0] >= 0)
1409  c = a;
1410  else
1411  unaryMinus(c,a);
1412 }
1413 
1414 template <typename ordinal_type, typename value_type>
1415 void
1420 {
1421  if (a[0] >= b[0])
1422  c = a;
1423  else
1424  c = b;
1425 }
1426 
1427 template <typename ordinal_type, typename value_type>
1428 void
1431  const value_type& a,
1433 {
1434  if (a >= b[0]) {
1436  c[0] = a;
1437  }
1438  else
1439  c = b;
1440 }
1441 
1442 template <typename ordinal_type, typename value_type>
1443 void
1447  const value_type& b)
1448 {
1449  if (a[0] >= b)
1450  c = a;
1451  else {
1453  c[0] = b;
1454  }
1455 }
1456 
1457 template <typename ordinal_type, typename value_type>
1458 void
1463 {
1464  if (a[0] <= b[0])
1465  c = a;
1466  else
1467  c = b;
1468 }
1469 
1470 template <typename ordinal_type, typename value_type>
1471 void
1474  const value_type& a,
1476 {
1477  if (a <= b[0]) {
1479  c[0] = a;
1480  }
1481  else
1482  c = b;
1483 }
1484 
1485 template <typename ordinal_type, typename value_type>
1486 void
1490  const value_type& b)
1491 {
1492  if (a[0] <= b)
1493  c = a;
1494  else {
1496  c[0] = b;
1497  }
1498 }
1499 
1500 template <typename ordinal_type, typename value_type>
1501 void
1505 {
1506  ordinal_type pc = a.size();
1507 
1508  const value_type* ca = a.coeff();
1509  value_type* cc = c.coeff();
1510 
1511  for (ordinal_type i=0; i<pc; i++) {
1512  cc[i] = 0.0;
1513  for (ordinal_type j=0; j<pc; j++)
1514  cc[i] += ca[j]*(*Bij)(i,j);
1515  cc[i] /= basis->norm_squared(i);
1516  }
1517 }
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)