FML
Typedefs | Functions
FML::HESSIAN Namespace Reference

This namespace deals with computing the Hessian of a grid and related algorithms that rely on this. More...

Typedefs

template<int N>
using FFTWGrid = FML::GRID::FFTWGrid< N >
 

Functions

template<int N>
void ComputeHessianWithFT (const FFTWGrid< N > &f_real, std::vector< FFTWGrid< N > > &hessian_real, double norm=1.0, bool hessian_of_potential_of_f=false)
 Computes the Hessian matrix of a grid [norm * f] via Fourier transforms. More...
 
template<int N>
void SymmetricTensorEigensystem (std::vector< FFTWGrid< N > > &tensor_real, std::vector< FFTWGrid< N > > &eigenvalues, std::vector< FFTWGrid< N > > &eigenvectors, bool compute_eigenvectors=false)
 For each point in the grid compute eigenvectors and eigenvalues of the tensor \( H_{ij} \) where tensor_real contains the \( N(N-1)/2 \) grids [ 00,01,02,..,11,12,...,NN ]. More...
 

Detailed Description

This namespace deals with computing the Hessian of a grid and related algorithms that rely on this.

Typedef Documentation

◆ FFTWGrid

template<int N>
using FML::HESSIAN::FFTWGrid = typedef FML::GRID::FFTWGrid<N>

Definition at line 17 of file Hessian.h.

Function Documentation

◆ ComputeHessianWithFT()

template<int N>
void FML::HESSIAN::ComputeHessianWithFT ( const FFTWGrid< N > &  f_real,
std::vector< FFTWGrid< N > > &  hessian_real,
double  norm = 1.0,
bool  hessian_of_potential_of_f = false 
)

Computes the Hessian matrix of a grid [norm * f] via Fourier transforms.

If hessian_of_potential_of_f is true then we compute the Hessian \( \phi_{ij} \) where \( \nabla^2 \phi = norm * f_{\rm real} \) Since \( f_{ij} = f_{ji} \) we only compute the elements for \( j \geq i \) and they are stored in the order fxx fxy ... fyy fyz ... etc. in hessian_real In 2D: [fxx fxy fyy] In 3D: [fxx fxy fxz fyy fyz fzz]

Template Parameters
NThe dimension we are working in
Parameters
[in]f_realThe grid we are to compute the hessian of
[out]hessian_realThe hessian of the grid (or its potential, see below)
[in]normA number to scale the grid by if needed (default is 1.0)
[in]hessian_of_potential_of_fCompute the hessian of the potential of the grid (default is false)

Definition at line 37 of file Hessian.h.

40 {
41
42 assert_mpi(f_real.get_nmesh() > 0, "[ComputeHessianWithFT] f_real grid is not allocated\n");
43
44 // Helper function to go from f(k) -> DiDj f
45 // Assuing we have f(k) in grid
46 auto ComputeSecondDerivative = [&](FFTWGrid<N> & grid, int i1, int i2) {
47 if (FML::ThisTask == 0)
48 std::cout << "[ComputeHessianWithFT::ComputeSecondDerivative] Computing phi_" << i1 << "," << i2
49 << "\n";
50
51 auto Local_nx = grid.get_local_nx();
52 auto Local_x_start = grid.get_local_x_start();
53#ifdef USE_OMP
54#pragma omp parallel for
55#endif
56 for (int islice = 0; islice < Local_nx; islice++) {
57 [[maybe_unused]] double kmag2;
58 [[maybe_unused]] std::array<double, N> kvec;
59 for (auto && fourier_index : grid.get_fourier_range(islice, islice + 1)) {
60 if(Local_x_start == 0 and fourier_index == 0)
61 continue; // DC mode (k=0)
62
63 grid.get_fourier_wavevector_and_norm2_by_index(fourier_index, kvec, kmag2);
64
65 // From f(k) -> -ika ikb f(k) / k^2 = (ka kb / k^2) f(k)
66 auto value = grid.get_fourier_from_index(fourier_index);
67 double factor = -norm * kvec[i1] * kvec[i2];
68 if (hessian_of_potential_of_f)
69 factor *= -1.0 / kmag2;
70 value *= factor;
71
72 grid.set_fourier_from_index(fourier_index, value);
73 }
74 }
75
76 // Deal with the DC mode
77 if (FML::ThisTask == 0)
78 grid.set_fourier_from_index(0, 0.0);
79
80 // Back to real space
81 grid.fftw_c2r();
82 };
83
84 // Take a copy and Fourier transform it
85 FFTWGrid<N> f_fourier = f_real;
86 f_fourier.fftw_r2c();
87
88 // Allocate grids
89 hessian_real.resize(N * (N - 1));
90
91 // Compute hessian matrix
92 int count = 0;
93 for (int idim = 0; idim < N; idim++) {
94 for (int idim2 = idim; idim2 < N; idim2++) {
95 hessian_real[count] = f_fourier;
96 ComputeSecondDerivative(hessian_real[count], idim, idim2);
97 count++;
98 }
99 }
100 }
Class for holding grids and performing real-to-complex and complex-to-real FFTs using FFTW with MPI.
Definition: FFTWGrid.h:88
int get_nmesh() const
Number of grid-cells per dimension.
Definition: FFTWGrid.h:1054
void fftw_r2c()
Perform real-to-complex fourier transform.
Definition: FFTWGrid.h:881
void grid(bool flag)
int ThisTask
The MPI task number.
Definition: Global.cpp:33

