anton/hf_again
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

some refactor towards fat struct

Browse Source
This commit is contained in:
antonl 2026-03-29 21:52:29 +02:00
parent 8c3f3a2947
commit 64e28327a1
4 changed files with 170 additions and 115 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

117
src/hf/bsplines_and_grid.c
View File

@ -1,16 +1,27 @@
global BSplineCtx g_bspline_ctx = {0};
global Grid g_grid = {0}; global BSplineCtx *g_bspline_ctx = 0;
global Grid *g_grid = 0;
global U32 g_debug_bspline_matrix = 0; global U32 g_debug_bspline_matrix = 0;
function void
bspline_ctx_assign(BSplineCtx *ctx) {
g_bspline_ctx = ctx;
}
function void
grid_assign(Grid *grid) {
g_grid = grid;
}
function F64 function F64
bspline_recursion(F64 x, U32 k, U32 i) bspline_recursion(F64 x, U32 k, U32 i)
{ {
F64 *t = g_bspline_ctx.knotpoints; F64 *t = g_bspline_ctx->knotpoints;
F64 tolerance = 1e-14; F64 tolerance = 1e-14;
if(k == 1) if(k == 1)
{ {
if(i == g_bspline_ctx.num_bsplines-1 && fabs(x -g_grid.end) < tolerance ) if(i == g_bspline_ctx->num_bsplines-1 && fabs(x - g_grid->end) < tolerance )
{ {
// TODO(anton): // TODO(anton):
// This is like a hack to get the last bspline to be 1 at the last point. // This is like a hack to get the last bspline to be 1 at the last point.
@ -18,7 +29,7 @@ bspline_recursion(F64 x, U32 k, U32 i)
// to unity at the last point, actually. I have to check this. // to unity at the last point, actually. I have to check this.
return 1.0; return 1.0;
} }
else if(i < g_bspline_ctx.num_bsplines && (x >= t[i] && x < t[i+1])) else if(i < g_bspline_ctx->num_bsplines && (x >= t[i] && x < t[i+1]))
{ {
return 1.0; return 1.0;
} else { } else {
@ -45,7 +56,7 @@ bspline_recursion(F64 x, U32 k, U32 i)
function F64 function F64
compute_bspline_F64(F64 x_coord, U32 index) compute_bspline_F64(F64 x_coord, U32 index)
{ {
U32 k = g_bspline_ctx.order; U32 k = g_bspline_ctx->order;
F64 out = bspline_recursion(x_coord, k, index); F64 out = bspline_recursion(x_coord, k, index);
return out; return out;
} }
@ -54,9 +65,9 @@ compute_bspline_F64(F64 x_coord, U32 index)
function F64 function F64
compute_dBspline_F64(F64 x_coord, U32 index) compute_dBspline_F64(F64 x_coord, U32 index)
{ {
U32 k = g_bspline_ctx.order; U32 k = g_bspline_ctx->order;
F64 prefac = (F64)(k - 1); F64 prefac = (F64)(k - 1);
F64 *t = g_bspline_ctx.knotpoints; F64 *t = g_bspline_ctx->knotpoints;
F64 term1_enum = bspline_recursion(x_coord, k-1, index); F64 term1_enum = bspline_recursion(x_coord, k-1, index);
F64 term2_enum = bspline_recursion(x_coord, k-1, index+1); F64 term2_enum = bspline_recursion(x_coord, k-1, index+1);
F64 term1_denom = t[index+k-1] - t[index]; F64 term1_denom = t[index+k-1] - t[index];
@ -70,17 +81,17 @@ compute_dBspline_F64(F64 x_coord, U32 index)
function void function void
set_up_grid(Arena *arena) set_up_grid(Arena *arena)
{ {
g_grid.start = GRID_START_POINT; g_grid->start = GRID_START_POINT;
g_grid.end = GRID_END_POINT; g_grid->end = GRID_END_POINT;
g_grid.num_steps = GRID_NUM_STEPS; g_grid->num_steps = GRID_NUM_STEPS;
g_grid.points = PushArray(arena, F64, g_grid.num_steps); g_grid->points = PushArray(arena, F64, g_grid->num_steps);
F64 step_size = (g_grid.end-g_grid.start)/(F64)g_grid.num_steps; F64 step_size = (g_grid->end-g_grid->start)/(F64)g_grid->num_steps;
g_grid.points[0] = g_grid.start; g_grid->points[0] = g_grid->start;
g_grid.points[g_grid.num_steps-1] = g_grid.end; g_grid->points[g_grid->num_steps-1] = g_grid->end;
for(U32 i = 1; i < g_grid.num_steps-1; i++) for(U32 i = 1; i < g_grid->num_steps-1; i++)
{ {
g_grid.points[i] = g_grid.points[i-1] + step_size; g_grid->points[i] = g_grid->points[i-1] + step_size;
} }
} }
@ -95,13 +106,13 @@ get_bspline_index_size(U32 size1, U32 i, U32 j)
function inline U32 function inline U32
get_bspline_index(U32 i, U32 j) get_bspline_index(U32 i, U32 j)
{ {
return g_bspline_ctx.num_knotpoints * j + i; return g_bspline_ctx->num_knotpoints * j + i;
} }
function inline U32 function inline U32
get_bspline_grid_index(U32 i, U32 j) get_bspline_grid_index(U32 i, U32 j)
{ {
return g_grid.num_steps * j + i; return g_grid->num_steps * j + i;
} }
@ -111,29 +122,29 @@ set_up_bspline_context(Arena* arena)
// Create knotpoint sequence. // Create knotpoint sequence.
U32 k = 4; U32 k = 4;
U32 bspl_N = 40; U32 bspl_N = 40;
g_bspline_ctx.order = k; g_bspline_ctx->order = k;
g_bspline_ctx.num_knotpoints = bspl_N; g_bspline_ctx->num_knotpoints = bspl_N;
g_bspline_ctx.num_bsplines = bspl_N-k; g_bspline_ctx->num_bsplines = bspl_N-k;
g_bspline_ctx.num_phys_points = bspl_N-(2*k)+2; // Remove k points at each end, and then add back the first and last points. g_bspline_ctx->num_phys_points = bspl_N-(2*k)+2; // Remove k points at each end, and then add back the first and last points.
g_bspline_ctx.arena = arena; g_bspline_ctx->arena = arena;
g_bspline_ctx.knotpoints = PushArray(arena, F64, g_bspline_ctx.num_knotpoints); g_bspline_ctx->knotpoints = PushArray(arena, F64, g_bspline_ctx->num_knotpoints);
// Set up physical points; // Set up physical points;
F64 delta = (g_grid.end-g_grid.start)/(g_bspline_ctx.num_phys_points-1); F64 delta = (g_grid->end-g_grid->start)/(g_bspline_ctx->num_phys_points-1);
// Set ghost points including first physical // Set ghost points including first physical
U32 phys_point_last_index = g_bspline_ctx.num_phys_points + k-1; U32 phys_point_last_index = g_bspline_ctx->num_phys_points + k-1;
for(U32 i = 0; i < k; i++) for(U32 i = 0; i < k; i++)
{ {
g_bspline_ctx.knotpoints[i] = g_grid.start; g_bspline_ctx->knotpoints[i] = g_grid->start;
} }
for(U32 i = k; i < phys_point_last_index; i++) for(U32 i = k; i < phys_point_last_index; i++)
{ {
g_bspline_ctx.knotpoints[i] = g_bspline_ctx.knotpoints[i-1] + delta; g_bspline_ctx->knotpoints[i] = g_bspline_ctx->knotpoints[i-1] + delta;
} }
// Set the last points // Set the last points
F64 last_physical = g_grid.end; F64 last_physical = g_grid->end;
for(U32 i = phys_point_last_index; i < g_bspline_ctx.num_knotpoints; i++) for(U32 i = phys_point_last_index; i < g_bspline_ctx->num_knotpoints; i++)
{ {
g_bspline_ctx.knotpoints[i] = last_physical; g_bspline_ctx->knotpoints[i] = last_physical;
} }
} }
@ -185,11 +196,11 @@ function void
set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena) set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena)
{ {
U64 num_bsplines = g_bspline_ctx.num_bsplines; U64 num_bsplines = g_bspline_ctx->num_bsplines;
U64 k = g_bspline_ctx.order; U64 k = g_bspline_ctx->order;
F64 *t = g_bspline_ctx.knotpoints; F64 *t = g_bspline_ctx->knotpoints;
U32 num_grid_points = g_grid.num_steps; U32 num_grid_points = g_grid->num_steps;
U32 num_knotpoints = g_bspline_ctx.num_knotpoints; U32 num_knotpoints = g_bspline_ctx->num_knotpoints;
// For sanity check we make the first 4 bsplines by hand. // For sanity check we make the first 4 bsplines by hand.
{ {
@ -205,7 +216,7 @@ set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena)
F64 *dBspl9 = PushArray(arena, F64, num_grid_points); F64 *dBspl9 = PushArray(arena, F64, num_grid_points);
for(U32 i = 0; i < num_grid_points; i++) for(U32 i = 0; i < num_grid_points; i++)
{ {
F64 x = g_grid.points[i]; F64 x = g_grid->points[i];
bspl0[i] = compute_bspline_F64(x, 0); bspl0[i] = compute_bspline_F64(x, 0);
bspl1[i] = compute_bspline_F64(x, 1); bspl1[i] = compute_bspline_F64(x, 1);
bspl2[i] = compute_bspline_F64(x, 2); bspl2[i] = compute_bspline_F64(x, 2);
@ -218,7 +229,7 @@ set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena)
dBspl9[i] = compute_dBspline_F64(x, 9); dBspl9[i] = compute_dBspline_F64(x, 9);
} }
F64 test = compute_bspline_F64(g_grid.points[num_grid_points-1], 9); F64 test = compute_bspline_F64(g_grid->points[num_grid_points-1], 9);
write_array_F64(str8_lit("E:\\dev\\hf_again\\out\\bspline0.dat"), bspl0, num_grid_points, "%13.6e\n"); write_array_F64(str8_lit("E:\\dev\\hf_again\\out\\bspline0.dat"), bspl0, num_grid_points, "%13.6e\n");
write_array_F64(str8_lit("E:\\dev\\hf_again\\out\\bspline1.dat"), bspl1, num_grid_points, "%13.6e\n"); write_array_F64(str8_lit("E:\\dev\\hf_again\\out\\bspline1.dat"), bspl1, num_grid_points, "%13.6e\n");
@ -235,18 +246,18 @@ set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena)
{ {
ArenaTemp scratch = scratch_get(0,0); ArenaTemp scratch = scratch_get(0,0);
g_bspline_ctx.bsplines_grid = PushArray(arena, F64, num_grid_points*num_bsplines); g_bspline_ctx->bsplines_grid = PushArray(arena, F64, num_grid_points*num_bsplines);
g_bspline_ctx.dBsplines_grid = PushArray(arena, F64, num_grid_points*num_bsplines); g_bspline_ctx->dBsplines_grid = PushArray(arena, F64, num_grid_points*num_bsplines);
for(U32 j = 0; j < num_bsplines; j++) for(U32 j = 0; j < num_bsplines; j++)
{ {
for(U32 i = 0; i < num_grid_points; i++) for(U32 i = 0; i < num_grid_points; i++)
{ {
U32 index = get_bspline_grid_index(i, j); U32 index = get_bspline_grid_index(i, j);
F64 bspline_value = compute_bspline_F64(g_grid.points[i], j); F64 bspline_value = compute_bspline_F64(g_grid->points[i], j);
F64 dBspline_value = compute_dBspline_F64(g_grid.points[i], j); F64 dBspline_value = compute_dBspline_F64(g_grid->points[i], j);
g_bspline_ctx.bsplines_grid[index] = bspline_value; g_bspline_ctx->bsplines_grid[index] = bspline_value;
g_bspline_ctx.dBsplines_grid[index] = dBspline_value; g_bspline_ctx->dBsplines_grid[index] = dBspline_value;
if(g_debug_bspline_matrix && j == 0) if(g_debug_bspline_matrix && j == 0)
{ {
String8 out = str8_pushf(scratch.arena, "%i %i \t %13.6e\n", i, index, bspline_value); String8 out = str8_pushf(scratch.arena, "%i %i \t %13.6e\n", i, index, bspline_value);
@ -258,35 +269,35 @@ set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena)
scratch_release(scratch); scratch_release(scratch);
} }
g_bspline_ctx.bsplines = PushArray(arena, F64, num_knotpoints*num_bsplines); g_bspline_ctx->bsplines = PushArray(arena, F64, num_knotpoints*num_bsplines);
g_bspline_ctx.dBsplines = PushArray(arena, F64, num_knotpoints*num_bsplines); g_bspline_ctx->dBsplines = PushArray(arena, F64, num_knotpoints*num_bsplines);
for(U32 i = 0; i < num_knotpoints; i++) for(U32 i = 0; i < num_knotpoints; i++)
{ {
for(U32 j = 0; j < num_bsplines; j++) for(U32 j = 0; j < num_bsplines; j++)
{ {
U32 index = get_bspline_index(i, j); U32 index = get_bspline_index(i, j);
g_bspline_ctx.bsplines[index] = compute_bspline_F64(t[i], j); g_bspline_ctx->bsplines[index] = compute_bspline_F64(t[i], j);
g_bspline_ctx.dBsplines[index] = compute_dBspline_F64(t[i], j); g_bspline_ctx->dBsplines[index] = compute_dBspline_F64(t[i], j);
} }
} }
write_bspline_matrix_F64(g_bspline_ctx.bsplines_grid, write_bspline_matrix_F64(g_bspline_ctx->bsplines_grid,
num_grid_points, num_grid_points,
num_bsplines, num_bsplines,
str8_lit(bspline_grid_array_file_path)); str8_lit(bspline_grid_array_file_path));
write_bspline_matrix_F64(g_bspline_ctx.dBsplines_grid, write_bspline_matrix_F64(g_bspline_ctx->dBsplines_grid,
num_grid_points, num_grid_points,
num_bsplines, num_bsplines,
str8_lit(dBspline_grid_array_file_path)); str8_lit(dBspline_grid_array_file_path));
write_bspline_matrix_F64(g_bspline_ctx.bsplines, write_bspline_matrix_F64(g_bspline_ctx->bsplines,
num_knotpoints, num_knotpoints,
num_bsplines, num_bsplines,
str8_lit(bspline_knots_array_file_path)); str8_lit(bspline_knots_array_file_path));
write_bspline_matrix_F64(g_bspline_ctx.dBsplines, write_bspline_matrix_F64(g_bspline_ctx->dBsplines,
num_knotpoints, num_knotpoints,
num_bsplines, num_bsplines,
str8_lit(dBspline_knots_array_file_path)); str8_lit(dBspline_knots_array_file_path));

3
src/hf/bsplines_and_grid.h
View File

@ -41,5 +41,6 @@ function inline U32 get_bspline_grid_index(U32 i, U32 j);
function void set_up_bspline_context(Arena* arena); function void set_up_bspline_context(Arena* arena);
function void write_bspline_matrix_F64(F64 *bsplines, U32 size1, U32 size2, String8 filename_path); function void write_bspline_matrix_F64(F64 *bsplines, U32 size1, U32 size2, String8 filename_path);
function void set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena); function void set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena);
function void bspline_ctx_assign(BSplineCtx *ctx);
function void grid_assign(Grid *grid);
#endif /* BSPLINES_AND_GRID_H */ #endif /* BSPLINES_AND_GRID_H */

3
src/hf/hf_base.c
View File

@ -1,6 +1,5 @@
//~ Globals //~ Globals
global GaussLegendre g_gauss_legendre = {0}; GaussLegendre g_gauss_legendre = {0};
function inline U32 function inline U32
mat_get_col_major_idx(U32 i, U32 j, U32 size1) { mat_get_col_major_idx(U32 i, U32 j, U32 size1) {

162
src/main.c
View File

@ -29,30 +29,75 @@
typedef struct Eigensolution_F64 Eigensolution_F64; typedef struct Eigensolution_F64 Eigensolution_F64;
struct Eigensolution_F64 { struct Eigensolution_F64 {
U32 l; U32 l;
F64 *eigenvalues_re; F64 *eigvals_re;
F64 *eigenvalues_im; F64 *eigvals_im;
Mat_F64 right_eigenvectors; Mat_F64 right_eigvecs;
Mat_F64 left_eigenvectors; Mat_F64 left_eigvecs;
};
typedef enum ElectronConfig ElectronConfig;
enum ElectronConfig {
ECFG_1s,
ECFG_2s,
ECFG_2p,
ECFG_3s,
ECFG_3p,
ECFG_3d,
ECFG_4s,
ECFG_4p,
ECFG_4d,
ECFG_4f,
ECFG_NUM_CONFIGS
};
typedef enum AngularMomenta AngularMomenta;
enum AngularMomenta {
ANGMOM_s,
ANGMOM_p,
ANGMOM_d,
ANGMOM_f,
ANGMOM_NUM_MOMENTA
}; };
#define NUM_ANGULAR_MOMENTA 3
typedef struct Atom Atom; typedef struct Atom Atom;
struct Atom { struct Atom {
Arena *arena; String8 name;
U32 Z; U32 Z;
Eigensolution_F64 eigensolutions[NUM_ANGULAR_MOMENTA]; U32 occupancy[ECFG_NUM_CONFIGS];
};
// We use a "fat struct" approach where everything just exists here
// in a single struct.
#define MAX_NUM_ANGULAR_MOMENTA 3
typedef struct Problem Problem;
struct Problem {
Arena *arena; // Just use a single arena to start with
Grid grid;
BSplineCtx bspline_ctx;
Atom atom;
F64 angular_momentum_l[ANGMOM_NUM_MOMENTA];
Eigensolution_F64 eigsols[MAX_NUM_ANGULAR_MOMENTA];
U32 num_eigsols;
Mat_F64 H;
Mat_F64 H_l;
}; };
////// //////
//~ //~
global Arena *g_base_arena = 0; function Problem problem_create() {
global Arena *g_filename_arena = 0; Problem out = {0};
global Atom g_atom = {0}; out.arena = m_make_arena();
global F64 g_angular_momenta[NUM_ANGULAR_MOMENTA] = {0.0, 1.0, 2.0}; out.atom.name = str8_lit("Hydrogen");
out.atom.Z = 1;
////// out.atom.occupancy[ECFG_1s] = 1;
//~ out.angular_momentum_l[ANGMOM_s] = 0.0;
out.angular_momentum_l[ANGMOM_p] = 1.0;
out.angular_momentum_l[ANGMOM_d] = 2.0;
out.angular_momentum_l[ANGMOM_f] = 3.0;
return out;
}
/* Auxiliary routine: printing a matrix */ /* Auxiliary routine: printing a matrix */
function void print_eigenvalues(S32 l, S32 n, F64 *wr, F64 *wi) { function void print_eigenvalues(S32 l, S32 n, F64 *wr, F64 *wi) {
@ -77,16 +122,16 @@ function void print_eigenvalues(S32 l, S32 n, F64 *wr, F64 *wi) {
} }
function void function void
set_up_first_matrices(Atom *atom, Mat_F64 *H, Mat_F64 *H_l, Mat_F64 *B_inv) { set_up_first_matrices(Problem *problem, Mat_F64 *H, Mat_F64 *H_l, Mat_F64 *B_inv) {
// We work in units hbar = 1, bohr radius a0 = 1, electron mass m_e = 1, and // We work in units hbar = 1, bohr radius a0 = 1, electron mass m_e = 1, and
// charge e = 1, and 1/(4piepsilon_0) = 1. Set up Hamiltonian: H = // charge e = 1, and 1/(4piepsilon_0) = 1. Set up Hamiltonian: H =
// -0.5*d^2/dr^2 + l(l+1)/(2r^2) - Z/r // -0.5*d^2/dr^2 + l(l+1)/(2r^2) - Z/r
{ {
ArenaTemp scratch = scratch_get(0, 0); ArenaTemp scratch = scratch_get(0, 0);
BSplineCtx *bspl_ctx = &problem->bspline_ctx;
F64 *t = g_bspline_ctx.knotpoints; F64 *t = bspl_ctx->knotpoints;
F64 Z = (F64)atom->Z; F64 Z = (F64)problem->atom.Z;
U32 k = g_bspline_ctx.order; U32 k = bspl_ctx->order;
// Skipping first bspline // Skipping first bspline
for (U32 i = 0; i < H->size1; i++) { for (U32 i = 0; i < H->size1; i++) {
@ -159,7 +204,7 @@ set_up_first_matrices(Atom *atom, Mat_F64 *H, Mat_F64 *H_l, Mat_F64 *B_inv) {
} }
LOG(str8_pushf(scratch.arena, "H.size1=N-k-2=%i, last bspline index=%i \n", LOG(str8_pushf(scratch.arena, "H.size1=N-k-2=%i, last bspline index=%i \n",
H->size1, g_bspline_ctx.num_bsplines - 1) H->size1, bspl_ctx->num_bsplines - 1)
.str); .str);
scratch_release(scratch); scratch_release(scratch);
//print_mat_F64(H); //print_mat_F64(H);
@ -174,9 +219,9 @@ function void compute_wf_norm_F64(F64 *coeffs, U64 coeff_size, U64 n, U64 l) {
// Gauss legendre integration // Gauss legendre integration
// //
F64 norm = 0.0; F64 norm = 0.0;
for (U64 i = 0; i < g_grid.num_steps - 1; i++) { for (U64 i = 0; i < g_grid->num_steps - 1; i++) {
F64 a = g_grid.points[i]; F64 a = g_grid->points[i];
F64 b = g_grid.points[i + 1]; F64 b = g_grid->points[i + 1];
F64 prefac = 0.5 * (b - a); F64 prefac = 0.5 * (b - a);
// Only integrate non-zero intervals // Only integrate non-zero intervals
@ -210,45 +255,44 @@ function void EntryPoint(void) {
OS_InitReceipt os_receipt = OS_init(); OS_InitReceipt os_receipt = OS_init();
OS_InitGfxReceipt os_gfx_receipt = OS_gfx_init(os_receipt); OS_InitGfxReceipt os_gfx_receipt = OS_gfx_init(os_receipt);
g_filename_arena = m_make_arena_reserve(Megabytes(2)); Problem problem = problem_create();
g_base_arena = m_make_arena(); LOG(str8_pushf(problem.arena, "Created Problem-struct for %s \n", problem.atom.name).str);
g_atom.Z = 1.0; // Hydrogen set_up_gauss_legendre_points(problem.arena);
g_atom.arena = m_make_arena();
set_up_gauss_legendre_points(g_base_arena);
//- Set up grid and write to file. //- Set up grid and write to file.
set_up_grid(g_base_arena); grid_assign(&problem.grid);
set_up_grid(problem.arena);
write_array_binary_F64(str8_lit(grid_file_path_bin), write_array_binary_F64(str8_lit(grid_file_path_bin),
g_grid.points, problem.grid.points,
g_grid.num_steps); problem.grid.num_steps);
write_array_F64(str8_lit(grid_file_path), g_grid.points, g_grid.num_steps, write_array_F64(str8_lit(grid_file_path), problem.grid.points,problem.grid.num_steps,
"%13.6e\n"); "%13.6e\n");
//- The BSpline context is the knotpoints and the BSpline order etc. //- The BSpline context is the knotpoints and the BSpline order etc.
set_up_bspline_context(g_base_arena); bspline_ctx_assign(&problem.bspline_ctx);
write_array_F64(str8_lit(knotpoints_file_path), g_bspline_ctx.knotpoints, set_up_bspline_context(problem.arena);
g_bspline_ctx.num_knotpoints, "%13.6e\n"); write_array_F64(str8_lit(knotpoints_file_path), problem.bspline_ctx.knotpoints,
problem.bspline_ctx.num_knotpoints, "%13.6e\n");
//- Then we generate the BSplines and save them off for reference and //- Then we generate the BSplines and save them off for reference and
// debugging. // debugging.
set_up_bsplines_at_points_and_write_matrix_F64(g_base_arena); set_up_bsplines_at_points_and_write_matrix_F64(problem.arena);
U32 N = g_bspline_ctx.num_knotpoints; U32 N = problem.bspline_ctx.num_knotpoints;
U32 k = g_bspline_ctx.order; U32 k = problem.bspline_ctx.order;
U32 mat_size1 = N - k - 2; U32 mat_size1 = N - k - 2;
U32 mat_size2 = mat_size1; U32 mat_size2 = mat_size1;
Mat_F64 H_base = mat_F64(g_base_arena, mat_size1, mat_size2); Mat_F64 H_base = mat_F64(problem.arena, mat_size1, mat_size2);
Mat_F64 H_l_base = mat_F64(g_base_arena, mat_size1, mat_size2); Mat_F64 H_l_base = mat_F64(problem.arena, mat_size1, mat_size2);
Mat_F64 H_l = mat_F64(g_base_arena, mat_size1, mat_size2); Mat_F64 H_l = mat_F64(problem.arena, mat_size1, mat_size2);
Mat_F64 H = mat_F64(g_base_arena, mat_size1, mat_size2); Mat_F64 H = mat_F64(problem.arena, mat_size1, mat_size2);
// This will be the inverse of B, but to start with we construct B. // This will be the inverse of B, but to start with we construct B.
Mat_F64 B_inv = mat_F64(g_base_arena, mat_size1, mat_size2); Mat_F64 B_inv = mat_F64(problem.arena, mat_size1, mat_size2);
// A is the actual matrix for each eigenvalue problem. // A is the actual matrix for each eigenvalue problem.
Mat_F64 A = mat_F64(g_base_arena, H.size1, H.size2); Mat_F64 A = mat_F64(problem.arena, H.size1, H.size2);
set_up_first_matrices(&g_atom, &H_base, &H_l_base, &B_inv); set_up_first_matrices(&problem, &H_base, &H_l_base, &B_inv);
// Our problem is Hc = EBc, but we want to solve B^-1Hc = Ec, // Our problem is Hc = EBc, but we want to solve B^-1Hc = Ec,
// so we invert the B matrix and compute the product A = B^-1H before calling // so we invert the B matrix and compute the product A = B^-1H before calling
@ -257,11 +301,11 @@ function void EntryPoint(void) {
// This arena is used to push results from f. ex eigenvalue computations. // This arena is used to push results from f. ex eigenvalue computations.
// For each angular momentum // For each angular momentum
for (U32 ang_mom_idx = 0; ang_mom_idx < NUM_ANGULAR_MOMENTA; ang_mom_idx++) { for (U32 ang_mom_idx = 0; ang_mom_idx < MAX_NUM_ANGULAR_MOMENTA; ang_mom_idx++) {
ArenaTemp scratch = scratch_get(0, 0); ArenaTemp scratch = scratch_get(0, 0);
mat_F64_copy_to_dst(&H, &H_base); mat_F64_copy_to_dst(&H, &H_base);
F64 l = g_angular_momenta[ang_mom_idx]; F64 l = problem.angular_momentum_l[ang_mom_idx];
Eigensolution_F64 *eigsol = &g_atom.eigensolutions[ang_mom_idx]; Eigensolution_F64 *eigsol = &problem.eigsols[ang_mom_idx];
eigsol->l = (U32)l; eigsol->l = (U32)l;
if (l > 1e-16) { if (l > 1e-16) {
@ -295,14 +339,14 @@ function void EntryPoint(void) {
F64 wkopt; F64 wkopt;
F64 *work; F64 *work;
eigsol->eigenvalues_re = PushArray(g_atom.arena, F64, size1); eigsol->eigvals_re= PushArray(problem.arena, F64, size1);
F64 *wr = eigsol->eigenvalues_re; F64 *wr = eigsol->eigvals_re;
eigsol->eigenvalues_im = PushArray(g_atom.arena, F64, size1); eigsol->eigvals_im = PushArray(problem.arena, F64, size1);
F64 *wi = eigsol->eigenvalues_im; F64 *wi = eigsol->eigvals_im;
eigsol->left_eigenvectors = mat_F64(g_atom.arena, ldvl, size1); eigsol->left_eigvecs = mat_F64(problem.arena, ldvl, size1);
F64 *vl = eigsol->left_eigenvectors.data; F64 *vl = eigsol->left_eigvecs.data;
eigsol->right_eigenvectors = mat_F64(g_atom.arena, size1, ldvr); eigsol->right_eigvecs = mat_F64(problem.arena, size1, ldvr);
F64 *vr = eigsol->right_eigenvectors.data; F64 *vr = eigsol->right_eigvecs.data;
lwork = -1; lwork = -1;
F64 *a = A.data; F64 *a = A.data;
@ -339,9 +383,9 @@ function void EntryPoint(void) {
// size1, n, ang_mom_idx); // size1, n, ang_mom_idx);
U64 eigvec_idx = mat_get_col_major_idx(0, energy_index, size1); U64 eigvec_idx = mat_get_col_major_idx(0, energy_index, size1);
F64 *eigvecs = &eigsol->right_eigenvectors.data[eigvec_idx]; F64 *eigvecs = &eigsol->right_eigvecs.data[eigvec_idx];
write_array_F64( write_array_F64(
get_eigenvector_filename(g_filename_arena, n, ang_mom_idx), get_eigenvector_filename(scratch.arena, n, ang_mom_idx),
eigvecs, size1, eigvecs, size1,
"%13.6e\n"); "%13.6e\n");