53 num_KL = solverParams.get<
int>(
"Number of KL Terms");
54 mean = solverParams.get<
double>(
"Mean");
55 std_dev = solverParams.get<
double>(
"Standard Deviation");
57 Teuchos::Array<double> domain_upper_bound_double;
58 Teuchos::Array<double> domain_lower_bound_double;
59 Teuchos::Array<double> correlation_length_double;
62 if (solverParams.isType<std::string>(
"Domain Upper Bounds"))
63 domain_upper_bound_double =
64 Teuchos::getArrayFromStringParameter<double>(
65 solverParams,
"Domain Upper Bounds");
67 domain_upper_bound_double =
68 solverParams.get< Teuchos::Array<double> >(
"Domain Upper Bounds");
71 Teuchos::Array<value_type> domain_upper_bound(domain_upper_bound_double.size());
72 for (
int i=0; i<domain_upper_bound.size(); i++)
73 domain_upper_bound[i]=domain_upper_bound_double[i];
76 if (solverParams.isType<std::string>(
"Domain Lower Bounds"))
77 domain_lower_bound_double =
78 Teuchos::getArrayFromStringParameter<double>(
79 solverParams,
"Domain Lower Bounds");
81 domain_lower_bound_double =
82 solverParams.get< Teuchos::Array<double> >(
"Domain Lower Bounds");
85 Teuchos::Array<value_type> domain_lower_bound(domain_lower_bound_double.size());
86 for (
int i=0; i<domain_lower_bound.size(); i++)
87 domain_lower_bound[i]=domain_lower_bound_double[i];
90 if (solverParams.isType<std::string>(
"Correlation Lengths"))
91 correlation_length_double =
92 Teuchos::getArrayFromStringParameter<double>(
93 solverParams,
"Correlation Lengths");
95 correlation_length_double =
96 solverParams.get< Teuchos::Array<double> >(
"Correlation Lengths");
99 Teuchos::Array<value_type> correlation_length(correlation_length_double.size());
100 for (
int i=0; i<correlation_length.size(); i++)
101 correlation_length[i]=correlation_length_double[i];
104 dim = domain_upper_bound.size();
105 Teuchos::Array< Teuchos::Array< one_d_eigen_pair_type > > eig_pairs(dim);
106 for (
int i=0; i<dim; i++) {
107 eig_pairs[i].resize(num_KL);
108 OneDExponentialCovarianceFunction<value_type> cov_func(
109 num_KL, domain_lower_bound[i], domain_upper_bound[i],
110 correlation_length[i], i, solverParams);
111 eig_pairs[i] = cov_func.getEigenPairs();
115 int num_prod =
static_cast<int>(std::pow(
static_cast<double>(num_KL),
116 static_cast<double>(dim)));
117 Teuchos::Array<product_eigen_pair_type> product_eig_pairs(num_prod);
118 Teuchos::Array<int> index(dim, 0);
120 Teuchos::Array<value_type> vals(dim);
121 Teuchos::Array<one_d_eigen_pair_type> eigs(dim);
122 while (cnt < num_prod) {
123 for (
int i=0; i<dim; i++) {
124 eigs[i] = eig_pairs[i][index[i]];
126 product_eig_pairs[cnt].set(eigs);
129 while (
j < dim-1 && index[
j] == num_KL) {
138 std::sort(product_eig_pairs.begin(), product_eig_pairs.end(),
139 ProductEigenPairGreater<one_d_eigen_func_type, execution_space>());
142 product_eigen_funcs =
144 product_eigen_values =
146 typename eigen_func_array_type::HostMirror host_product_eigen_funcs =
148 typename eigen_value_array_type::HostMirror host_product_eigen_values =
150 for (
int i=0; i<num_prod; ++i) {
151 host_product_eigen_values(i) = 1.0;
152 for (
int j=0;
j<dim; ++
j) {
153 host_product_eigen_values(i) *= product_eig_pairs[i].eig_pairs[
j].eig_val;
154 host_product_eigen_funcs(i,
j) =
155 product_eig_pairs[i].eig_pairs[
j].eig_func;
162template <
typename value_type,
typename execution_space>
163template <
typename po
int_type,
typename rv_type>
164KOKKOS_INLINE_FUNCTION
165typename Teuchos::PromotionTraits<typename rv_type::value_type, value_type>::promote
167evaluate(
const point_type& point,
const rv_type& random_variables)
const
169 typedef typename Teuchos::PromotionTraits<typename rv_type::value_type, value_type>::promote result_type;
170 result_type result = mean;
171 for (
int i=0; i<num_KL; ++i) {
172 result_type t = std_dev*std::sqrt(product_eigen_values(i));
173 for (
int j=0;
j<dim; ++
j)
174 t *= product_eigen_funcs(i,
j).evaluate(point[
j]);
175 result += random_variables[i]*t;
180template <
typename value_type,
typename execution_space>
181template <
typename po
int_type>
182KOKKOS_INLINE_FUNCTION
183typename Teuchos::PromotionTraits<typename point_type::value_type, value_type>::promote
187 typedef typename Teuchos::PromotionTraits<typename point_type::value_type, value_type>::promote result_type;
188 result_type result = 0.0;
189 for (
int i=0; i<num_KL; i++) {
191 for (
int j=0;
j<dim; ++
j)
192 t *= product_eigen_funcs(i,
j).evaluate(point[
j]);
193 result += product_eigen_values(i).eig_val*t*t;
195 return std::sqrt(result);
198template <
typename value_type,
typename execution_space>
199template <
typename po
int_type>
200KOKKOS_INLINE_FUNCTION
201typename Teuchos::PromotionTraits<typename point_type::value_type, value_type>::promote
205 typedef typename Teuchos::PromotionTraits<typename point_type::value_type, value_type>::promote result_type;
206 result_type t = std_dev*std::sqrt(product_eigen_values(i));
207 for (
int j=0;
j<dim; ++
j)
208 t *= product_eigen_funcs(i,
j).evaluate(point[
j]);
212template <
typename value_type,
typename execution_space>
215print(std::ostream& os)
const
217 os <<
"KL expansion using " << num_KL <<
" terms in " << dim
218 <<
" dimensions:" << std::endl;
219 os <<
"\tMean = " << mean <<
", standard deviation = " << std_dev
221 os <<
"\tEigenvalues, Eigenfunctions:" << std::endl;
222 for (
int i=0; i<num_KL; i++) {
223 os <<
"\t\t" << product_eigen_values(i) <<
", ";
224 for (
int j=0;
j<dim-1; ++
j) {
226 product_eigen_funcs(i,
j).print(os);
230 product_eigen_funcs(i,dim-1).print(os);
231 os <<
")" << std::endl;
void deep_copy(const Stokhos::CrsMatrix< ValueType, DstDevice, Layout > &dst, const Stokhos::CrsMatrix< ValueType, SrcDevice, Layout > &src)