44#include "Teuchos_Assert.hpp"
47template <
typename ordinal_type,
typename value_type>
51 const Teuchos::RCP<
const Teuchos::SerialDenseMatrix<ordinal_type, value_type> >& Bij_,
67 double dlange_(
char*,
int*,
int*,
double*,
int*,
double*);
70template <
typename ordinal_type,
typename value_type>
73solve(ordinal_type s, ordinal_type nrhs)
75 if (s == 0 || nrhs == 0)
81 lapack.GETRF(s, s, A.values(), A.numRows(), &(piv[0]), &info);
82 value_type norm, rcond;
83 Teuchos::Array<ordinal_type> iwork(4*s);
84 Teuchos::Array<value_type> work(4*s);
86 ordinal_type n = A.numRows();
88 norm =
dlange_(&t, &s, &s, A.values(), &n, &work[0]);
89 lapack.GECON(
'1', s, A.values(), A.numRows(), norm, &rcond, &work[0],
92 lapack.GETRS(
'N', s, nrhs, A.values(), A.numRows(), &(piv[0]), B.values(),
97template <
typename ordinal_type,
typename value_type>
104 ordinal_type pc = a.
size();
108 value_type* cc = c.
coeff();
109 const value_type* ca = a.
coeff();
111 for (ordinal_type i=0; i<pc; i++)
115template <
typename ordinal_type,
typename value_type>
119 const value_type&
val)
124template <
typename ordinal_type,
typename value_type>
128 const value_type&
val)
133template <
typename ordinal_type,
typename value_type>
137 const value_type&
val)
139 ordinal_type pc = c.
size();
140 value_type* cc = c.
coeff();
141 for (ordinal_type i=0; i<pc; i++)
145template <
typename ordinal_type,
typename value_type>
149 const value_type&
val)
151 ordinal_type pc = c.
size();
152 value_type* cc = c.
coeff();
153 for (ordinal_type i=0; i<pc; i++)
157template <
typename ordinal_type,
typename value_type>
164 ordinal_type xp = x.size();
168 value_type* cc = c.
coeff();
169 const value_type* xc = x.coeff();
170 for (ordinal_type i=0; i<xp; i++)
174template <
typename ordinal_type,
typename value_type>
181 ordinal_type xp = x.size();
185 value_type* cc = c.
coeff();
186 const value_type* xc = x.coeff();
187 for (ordinal_type i=0; i<xp; i++)
191template <
typename ordinal_type,
typename value_type>
198 ordinal_type p = c.
size();
199 ordinal_type xp = x.size();
205 TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
206 "Stokhos::DerivOrthogPolyExpansion::timesEqual()" <<
207 ": Expansion size (" << sz <<
208 ") is too small for computation.");
212 value_type* cc = c.
coeff();
213 const value_type* xc = x.coeff();
215 if (p > 1 && xp > 1) {
218 value_type tmp, cijk;
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);
226 tmp += cijk*tc[i]*xc[
j];
228 cc[k] = tmp / basis->norm_squared(k);
232 for (ordinal_type i=0; i<p; i++)
236 for (ordinal_type i=1; i<xp; i++)
245template <
typename ordinal_type,
typename value_type>
249 const OrthogPolyApprox<ordinal_type, value_type, node_type>& x)
251 const char* func =
"Stokhos::DerivOrthogPolyExpansion::divide()";
252 ordinal_type p = c.
size();
253 ordinal_type xp = x.size();
259 TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
260 "Stokhos::DerivOrthogPolyExpansion::divideEqual()" <<
261 ": Expansion size (" << sz <<
262 ") is too small for computation.");
266 value_type* cc = c.
coeff();
267 const value_type* xc = x.coeff();
274 for (ordinal_type i=0; i<pc; i++)
275 for (ordinal_type
j=0;
j<pc;
j++)
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);
286 for (ordinal_type i=0; i<p; i++)
288 for (ordinal_type i=p; i<pc; i++)
289 B(i,0) = value_type(0.0);
292 int info = solve(pc, 1);
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");
302 for (ordinal_type i=0; i<pc; i++)
306 for (ordinal_type i=0; i<p; i++)
311template <
typename ordinal_type,
typename value_type>
318 ordinal_type pa = a.
size();
319 ordinal_type pb = b.
size();
320 ordinal_type pc = pa > pb ? pa : pb;
324 const value_type* ca = a.
coeff();
325 const value_type* cb = b.
coeff();
326 value_type* cc = c.
coeff();
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++)
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++)
342template <
typename ordinal_type,
typename value_type>
349 ordinal_type pc = b.
size();
353 const value_type* cb = b.
coeff();
354 value_type* cc = c.
coeff();
357 for (ordinal_type i=1; i<pc; i++)
361template <
typename ordinal_type,
typename value_type>
368 ordinal_type pc = a.
size();
372 const value_type* ca = a.
coeff();
373 value_type* cc = c.
coeff();
376 for (ordinal_type i=1; i<pc; i++)
380template <
typename ordinal_type,
typename value_type>
387 ordinal_type pa = a.
size();
388 ordinal_type pb = b.
size();
389 ordinal_type pc = pa > pb ? pa : pb;
393 const value_type* ca = a.
coeff();
394 const value_type* cb = b.
coeff();
395 value_type* cc = c.
coeff();
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++)
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++)
411template <
typename ordinal_type,
typename value_type>
418 ordinal_type pc = b.
size();
422 const value_type* cb = b.
coeff();
423 value_type* cc = c.
coeff();
426 for (ordinal_type i=1; i<pc; i++)
430template <
typename ordinal_type,
typename value_type>
437 ordinal_type pc = a.
size();
441 const value_type* ca = a.
coeff();
442 value_type* cc = c.
coeff();
445 for (ordinal_type i=1; i<pc; i++)
449template <
typename ordinal_type,
typename value_type>
456 ordinal_type pa = a.
size();
457 ordinal_type pb = b.
size();
459 if (pa > 1 && pb > 1)
463 TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
464 "Stokhos::DerivOrthogPolyExpansion::times()" <<
465 ": Expansion size (" << sz <<
466 ") is too small for computation.");
470 const value_type* ca = a.
coeff();
471 const value_type* cb = b.
coeff();
472 value_type* cc = c.
coeff();
474 if (pa > 1 && pb > 1) {
475 value_type tmp, cijk;
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];
485 cc[k] = tmp / basis->norm_squared(k);
489 for (ordinal_type i=0; i<pc; i++)
493 for (ordinal_type i=0; i<pc; i++)
501template <
typename ordinal_type,
typename value_type>
508 ordinal_type pc = b.
size();
512 const value_type* cb = b.
coeff();
513 value_type* cc = c.
coeff();
515 for (ordinal_type i=0; i<pc; i++)
519template <
typename ordinal_type,
typename value_type>
526 ordinal_type pc = a.
size();
530 const value_type* ca = a.
coeff();
531 value_type* cc = c.
coeff();
533 for (ordinal_type i=0; i<pc; i++)
537template <
typename ordinal_type,
typename value_type>
544 const char* func =
"Stokhos::DerivOrthogPolyExpansion::divide()";
545 ordinal_type pa = a.
size();
546 ordinal_type pb = b.
size();
552 TEUCHOS_TEST_FOR_EXCEPTION(sz < pc, std::logic_error,
553 "Stokhos::DerivOrthogPolyExpansion::divide()" <<
554 ": Expansion size (" << sz <<
555 ") is too small for computation.");
559 const value_type* ca = a.
coeff();
560 const value_type* cb = b.
coeff();
561 value_type* cc = c.
coeff();
568 for (ordinal_type i=0; i<pc; i++)
569 for (ordinal_type
j=0;
j<pc;
j++)
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);
581 for (ordinal_type i=0; i<pa; i++)
583 for (ordinal_type i=pa; i<pc; i++)
584 B(i,0) = value_type(0.0);
587 int info = solve(pc, 1);
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");
597 for (ordinal_type i=0; i<pc; i++)
601 for (ordinal_type i=0; i<pa; i++)
606template <
typename ordinal_type,
typename value_type>
613 const char* func =
"Stokhos::DerivOrthogPolyExpansion::divide()";
614 ordinal_type pb = b.
size();
623 const value_type* cb = b.
coeff();
624 value_type* cc = c.
coeff();
631 for (ordinal_type i=0; i<pc; i++)
632 for (ordinal_type
j=0;
j<pc;
j++)
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);
645 for (ordinal_type i=1; i<pc; i++)
646 B(i,0) = value_type(0.0);
649 int info = solve(pc, 1);
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");
659 for (ordinal_type i=0; i<pc; i++)
666template <
typename ordinal_type,
typename value_type>
673 ordinal_type pc = a.
size();
677 const value_type* ca = a.
coeff();
678 value_type* cc = c.
coeff();
680 for (ordinal_type i=0; i<pc; i++)
684template <
typename ordinal_type,
typename value_type>
690 const char* func =
"Stokhos::DerivOrthogPolyExpansion::exp()";
691 ordinal_type pa = a.
size();
700 const value_type* ca = a.
coeff();
701 value_type* cc = c.
coeff();
704 value_type psi_0 = basis->evaluateZero(0);
707 for (ordinal_type i=1; i<pc; i++) {
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);
719 int info = solve(pc-1, 1);
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");
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);
739 for (ordinal_type i=1; i<pc; i++)
740 cc[i] = B(i-1,0) * cc[0];
743 cc[0] = std::exp(ca[0]);
746template <
typename ordinal_type,
typename value_type>
752 const char* func =
"Stokhos::DerivOrthogPolyExpansion::log()";
753 ordinal_type pa = a.
size();
762 const value_type* ca = a.
coeff();
763 value_type* cc = c.
coeff();
766 value_type psi_0 = basis->evaluateZero(0);
769 for (ordinal_type i=1; i<pc; i++) {
771 for (ordinal_type
j=1;
j<pc;
j++) {
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);
780 int info = solve(pc-1, 1);
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");
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);
796 cc[0] = (std::log(s) - t) / psi_0;
799 for (ordinal_type i=1; i<pc; i++)
803 cc[0] = std::log(ca[0]);
806template <
typename ordinal_type,
typename value_type>
814 divide(c,c,std::log(10.0));
819 c[0] = std::log10(a[0]);
823template <
typename ordinal_type,
typename value_type>
831 timesEqual(c,value_type(0.5));
837 c[0] = std::sqrt(a[0]);
841template <
typename ordinal_type,
typename value_type>
849 timesEqual(c,value_type(1.0/3.0));
855 c[0] = std::cbrt(a[0]);
859template <
typename ordinal_type,
typename value_type>
874 c[0] = std::pow(a[0], b[0]);
878template <
typename ordinal_type,
typename value_type>
886 times(c,std::log(a),b);
892 c[0] = std::pow(a, b[0]);
896template <
typename ordinal_type,
typename value_type>
911 c[0] = std::pow(a[0], b);
915template <
typename ordinal_type,
typename value_type>
922 const char* func =
"Stokhos::DerivOrthogPolyExpansion::sincos()";
923 ordinal_type pa = a.
size();
934 const value_type* ca = a.
coeff();
935 value_type* cs = s.
coeff();
936 value_type* cc = c.
coeff();
939 value_type psi_0 = basis->evaluateZero(0);
940 ordinal_type offset = pc-1;
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++) {
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);
958 B(i-1,0) = tmp2*psi_0;
959 B(i-1+offset,1) = tmp2*psi_0;
963 int info = solve(2*pc-2, 2);
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");
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);
991 B(0,0) = std::sin(t);
992 B(1,0) = std::cos(t);
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");
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);
1014 cs[0] = std::sin(ca[0]);
1015 cc[0] = std::cos(ca[0]);
1019template <
typename ordinal_type,
typename value_type>
1026 OrthogPolyApprox<ordinal_type, value_type, node_type> c(s);
1032 s[0] = std::sin(a[0]);
1036template <
typename ordinal_type,
typename value_type>
1043 OrthogPolyApprox<ordinal_type, value_type, node_type> s(c);
1049 c[0] = std::cos(a[0]);
1053template <
typename ordinal_type,
typename value_type>
1060 OrthogPolyApprox<ordinal_type, value_type, node_type> c(t);
1067 t[0] = std::tan(a[0]);
1071template <
typename ordinal_type,
typename value_type>
1078 const char* func =
"Stokhos::DerivOrthogPolyExpansion::sinhcosh()";
1079 ordinal_type pa = a.
size();
1090 const value_type* ca = a.
coeff();
1091 value_type* cs = s.
coeff();
1092 value_type* cc = c.
coeff();
1095 value_type psi_0 = basis->evaluateZero(0);
1096 ordinal_type offset = pc-1;
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);
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);
1114 B(i-1,0) = tmp2*psi_0;
1115 B(i-1+offset,1) = tmp2*psi_0;
1119 int info = solve(2*pc-2, 2);
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");
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);
1147 B(0,0) = std::sinh(t);
1148 B(1,0) = std::cosh(t);
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");
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);
1168 cs[0] = std::sinh(ca[0]);
1169 cc[0] = std::cosh(ca[0]);
1173template <
typename ordinal_type,
typename value_type>
1180 OrthogPolyApprox<ordinal_type, value_type, node_type> c(s);
1186 s[0] = std::sinh(a[0]);
1190template <
typename ordinal_type,
typename value_type>
1197 OrthogPolyApprox<ordinal_type, value_type, node_type> s(c);
1203 c[0] = std::cosh(a[0]);
1207template <
typename ordinal_type,
typename value_type>
1214 OrthogPolyApprox<ordinal_type, value_type, node_type> c(t);
1221 t[0] = std::tanh(a[0]);
1225template <
typename ordinal_type,
typename value_type>
1226template <
typename OpT>
1229quad(
const OpT& quad_func,
1234 const char* func =
"Stokhos::DerivOrthogPolyExpansion::quad()";
1235 ordinal_type pa = a.
size();
1236 ordinal_type pb = b.
size();
1237 ordinal_type pc = sz;
1241 const value_type* ca = a.
coeff();
1242 const value_type* cb = b.
coeff();
1243 value_type* cc = c.
coeff();
1245 value_type psi_0 = basis->evaluateZero(0);
1248 for (ordinal_type i=1; i<pc; i++) {
1250 for (ordinal_type
j=1;
j<pc;
j++) {
1252 for (ordinal_type k=0; k<pb; k++)
1253 A(i-1,
j-1) += cb[k]*(*Dijk)(i-1,k,
j);
1255 for (ordinal_type
j=1;
j<pa;
j++) {
1256 B(i-1,0) += ca[
j]*(*Bij)(i-1,
j);
1261 int info = solve(pc-1, 1);
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");
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);
1277 cc[0] = (quad_func(s) - t) / psi_0;
1280 for (ordinal_type i=1; i<pc; i++)
1284template <
typename ordinal_type,
typename value_type>
1292 minus(c,value_type(1.0),c);
1294 timesEqual(c,value_type(-1.0));
1300 c[0] = std::acos(a[0]);
1304template <
typename ordinal_type,
typename value_type>
1312 minus(c,value_type(1.0),c);
1319 c[0] = std::asin(a[0]);
1323template <
typename ordinal_type,
typename value_type>
1331 plusEqual(c,value_type(1.0));
1337 c[0] = std::atan(a[0]);
1368template <
typename ordinal_type,
typename value_type>
1376 minusEqual(c,value_type(1.0));
1383 c[0] = std::log(a[0]+std::sqrt(a[0]*a[0]-value_type(1.0)));
1387template <
typename ordinal_type,
typename value_type>
1395 plusEqual(c,value_type(1.0));
1402 c[0] = std::log(a[0]+std::sqrt(a[0]*a[0]+value_type(1.0)));
1406template <
typename ordinal_type,
typename value_type>
1414 minus(c,value_type(1.0),c);
1420 c[0] = 0.5*std::log((value_type(1.0)+a[0])/(value_type(1.0)-a[0]));
1424template <
typename ordinal_type,
typename value_type>
1436template <
typename ordinal_type,
typename value_type>
1448template <
typename ordinal_type,
typename value_type>
1461template <
typename ordinal_type,
typename value_type>
1465 const value_type& a,
1469 c = OrthogPolyApprox<ordinal_type, value_type, node_type>(basis);
1476template <
typename ordinal_type,
typename value_type>
1481 const value_type& b)
1486 c = OrthogPolyApprox<ordinal_type, value_type, node_type>(basis);
1491template <
typename ordinal_type,
typename value_type>
1504template <
typename ordinal_type,
typename value_type>
1508 const value_type& a,
1512 c = OrthogPolyApprox<ordinal_type, value_type, node_type>(basis);
1519template <
typename ordinal_type,
typename value_type>
1524 const value_type& b)
1529 c = OrthogPolyApprox<ordinal_type, value_type, node_type>(basis);
1534template <
typename ordinal_type,
typename value_type>
1540 ordinal_type pc = a.
size();
1542 const value_type* ca = a.
coeff();
1543 value_type* cc = c.
coeff();
1545 for (ordinal_type i=0; i<pc; i++) {
1547 for (ordinal_type
j=0;
j<pc;
j++)
1548 cc[i] += ca[
j]*(*Bij)(i,
j);
1549 cc[i] /= basis->norm_squared(i);
double dlange_(char *, int *, int *, double *, int *, double *)
Data structure storing a dense 3-tensor C(i,j,k).
Abstract base class for multivariate orthogonal polynomials that support computing double and triple ...
void tan(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)
void cbrt(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 exp(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)
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)
void cosh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
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 unaryMinus(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void plusEqual(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)
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 atan(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void cos(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void acos(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void atanh(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 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)
void minusEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
void timesEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
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 log10(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)
void sin(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
ordinal_type solve(ordinal_type s, ordinal_type nrhs)
Solve linear system.
void sinh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &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 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)
void divideEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
void fabs(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void asinh(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 sqrt(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void tanh(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void log(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
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)
Class to store coefficients of a projection onto an orthogonal polynomial basis.
void resize(ordinal_type sz)
Resize coefficient array (coefficients are preserved)
pointer coeff()
Return coefficient array.
ordinal_type size() const
Return size.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
static T * get_and_fill(int sz)
Get memory for new array of length sz and fill with zeros.