From 588f1f1e4689f7ff90f4be4c8635c73c03b964c7 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 27 Mar 2019 16:54:02 -0600 Subject: [PATCH] Make all but one of the changes requested by Timo. They are all no-ops -- the output is the exact same as before. --- MCMC-Laplace/mcmc-laplace.cc | 55 +++++++++++++++++------------------- 1 file changed, 26 insertions(+), 29 deletions(-) diff --git a/MCMC-Laplace/mcmc-laplace.cc b/MCMC-Laplace/mcmc-laplace.cc index 179fe9b..9d59373 100644 --- a/MCMC-Laplace/mcmc-laplace.cc +++ b/MCMC-Laplace/mcmc-laplace.cc @@ -91,24 +91,24 @@ namespace ForwardSimulator void solve(); void output_results(const Vector &coefficients) const; - Triangulation triangulation; - FE_Q fe; - DoFHandler dof_handler; + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; - FullMatrix cell_matrix; - Vector cell_rhs; - std::map boundary_values; + FullMatrix cell_matrix; + Vector cell_rhs; + std::map boundary_values; - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; - Vector solution; - Vector system_rhs; + Vector solution; + Vector system_rhs; - std::vector> measurement_points; + std::vector> measurement_points; - SparsityPattern measurement_sparsity; - SparseMatrix measurement_matrix; + SparsityPattern measurement_sparsity; + SparseMatrix measurement_matrix; TimerOutput timer; unsigned int nth_evaluation; @@ -177,11 +177,6 @@ namespace ForwardSimulator const unsigned int n_points_per_direction = 13; const double dx = 1. / (n_points_per_direction + 1); - Vector weights(dof_handler.n_dofs()); - FullMatrix full_measurement_matrix(n_points_per_direction * - n_points_per_direction, - dof_handler.n_dofs()); - for (unsigned int x = 1; x <= n_points_per_direction; ++x) for (unsigned int y = 1; y <= n_points_per_direction; ++y) measurement_points.emplace_back(x * dx, y * dx); @@ -191,6 +186,11 @@ namespace ForwardSimulator // which entries are nonzero. Later, the `copy_from()` function // calls build a sparsity pattern and a sparse matrix from // the dense matrix. + Vector weights(dof_handler.n_dofs()); + FullMatrix full_measurement_matrix(n_points_per_direction * + n_points_per_direction, + dof_handler.n_dofs()); + for (unsigned int index = 0; index < measurement_points.size(); ++index) { VectorTools::create_point_source_vector(dof_handler, @@ -271,27 +271,24 @@ namespace ForwardSimulator const unsigned int dofs_per_cell = fe.dofs_per_cell; - FullMatrix this_cell_matrix(dofs_per_cell, dofs_per_cell); - std::vector local_dof_indices(dofs_per_cell); for (const auto &cell : dof_handler.active_cell_iterators()) { - this_cell_matrix = cell_matrix; - this_cell_matrix *= coefficients(cell->user_index()); - - cell->get_dof_indices(local_dof_indices); + const double coefficient = coefficients(cell->user_index()); + + cell->get_dof_indices(local_dof_indices); for (unsigned int i = 0; i < dofs_per_cell; ++i) { for (unsigned int j = 0; j < dofs_per_cell; ++j) system_matrix.add(local_dof_indices[i], local_dof_indices[j], - this_cell_matrix(i, j)); + coefficient * cell_matrix(i, j)); system_rhs(local_dof_indices[i]) += cell_rhs(i); } } - + MatrixTools::apply_boundary_values(boundary_values, system_matrix, solution, @@ -353,8 +350,8 @@ namespace ForwardSimulator // and creating output for each one of them is surely more data // than you ever want to see! // - // At the end of the function, we output some timing information for - // the computations once for every 10,000 samples. + // At the end of the function, we output some timing information + // every 10,000 samples. template Vector PoissonSolver::evaluate(const Vector &coefficients) @@ -397,7 +394,7 @@ namespace ForwardSimulator // of the `Gaussian` class and standard deviation $\sigma$. // // For reasons of numerical accuracy, it is useful to not return the -// actual likelihood, but it's logarithm. This is because these +// actual likelihood, but its logarithm. This is because these // values can be very small, occasionally on the order of $e^{-100}$, // for which it becomes very difficult to compute accurate // values. -- 2.39.5