47 template <
typename ordinal_type,
typename value_type>
67 double dlange_(
char*,
int*,
int*,
double*,
int*,
double*);
70 template <
typename ordinal_type,
typename value_type>
75 if (s == 0 || nrhs == 0)
81 lapack.GETRF(s, s,
A.values(),
A.numRows(), &(piv[0]), &info);
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(),
97 template <
typename ordinal_type,
typename value_type>
115 template <
typename ordinal_type,
typename value_type>
124 template <
typename ordinal_type,
typename value_type>
133 template <
typename ordinal_type,
typename value_type>
145 template <
typename ordinal_type,
typename value_type>
157 template <
typename ordinal_type,
typename value_type>
174 template <
typename ordinal_type,
typename value_type>
191 template <
typename ordinal_type,
typename value_type>
206 "Stokhos::DerivOrthogPolyExpansion::timesEqual()" <<
207 ": Expansion size (" << sz <<
208 ") is too small for computation.");
215 if (p > 1 && xp > 1) {
224 Cijk->value(k,l,i,j,cijk);
226 tmp += cijk*tc[i]*xc[
j];
228 cc[k] = tmp / basis->norm_squared(k);
245 template <
typename ordinal_type,
typename value_type>
251 const char* func =
"Stokhos::DerivOrthogPolyExpansion::divide()";
260 "Stokhos::DerivOrthogPolyExpansion::divideEqual()" <<
261 ": Expansion size (" << sz <<
262 ") is too small for computation.");
280 Cijk->value(k,l,i,j,cijk);
281 A(i,j) += cijk*xc[k]/basis->norm_squared(i);
292 int info = solve(pc, 1);
295 func <<
": Argument " << info
296 <<
" for solve had illegal value");
298 func <<
": Diagonal entry " << info
299 <<
" in LU factorization is exactly zero");
311 template <
typename ordinal_type,
typename value_type>
330 cc[i] = ca[i] + cb[i];
336 cc[i] = ca[i] + cb[i];
342 template <
typename ordinal_type,
typename value_type>
361 template <
typename ordinal_type,
typename value_type>
380 template <
typename ordinal_type,
typename value_type>
399 cc[i] = ca[i] - cb[i];
405 cc[i] = ca[i] - cb[i];
411 template <
typename ordinal_type,
typename value_type>
430 template <
typename ordinal_type,
typename value_type>
449 template <
typename ordinal_type,
typename value_type>
459 if (pa > 1 && pb > 1)
464 "Stokhos::DerivOrthogPolyExpansion::times()" <<
465 ": Expansion size (" << sz <<
466 ") is too small for computation.");
474 if (pa > 1 && pb > 1) {
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);
501 template <
typename ordinal_type,
typename value_type>
519 template <
typename ordinal_type,
typename value_type>
537 template <
typename ordinal_type,
typename value_type>
544 const char* func =
"Stokhos::DerivOrthogPolyExpansion::divide()";
553 "Stokhos::DerivOrthogPolyExpansion::divide()" <<
554 ": Expansion size (" << sz <<
555 ") is too small for computation.");
575 Cijk->value(k,l,i,j,cijk);
576 A(i,j) += cijk*cb[k] / basis->norm_squared(i);
587 int info = solve(pc, 1);
590 func <<
": Argument " << info
591 <<
" for solve had illegal value");
593 func <<
": Diagonal entry " << info
594 <<
" in LU factorization is exactly zero");
606 template <
typename ordinal_type,
typename value_type>
613 const char* func =
"Stokhos::DerivOrthogPolyExpansion::divide()";
638 Cijk->value(k,l,i,j,cijk);
639 A(i,j) += cijk*cb[k] / basis->norm_squared(i);
649 int info = solve(pc, 1);
652 func <<
": Argument " << info
653 <<
" for solve had illegal value");
655 func <<
": Diagonal entry " << info
656 <<
" in LU factorization is exactly zero");
666 template <
typename ordinal_type,
typename value_type>
684 template <
typename ordinal_type,
typename value_type>
690 const char* func =
"Stokhos::DerivOrthogPolyExpansion::exp()";
710 A(i-1,
j-1) = (*Bij)(i-1,
j);
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);
722 func <<
": Argument " << info
723 <<
" for solve had illegal value");
725 func <<
": Diagonal entry " << info
726 <<
" in LU factorization is exactly zero");
732 s += basis->evaluateZero(i) * ca[i];
733 t += basis->evaluateZero(i) *
B(i-1,0);
740 cc[i] =
B(i-1,0) * cc[0];
746 template <
typename ordinal_type,
typename value_type>
752 const char* func =
"Stokhos::DerivOrthogPolyExpansion::log()";
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);
783 func <<
": Argument " << info
784 <<
" for solve had illegal value");
786 func <<
": Diagonal entry " << info
787 <<
" in LU factorization is exactly zero");
793 s += basis->evaluateZero(i) * ca[i];
794 t += basis->evaluateZero(i) *
B(i-1,0);
806 template <
typename ordinal_type,
typename value_type>
823 template <
typename ordinal_type,
typename value_type>
841 template <
typename ordinal_type,
typename value_type>
859 template <
typename ordinal_type,
typename value_type>
878 template <
typename ordinal_type,
typename value_type>
896 template <
typename ordinal_type,
typename value_type>
915 template <
typename ordinal_type,
typename value_type>
922 const char* func =
"Stokhos::DerivOrthogPolyExpansion::sincos()";
950 A(i-1+offset,
j-1+offset) = tmp;
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);
966 func <<
": Argument " << info
967 <<
" for solve had illegal value");
969 func <<
": Diagonal entry " << info
970 <<
" in LU factorization is exactly zero");
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);
997 func <<
": Argument " << info
998 <<
" for (2x2) solve had illegal value");
1000 func <<
": Diagonal entry " << info
1001 <<
" in (2x2) LU factorization is exactly zero");
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);
1019 template <
typename ordinal_type,
typename value_type>
1036 template <
typename ordinal_type,
typename value_type>
1053 template <
typename ordinal_type,
typename value_type>
1071 template <
typename ordinal_type,
typename value_type>
1078 const char* func =
"Stokhos::DerivOrthogPolyExpansion::sinhcosh()";
1104 tmp = (*Bij)(i-1,
j);
1106 A(i-1+offset,
j-1+offset) = tmp;
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);
1122 func <<
": Argument " << info
1123 <<
" for solve had illegal value");
1125 func <<
": Diagonal entry " << info
1126 <<
" in LU factorization is exactly zero");
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);
1151 func <<
": Argument " << info
1152 <<
" for (2x2) solve had illegal value");
1154 func <<
": Diagonal entry " << info
1155 <<
" in (2x2) LU factorization is exactly zero");
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);
1173 template <
typename ordinal_type,
typename value_type>
1190 template <
typename ordinal_type,
typename value_type>
1207 template <
typename ordinal_type,
typename value_type>
1225 template <
typename ordinal_type,
typename value_type>
1226 template <
typename OpT>
1234 const char* func =
"Stokhos::DerivOrthogPolyExpansion::quad()";
1253 A(i-1,
j-1) += cb[k]*(*Dijk)(i-1,k,
j);
1256 B(i-1,0) += ca[
j]*(*Bij)(i-1,
j);
1261 int info = solve(pc-1, 1);
1264 func <<
": Argument " << info
1265 <<
" for solve had illegal value");
1267 func <<
": Diagonal entry " << info
1268 <<
" in LU factorization is exactly zero");
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;
1284 template <
typename ordinal_type,
typename value_type>
1304 template <
typename ordinal_type,
typename value_type>
1323 template <
typename ordinal_type,
typename value_type>
1368 template <
typename ordinal_type,
typename value_type>
1387 template <
typename ordinal_type,
typename value_type>
1406 template <
typename ordinal_type,
typename value_type>
1424 template <
typename ordinal_type,
typename value_type>
1436 template <
typename ordinal_type,
typename value_type>
1448 template <
typename ordinal_type,
typename value_type>
1461 template <
typename ordinal_type,
typename value_type>
1476 template <
typename ordinal_type,
typename value_type>
1491 template <
typename ordinal_type,
typename value_type>
1504 template <
typename ordinal_type,
typename value_type>
1519 template <
typename ordinal_type,
typename value_type>
1534 template <
typename ordinal_type,
typename value_type>
1548 cc[i] += ca[
j]*(*Bij)(i,
j);
1549 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)