215#ifdef HAVE_STOKHOS_ANASAZI
216#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
217 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::KLReducedMatrixFreeOperator -- Calculation/setup of KL opeator");
220 mean = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(
221 block_ops->getCoeffPtr(0));
224 for (
int coeff=0; coeff<num_blocks; coeff++) {
225 Teuchos::RCP<const Epetra_CrsMatrix> block_coeff =
226 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>
227 (block_ops->getCoeffPtr(coeff));
229 for (
int i=0; i<mean->NumMyRows(); i++) {
231 mean->NumMyRowEntries(i, num_col);
232 for (
int j=0;
j<num_col;
j++)
233 (*block_vec_poly)[coeff][row++] = (*block_coeff)[i][
j];
237 int myPID = sg_comm->MyPID();
240 Stokhos::PCEAnasaziKL pceKL(*block_vec_poly, num_KL);
241 Teuchos::ParameterList anasazi_params = pceKL.getDefaultParams();
242 bool result = pceKL.computeKL(anasazi_params);
243 if (!result && myPID == 0)
244 std::cout <<
"KL Eigensolver did not converge!" << std::endl;
245 Teuchos::RCP<Epetra_MultiVector> evecs = pceKL.getEigenvectors();
246 Teuchos::Array<double> evals = pceKL.getEigenvalues();
249 std::cout <<
"num computed eigenvectors = "
250 << evecs->NumVectors() << std::endl;
251 double kl_tol = params->get(
"KL Tolerance", 1e-6);
253 while (num_KL_computed < evals.size() &&
254 std::sqrt(evals[num_KL_computed]/evals[0]) > kl_tol)
256 if (num_KL_computed == evals.size() && myPID == 0)
257 std::cout <<
"Can't achieve KL tolerance " << kl_tol
258 <<
". Smallest eigenvalue / largest eigenvalue = "
259 << std::sqrt(evals[num_KL_computed-1]/evals[0]) << std::endl;
261 std::cout <<
"num KL eigenvectors = " << num_KL_computed << std::endl;
264 dot_products.resize(num_KL_computed);
265 for (
int rv=0; rv < num_KL_computed; rv++) {
266 dot_products[rv].resize(num_blocks-1);
267 for (
int coeff=1; coeff < num_blocks; coeff++) {
269 (*block_vec_poly)[coeff].Dot(*((*evecs)(rv)), &dot);
270 dot_products[rv][coeff-1] = dot;
275 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
279 i_it!=Cijk->i_end(); ++i_it) {
280 int i = epetraCijk->GRID(index(i_it));
281 sparse_kl_coeffs->sum_term(i, i, 0, norms[i]);
288 j_it != Cijk->j_end(l_it); ++j_it) {
289 int j = epetraCijk->GCID(index(j_it));
291 i_it != Cijk->i_end(j_it); ++i_it) {
292 int i = epetraCijk->GRID(index(i_it));
293 double c = value(i_it);
294 for (
int k=1; k<num_KL_computed+1; k++) {
295 double dp = dot_products[k-1][l-1];
297 if (std::abs(v) > drop_tolerance)
298 sparse_kl_coeffs->sum_term(i,
j,k,v);
303 sparse_kl_coeffs->fillComplete();
305 bool save_tensor = params->get(
"Save Sparse 3 Tensor To File",
false);
308 std::string basename = params->get(
"Sparse 3 Tensor Base Filename",
310 std::stringstream ss;
311 ss << basename <<
"_" << idx++ <<
".mm";
313 *(epetraCijk->getStochasticRowMap()), ss.str());
317 kl_blocks.resize(num_KL_computed);
318 Teuchos::RCP<Epetra_BlockMap> kl_map =
320 sg_comm->TimeDomainComm()));
323 sg_basis, kl_map, domain_base_map, range_base_map,
324 range_sg_map, sg_comm));
325 kl_ops->setCoeffPtr(0, mean);
326 for (
int rv=0; rv<num_KL_computed; rv++) {
327 if (kl_blocks[rv] == Teuchos::null)
330 for (
int i=0; i<mean->NumMyRows(); i++) {
332 mean->NumMyRowEntries(i, num_col);
333 for (
int j=0;
j<num_col;
j++)
334 (*kl_blocks[rv])[i][
j] = (*evecs)[rv][row++];
336 kl_ops->setCoeffPtr(rv+1, kl_blocks[rv]);
339 Teuchos::RCP<Stokhos::EpetraSparse3Tensor> reducedEpetraCijk =
341 sg_basis, sparse_kl_coeffs, sg_comm,
342 epetraCijk->getStochasticRowMap(), sparse_kl_coeffs,
344 reducedEpetraCijk->transformToLocal();
348 sg_comm, sg_basis, reducedEpetraCijk,
349 domain_base_map, range_base_map,
350 domain_sg_map, range_sg_map, params));
351 kl_mat_free_op->setupOperator(kl_ops);
354 if (do_error_tests) {
355 Teuchos::Array<double> point(sg_basis->dimension());
356 for (
int i=0; i<sg_basis->dimension(); i++)
358 Teuchos::Array<double> basis_vals(sg_basis->size());
359 sg_basis->evaluateBases(point, basis_vals);
363 block_vec_poly->evaluate(basis_vals,
val);
364 val_kl.Update(1.0, (*block_vec_poly)[0], 0.0);
365 Teuchos::Array< Stokhos::OrthogPolyApprox<int,double> > rvs(num_KL_computed);
366 Teuchos::Array<double> val_rvs(num_KL_computed);
367 for (
int rv=0; rv<num_KL_computed; rv++) {
368 rvs[rv].reset(sg_basis);
370 for (
int coeff=1; coeff<num_blocks; coeff++)
371 rvs[rv][coeff] = dot_products[rv][coeff-1];
372 val_rvs[rv] = rvs[rv].evaluate(point, basis_vals);
373 val_kl.Update(val_rvs[rv], *((*evecs)(rv)), 1.0);
377 val.Update(-1.0, val_kl, 1.0);
381 std::cout <<
"Infinity norm of random field difference = " << diff/nrm
385 Epetra_Vector op_input(*domain_sg_map), op_result(*range_sg_map), op_kl_result(*range_sg_map);
386 op_input.PutScalar(1.0);
388 domain_base_map, range_base_map,
389 domain_sg_map, range_sg_map, params);
391 op.
Apply(op_input, op_result);
392 this->Apply(op_input, op_kl_result);
393 op_result.NormInf(&nrm);
394 op_result.Update(-1.0, op_kl_result, 1.0);
395 op_result.NormInf(&diff);
397 std::cout <<
"Infinity norm of operator difference = " << diff/nrm
402 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
403 "Stokhos::KLReducedMatrixFreeOperator is available " <<
404 "only when configured with Anasazi support!")