◆ SymmetricTensorEigensystem()

template<int N>
void FML::HESSIAN::SymmetricTensorEigensystem ( std::vector< FFTWGrid< N > > &  tensor_real,
std::vector< FFTWGrid< N > > &  eigenvalues,
std::vector< FFTWGrid< N > > &  eigenvectors,
bool  compute_eigenvectors = false 
)

For each point in the grid compute eigenvectors and eigenvalues of the tensor \( H_{ij} \) where tensor_real contains the \( N(N-1)/2 \) grids [ 00,01,02,..,11,12,...,NN ].

Eigenvalues are ordered in descending order

Eigenvectors are stored in row major order in the grid vector

This allocates N grids if compute_eigenvectors = false and N(N+1) grids otherwise

Definition at line 115 of file Hessian.h.

118 {
119
120 assert_mpi(tensor_real.size() > 0, "[SymmetricTensorEigensystem] tensor_real is not allocated\n");
121 assert_mpi(tensor_real[0].get_nmesh() > 0,
122 "[SymmetricTensorEigensystem] tensor_real[0] is not allocated\n");
123 for (size_t i = 1; i < tensor_real.size(); i++)
124 assert_mpi(tensor_real[i].get_nmesh() == tensor_real[i - 1].get_nmesh(),
125 "[SymmetricTensorEigensystem] tensor_real[i] is not allocated\n");
126
127 // N eigenvalues
128 eigenvalues.resize(N);
129 for (int idim = 0; idim < N; idim++)
130 eigenvalues[idim] = tensor_real[0];
131
132 // N eigenvectors with N components
133 // We store the components in the same (row major) order as GSL uses
134 if (compute_eigenvectors) {
135 eigenvectors.resize(N * N);
136 for (int i = 0; i < N * N; i++)
137 eigenvectors[i] = tensor_real[0];
138 }
139
140 // Set up the GSL stuff we need
141 gsl_matrix * matrix = gsl_matrix_alloc(N, N);
142 gsl_matrix * evec = gsl_matrix_alloc(N, N);
143 gsl_vector * eval = gsl_vector_alloc(N);
144 gsl_eigen_symm_workspace * workspace = gsl_eigen_symm_alloc(N);
145 gsl_eigen_symmv_workspace * workspacev = gsl_eigen_symmv_alloc(N);
146
147 // Solves the full eigensystem
148 auto SolveEigensystem = [&](gsl_matrix * _matrix,
149 gsl_vector * _eval,
150 gsl_matrix * _evec,
151 gsl_eigen_symmv_workspace * _workspace) {
152 // Compute eigenvalues and eigenvectors
153 gsl_eigen_symmv(_matrix, _eval, _evec, _workspace);
154 // Sort in descending order
155 gsl_eigen_symmv_sort(_eval, _evec, GSL_EIGEN_SORT_VAL_DESC);
156 };
157
158 // Solves just for eigenvalues
159 auto SolveEigenvalues =
160 [&](gsl_matrix * _matrix, gsl_vector * _eval, gsl_eigen_symm_workspace * _workspace) {
161 // Compute eigenvalues and eigenvectors
162 gsl_eigen_symm(_matrix, _eval, _workspace);
163 // Order the eigenvalues in descending order
164 std::sort(_eval->data, _eval->data + _matrix->size1, std::greater<double>());
165 };
166
167 // Loop over all cells
168 auto Local_nx = tensor_real[0].get_local_nx();
169#ifdef USE_OMP
170#pragma omp parallel for
171#endif
172 for (int islice = 0; islice < Local_nx; islice++) {
173 for (auto && real_index : tensor_real[0].get_real_range(islice,islice+1)) {
174
175 // Set the matrix
176 int count = 0;
177 for (int idim = 0; idim < N; idim++) {
178 auto value = tensor_real[count].get_real_from_index(real_index);
179 gsl_matrix_set(matrix, idim, idim, value);
180 count++;
181 for (int idim2 = idim + 1; idim2 < N; idim2++) {
182 value = tensor_real[count].get_real_from_index(real_index);
183 gsl_matrix_set(matrix, idim, idim2, value);
184 gsl_matrix_set(matrix, idim2, idim, value);
185 count++;
186 }
187 }
188
189 // Compute eigenvectors+eigenvalues or just eigenvalues
190 // In the latter case we sort the eigenvalues
191 if (compute_eigenvectors) {
192 SolveEigensystem(matrix, eval, evec, workspacev);
193
194 // Set eigenvectors
195 for (int i = 0; i < N * N; i++) {
196 eigenvectors[i].set_real_from_index(real_index, evec->data[i]);
197 // For column major order: gsl_matrix_get(evec, i / N, i % N);
198 }
199
200 } else {
201 SolveEigenvalues(matrix, eval, workspace);
202 }
203
204 // Store the eigenvalues
205 for (int idim = 0; idim < N; idim++)
206 eigenvalues[idim].set_real_from_index(real_index, eval->data[idim]);
207 }
208 }
209
210 // Free up GSL allocations
211 gsl_matrix_free(matrix);
212 gsl_matrix_free(evec);
213 gsl_vector_free(eval);
214 gsl_eigen_symm_free(workspace);
215 gsl_eigen_symmv_free(workspacev);
216 }