13 template <
typename ordinal_type,
typename value_type>
33 double dlange_(
char*,
int*,
int*,
double*,
int*,
double*);
36 template <
typename ordinal_type,
typename value_type>
41 if (s == 0 || nrhs == 0)
47 lapack.GETRF(s, s,
A.values(),
A.numRows(), &(piv[0]), &info);
54 norm =
dlange_(&t, &s, &s,
A.values(), &n, &work[0]);
55 lapack.GECON(
'1', s,
A.values(),
A.numRows(), norm, &rcond, &work[0],
58 lapack.GETRS(
'N', s, nrhs,
A.values(),
A.numRows(), &(piv[0]),
B.values(),
63 template <
typename ordinal_type,
typename value_type>
81 template <
typename ordinal_type,
typename value_type>
90 template <
typename ordinal_type,
typename value_type>
99 template <
typename ordinal_type,
typename value_type>
111 template <
typename ordinal_type,
typename value_type>
123 template <
typename ordinal_type,
typename value_type>
140 template <
typename ordinal_type,
typename value_type>
157 template <
typename ordinal_type,
typename value_type>
172 "Stokhos::DerivOrthogPolyExpansion::timesEqual()" <<
173 ": Expansion size (" << sz <<
174 ") is too small for computation.");
181 if (p > 1 && xp > 1) {
190 Cijk->value(k,l,i,j,cijk);
192 tmp += cijk*tc[i]*xc[
j];
194 cc[k] = tmp / basis->norm_squared(k);
211 template <
typename ordinal_type,
typename value_type>
217 const char* func =
"Stokhos::DerivOrthogPolyExpansion::divide()";
226 "Stokhos::DerivOrthogPolyExpansion::divideEqual()" <<
227 ": Expansion size (" << sz <<
228 ") is too small for computation.");
246 Cijk->value(k,l,i,j,cijk);
247 A(i,j) += cijk*xc[k]/basis->norm_squared(i);
258 int info = solve(pc, 1);
261 func <<
": Argument " << info
262 <<
" for solve had illegal value");
264 func <<
": Diagonal entry " << info
265 <<
" in LU factorization is exactly zero");
277 template <
typename ordinal_type,
typename value_type>
296 cc[i] = ca[i] + cb[i];
302 cc[i] = ca[i] + cb[i];
308 template <
typename ordinal_type,
typename value_type>
327 template <
typename ordinal_type,
typename value_type>
346 template <
typename ordinal_type,
typename value_type>
365 cc[i] = ca[i] - cb[i];
371 cc[i] = ca[i] - cb[i];
377 template <
typename ordinal_type,
typename value_type>
396 template <
typename ordinal_type,
typename value_type>
415 template <
typename ordinal_type,
typename value_type>
425 if (pa > 1 && pb > 1)
430 "Stokhos::DerivOrthogPolyExpansion::times()" <<
431 ": Expansion size (" << sz <<
432 ") is too small for computation.");
440 if (pa > 1 && pb > 1) {
447 Cijk->value(k,l,i,j,cijk);
448 if (i < pa && j < pb)
449 tmp += cijk*ca[i]*cb[
j];
451 cc[k] = tmp / basis->norm_squared(k);
467 template <
typename ordinal_type,
typename value_type>
485 template <
typename ordinal_type,
typename value_type>
503 template <
typename ordinal_type,
typename value_type>
510 const char* func =
"Stokhos::DerivOrthogPolyExpansion::divide()";
519 "Stokhos::DerivOrthogPolyExpansion::divide()" <<
520 ": Expansion size (" << sz <<
521 ") is too small for computation.");
541 Cijk->value(k,l,i,j,cijk);
542 A(i,j) += cijk*cb[k] / basis->norm_squared(i);
553 int info = solve(pc, 1);
556 func <<
": Argument " << info
557 <<
" for solve had illegal value");
559 func <<
": Diagonal entry " << info
560 <<
" in LU factorization is exactly zero");
572 template <
typename ordinal_type,
typename value_type>
579 const char* func =
"Stokhos::DerivOrthogPolyExpansion::divide()";
604 Cijk->value(k,l,i,j,cijk);
605 A(i,j) += cijk*cb[k] / basis->norm_squared(i);
615 int info = solve(pc, 1);
618 func <<
": Argument " << info
619 <<
" for solve had illegal value");
621 func <<
": Diagonal entry " << info
622 <<
" in LU factorization is exactly zero");
632 template <
typename ordinal_type,
typename value_type>
650 template <
typename ordinal_type,
typename value_type>
656 const char* func =
"Stokhos::DerivOrthogPolyExpansion::exp()";
676 A(i-1,
j-1) = (*Bij)(i-1,
j);
678 A(i-1,
j-1) -= ca[k]*(*Dijk)(i-1,
j,k);
679 B(i-1,0) += ca[
j]*(*Bij)(i-1,
j);
685 int info = solve(pc-1, 1);
688 func <<
": Argument " << info
689 <<
" for solve had illegal value");
691 func <<
": Diagonal entry " << info
692 <<
" in LU factorization is exactly zero");
698 s += basis->evaluateZero(i) * ca[i];
699 t += basis->evaluateZero(i) *
B(i-1,0);
706 cc[i] =
B(i-1,0) * cc[0];
712 template <
typename ordinal_type,
typename value_type>
718 const char* func =
"Stokhos::DerivOrthogPolyExpansion::log()";
740 A(i-1,
j-1) += ca[k]*(*Dijk)(i-1,k,
j);
741 B(i-1,0) += ca[
j]*(*Bij)(i-1,
j);
746 int info = solve(pc-1, 1);
749 func <<
": Argument " << info
750 <<
" for solve had illegal value");
752 func <<
": Diagonal entry " << info
753 <<
" in LU factorization is exactly zero");
759 s += basis->evaluateZero(i) * ca[i];
760 t += basis->evaluateZero(i) *
B(i-1,0);
772 template <
typename ordinal_type,
typename value_type>
789 template <
typename ordinal_type,
typename value_type>
807 template <
typename ordinal_type,
typename value_type>
825 template <
typename ordinal_type,
typename value_type>
844 template <
typename ordinal_type,
typename value_type>
862 template <
typename ordinal_type,
typename value_type>
881 template <
typename ordinal_type,
typename value_type>
888 const char* func =
"Stokhos::DerivOrthogPolyExpansion::sincos()";
916 A(i-1+offset,
j-1+offset) = tmp;
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);
924 B(i-1,0) = tmp2*psi_0;
925 B(i-1+offset,1) = tmp2*psi_0;
929 int info = solve(2*pc-2, 2);
932 func <<
": Argument " << info
933 <<
" for solve had illegal value");
935 func <<
": Diagonal entry " << info
936 <<
" in LU factorization is exactly zero");
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);
963 func <<
": Argument " << info
964 <<
" for (2x2) solve had illegal value");
966 func <<
": Diagonal entry " << info
967 <<
" in (2x2) LU factorization is exactly zero");
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);
985 template <
typename ordinal_type,
typename value_type>
1002 template <
typename ordinal_type,
typename value_type>
1019 template <
typename ordinal_type,
typename value_type>
1037 template <
typename ordinal_type,
typename value_type>
1044 const char* func =
"Stokhos::DerivOrthogPolyExpansion::sinhcosh()";
1070 tmp = (*Bij)(i-1,
j);
1072 A(i-1+offset,
j-1+offset) = tmp;
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);
1080 B(i-1,0) = tmp2*psi_0;
1081 B(i-1+offset,1) = tmp2*psi_0;
1085 int info = solve(2*pc-2, 2);
1088 func <<
": Argument " << info
1089 <<
" for solve had illegal value");
1091 func <<
": Diagonal entry " << info
1092 <<
" in LU factorization is exactly zero");
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);
1117 func <<
": Argument " << info
1118 <<
" for (2x2) solve had illegal value");
1120 func <<
": Diagonal entry " << info
1121 <<
" in (2x2) LU factorization is exactly zero");
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);
1139 template <
typename ordinal_type,
typename value_type>
1156 template <
typename ordinal_type,
typename value_type>
1173 template <
typename ordinal_type,
typename value_type>
1191 template <
typename ordinal_type,
typename value_type>
1192 template <
typename OpT>
1200 const char* func =
"Stokhos::DerivOrthogPolyExpansion::quad()";
1219 A(i-1,
j-1) += cb[k]*(*Dijk)(i-1,k,
j);
1222 B(i-1,0) += ca[
j]*(*Bij)(i-1,
j);
1227 int info = solve(pc-1, 1);
1230 func <<
": Argument " << info
1231 <<
" for solve had illegal value");
1233 func <<
": Diagonal entry " << info
1234 <<
" in LU factorization is exactly zero");
1240 s += basis->evaluateZero(i) * ca[i];
1241 t += basis->evaluateZero(i) *
B(i-1,0);
1243 cc[0] = (quad_func(s) - t) / psi_0;
1250 template <
typename ordinal_type,
typename value_type>
1270 template <
typename ordinal_type,
typename value_type>
1289 template <
typename ordinal_type,
typename value_type>
1334 template <
typename ordinal_type,
typename value_type>
1353 template <
typename ordinal_type,
typename value_type>
1372 template <
typename ordinal_type,
typename value_type>
1390 template <
typename ordinal_type,
typename value_type>
1402 template <
typename ordinal_type,
typename value_type>
1414 template <
typename ordinal_type,
typename value_type>
1427 template <
typename ordinal_type,
typename value_type>
1442 template <
typename ordinal_type,
typename value_type>
1457 template <
typename ordinal_type,
typename value_type>
1470 template <
typename ordinal_type,
typename value_type>
1485 template <
typename ordinal_type,
typename value_type>
1500 template <
typename ordinal_type,
typename value_type>
1514 cc[i] += ca[
j]*(*Bij)(i,
j);
1515 cc[i] /= basis->norm_squared(i);
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.
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)
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)