From 261a1030ea24f6dfd11329d012e01f62e48e6e3e Mon Sep 17 00:00:00 2001 From: antonl Date: Tue, 24 Mar 2026 16:37:04 +0100 Subject: [PATCH] arenas for matrices and some cleanup --- TODO.txt | 2 +- out/eigvec_n1_l0.dat | 26 ++++++------ out/eigvec_n2_l2.dat | 68 +++++++++++++++--------------- src/base/base_memory.c | 1 + src/hf/file_io.c | 27 ++++++------ src/hf/hf_base.c | 75 ++++++++++----------------------- src/hf/hf_base.h | 11 ++--- src/main.c | 82 +++++++++++++++--------------------- src/os/win32/os_core_win32.c | 74 ++++++++++++++++---------------- 9 files changed, 163 insertions(+), 203 deletions(-) diff --git a/TODO.txt b/TODO.txt index d570b1f..125b0d7 100644 --- a/TODO.txt +++ b/TODO.txt @@ -3,7 +3,7 @@ //- /** * - [-] Use arena for matrices instead of malloc maybe + [DONE] Use arena for matrices instead of malloc maybe [DONE] Factor l-dependent matrix [DONE] Print eigenfunctions for plotting [DONE] Solve for more shells diff --git a/out/eigvec_n1_l0.dat b/out/eigvec_n1_l0.dat index dd3943a..3c519f1 100644 --- a/out/eigvec_n1_l0.dat +++ b/out/eigvec_n1_l0.dat @@ -19,16 +19,16 @@ 2.445164e-08 7.783642e-09 2.403077e-09 - 7.695295e-10 - 2.325414e-10 - 7.588069e-11 - 2.207375e-11 - 7.527071e-12 - 2.032901e-12 - 7.619167e-13 - 1.762993e-13 - 8.066878e-14 - 1.303573e-14 - 9.445442e-15 --2.760762e-16 - 1.847288e-15 + 7.695297e-10 + 2.325411e-10 + 7.588094e-11 + 2.207388e-11 + 7.526951e-12 + 2.032667e-12 + 7.620890e-13 + 1.763452e-13 + 8.070080e-14 + 1.321055e-14 + 9.522466e-15 +-3.487734e-16 + 1.639149e-15 diff --git a/out/eigvec_n2_l2.dat b/out/eigvec_n2_l2.dat index 38aa2e6..e63ac94 100644 --- a/out/eigvec_n2_l2.dat +++ b/out/eigvec_n2_l2.dat @@ -1,34 +1,34 @@ --2.548436e-03 --1.768047e-02 --6.801789e-02 --1.434030e-01 --2.235356e-01 --2.907379e-01 --3.357788e-01 --3.567533e-01 --3.563788e-01 --3.395486e-01 --3.116223e-01 --2.774422e-01 --2.408901e-01 --2.047894e-01 --1.709988e-01 --1.405887e-01 --1.140369e-01 --9.140819e-02 --7.250238e-02 --5.696782e-02 --4.438296e-02 --3.431126e-02 --2.633551e-02 --2.007697e-02 --1.520409e-02 --1.143438e-02 --8.531931e-03 --6.302777e-03 --4.589046e-03 --3.262862e-03 --2.220364e-03 --1.376048e-03 --5.929916e-04 --1.933374e-04 + 2.548436e-03 + 1.768047e-02 + 6.801789e-02 + 1.434030e-01 + 2.235356e-01 + 2.907379e-01 + 3.357788e-01 + 3.567533e-01 + 3.563788e-01 + 3.395486e-01 + 3.116223e-01 + 2.774422e-01 + 2.408901e-01 + 2.047894e-01 + 1.709988e-01 + 1.405887e-01 + 1.140369e-01 + 9.140819e-02 + 7.250238e-02 + 5.696782e-02 + 4.438296e-02 + 3.431126e-02 + 2.633551e-02 + 2.007697e-02 + 1.520409e-02 + 1.143438e-02 + 8.531931e-03 + 6.302777e-03 + 4.589046e-03 + 3.262862e-03 + 2.220364e-03 + 1.376048e-03 + 5.929916e-04 + 1.933374e-04 diff --git a/src/base/base_memory.c b/src/base/base_memory.c index 1482517..52b58b7 100644 --- a/src/base/base_memory.c +++ b/src/base/base_memory.c @@ -147,6 +147,7 @@ m_arena_align(Arena *arena, U64 pow2_alignment) { root_function void m_arena_release(Arena* arena) { m_release(arena, arena->capacity); + arena = 0; } root_function ArenaTemp diff --git a/src/hf/file_io.c b/src/hf/file_io.c index e090535..908ed37 100644 --- a/src/hf/file_io.c +++ b/src/hf/file_io.c @@ -1,4 +1,5 @@ - +#define DEBUG_LOG_WRITE_STRING_LIST_TO_FILE 0 +#define DEBUG_LOG_WRITE_ARRAY_BINARY 0 function void write_array_binary_F64(String8 path_to_file, F64 *values, U32 array_size) { @@ -13,6 +14,7 @@ function void write_array_binary_F64(String8 path_to_file, F64 *values, str8_list_push(scratch.arena, &list, temp); OS_file_write(scratch.arena, file_handle, 0, list, 0); +#if DEBUG_LOG_WRITE_ARRAY_BINARY String8List log_list = {0}; str8_list_push(scratch.arena, &log_list, str8_lit("Wrote binary array data to")); @@ -22,6 +24,7 @@ function void write_array_binary_F64(String8 path_to_file, F64 *values, join.post = str8_lit("\n"); String8 log_msg = str8_list_join(scratch.arena, log_list, &join); LOG(log_msg.str); +#endif scratch_release(scratch); } OS_file_close(file_handle); @@ -34,17 +37,17 @@ function void write_string_list_to_file(Arena *arena, String8 path, OS_file_open(OS_AccessFlag_Write | OS_AccessFlag_CreateNew, path); OS_file_write(arena, file_handle, 0, *list, 0); - U32 debug = 1; - if (debug) { - String8List log_list = {0}; - str8_list_push(arena, &log_list, str8_lit("Wrote array to")); - str8_list_push(arena, &log_list, path); - StringJoin join = {0}; - join.sep = str8_lit(" "); - join.post = str8_lit("\n"); - String8 log_msg = str8_list_join(arena, log_list, &join); - LOG(log_msg.str); - } +#if DEBUG_WRITE_STRING_LIST_TO_FILE + String8List log_list = {0}; + str8_list_push(arena, &log_list, str8_lit("Wrote array to")); + str8_list_push(arena, &log_list, path); + StringJoin join = {0}; + join.sep = str8_lit(" "); + join.post = str8_lit("\n"); + String8 log_msg = str8_list_join(arena, log_list, &join); + LOG(log_msg.str); +#endif + OS_file_close(file_handle); } diff --git a/src/hf/hf_base.c b/src/hf/hf_base.c index 1c2044a..2bdd342 100644 --- a/src/hf/hf_base.c +++ b/src/hf/hf_base.c @@ -2,78 +2,43 @@ //~ Globals global GaussLegendre g_gauss_legendre = {0}; +function inline U32 +mat_get_col_major_idx(U32 i, U32 j, U32 size1) { + return j * size1 + i; +} function Mat_F64 -mat_F64(U32 size1, U32 size2) +mat_F64(Arena *arena, U32 size1, U32 size2) { Mat_F64 out = {0}; + out.arena = arena; out.size1 = size1; out.size2 = size2; - out.matrix = (F64 **)malloc(size1 * sizeof(F64 *)); - U64 data_byte_size = size1 * size2 * sizeof(F64); - out.data = (F64 *)malloc(data_byte_size); - MemoryZero(out.data, data_byte_size); - - for(U32 i = 0; i < size1; i++) - { - out.matrix[i] = &out.data[i * size2]; - } - + out.data = PushArray(arena, F64, size1*size2); + MemoryZero(out.data, size1*size2*sizeof(F64)); return out; } -function Mat_F64 -mat_F64_from_data(U32 size1, U32 size2, F64 *data) -{ - Mat_F64 out = {0}; - out.size1 = size1; - out.size2 = size2; - out.matrix = (F64 **)malloc(size1 * sizeof(F64 *)); - out.data = data; - - for(U32 i = 0; i < size1; i++) - { - out.matrix[i] = &out.data[i * size2]; - } - - return out; -} - - -function void -mat_F64_free(Mat_F64 *mat) -{ - free(mat->matrix); - free(mat->data); - mat->size1 = 0; - mat->size2 = 0; - mat->matrix = 0; - mat->data = 0; -} - - function inline void mat_F64_set(Mat_F64 *mat, U32 i, U32 j, F64 val) { - U32 row = i < mat->size1 ? i : mat->size1-1; - U32 col = j < mat->size2 ? j : mat->size2-1; - mat->matrix[row][col] = val; + U32 index = mat_get_col_major_idx(i, j, mat->size1); + mat->data[index] = val; return; } function inline F64 mat_F64_get(Mat_F64 *mat, U32 i, U32 j) { - U32 row = i < mat->size1 ? i : mat->size1-1; - U32 col = j < mat->size2 ? j : mat->size2-1; - return mat->matrix[row][col]; + U32 index = mat_get_col_major_idx(i, j, mat->size1); + return mat->data[index]; } function Mat_F64 -mat_F64_copy(Mat_F64 *src) +mat_F64_copy(Arena *arena, Mat_F64 *src) { - Mat_F64 out = mat_F64(src->size1, src->size2); + Mat_F64 out = mat_F64(arena, src->size1, src->size2); U64 data_byte_size = src->size1 * src->size2 * sizeof(F64); MemoryCopy(out.data, src->data, data_byte_size); return out; @@ -82,8 +47,11 @@ mat_F64_copy(Mat_F64 *src) function void mat_F64_copy_to_dst(Mat_F64 *dst, Mat_F64 *src) { - U64 data_byte_size = src->size1 * src->size2 * sizeof(F64); - MemoryCopy(dst->data, src->data, data_byte_size); + // We assume that the dst matrix has been initialised and has + // memory for available for the src data. + // TODO(anton): Assert that dst->data != 0 or something? + U64 matsize = src->size1 * src->size2; + MemoryCopy(dst->data, src->data, matsize*sizeof(F64)); } function void @@ -225,7 +193,10 @@ mat_invert_F64(Mat_F64 *mat) } - +function F64 +hartree2eV(F64 energy_hartree) { + return 27.2114*energy_hartree; +} diff --git a/src/hf/hf_base.h b/src/hf/hf_base.h index 49acec5..14c3ccf 100644 --- a/src/hf/hf_base.h +++ b/src/hf/hf_base.h @@ -13,9 +13,9 @@ struct Z64 { typedef struct Mat_F64 Mat_F64; struct Mat_F64 { + Arena *arena; U32 size1; U32 size2; - F64 **matrix; F64 *data; }; @@ -29,12 +29,11 @@ struct GaussLegendre { //~ Base math and utility functions // Mat_F64 functions -function Mat_F64 mat_F64(U32 size1, U32 size2); -function Mat_F64 mat_F64_from_data(U32 size1, U32 size2, F64 *data); -function void mat_F64_free(Mat_F64 *mat); +function inline U32 mat_get_col_major_idx(U32 i, U32 j, U32 size1); +function Mat_F64 mat_F64(Arena* arena, U32 size1, U32 size2); function inline void mat_F64_set(Mat_F64 *mat, U32 i, U32 j, F64 val); function inline F64 mat_F64_get(Mat_F64 *mat, U32 i, U32 j); -function Mat_F64 mat_F64_copy(Mat_F64 *src); +function Mat_F64 mat_F64_copy(Arena *arena, Mat_F64 *src); function void print_mat_F64(Mat_F64 *mat); function void mat_invert_F64(Mat_F64 *mat); @@ -46,4 +45,6 @@ function void set_up_gauss_legendre_points(Arena *arena); function void print_matrix_Z64(char *desc, int m, int n, Z64 *a, int lda); function void print_matrix_F64(char *desc, int m, int n, F64 *a, int lda); +function F64 hartree2eV(F64 energy_hartree); + #endif /* HF_BASE_H */ diff --git a/src/main.c b/src/main.c index 8cbad1b..87f249c 100644 --- a/src/main.c +++ b/src/main.c @@ -1,3 +1,4 @@ +//~ //---- // Header includes #include "base/base_inc.h" @@ -31,40 +32,27 @@ struct Eigensolution_F64 { Mat_F64 left_eigenvectors; }; -typedef struct Orbital Orbital; -struct Orbital { - U32 n; - U32 l; - U32 j; - Eigensolution_F64 eigensolution; -}; - -typedef struct Atom Atom; -struct Atom { - U32 N; - Orbital *orbitals; -}; - global Arena *g_base_arena = 0; global Arena *g_filename_arena = 0; #define NUM_ANGULAR_MOMENTA 3 -global F64 angular_momenta[NUM_ANGULAR_MOMENTA] = {0.0, 1.0, 2.0}; -global F64 *temp_wavefunction_F64; +global F64 g_angular_momenta[NUM_ANGULAR_MOMENTA] = {0.0, 1.0, 2.0}; ////// //~ /* Auxiliary routine: printing a matrix */ -function void print_eigenvalues(char *desc, int n, F64 *wr, F64 *wi) { +function void print_eigenvalues(S32 l, S32 n, F64 *wr, F64 *wi) { ArenaTemp scratch = scratch_get(0, 0); - int i, j; + S32 i, j; String8 newline = str8_lit("\n"); - String8 header = str8_pushf(scratch.arena, "\n %s\n", desc); + String8 header = str8_pushf(scratch.arena, "\n ------- \n" + "Sorted eigenvalues for l = %i\n", l); LOG(header.str); // printf("\n %s \n", desc); for (j = 0; j < n; j++) { String8 outstr = - str8_pushf(scratch.arena, " (%4.5f, %4.5f)\n", wr[j], wi[j]); + str8_pushf(scratch.arena, " (%4.5f, %4.5f) Hartree, %4.5f eV \n", + wr[j], wi[j], hartree2eV(wr[j])); LOG(outstr.str); // printf(" (%6.2f,%6.2f)", a[i+j*lda].real, a[i+j*lda].imag ); } @@ -74,20 +62,13 @@ function void print_eigenvalues(char *desc, int n, F64 *wr, F64 *wi) { scratch_release(scratch); } -function void set_up_first_matrices(Mat_F64 *H, Mat_F64 *H_l, Mat_F64 *B_inv) { +function void +set_up_first_matrices(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 // 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 { ArenaTemp scratch = scratch_get(0, 0); - /* { */ - /* LOG("Setting up matrices.\n"); */ - /* String8 setuplog = str8_pushf(scratch.arena, */ - /* "Num knotpoints: %i, Bspline order k = %i, Matrix size1 = N-k-2 - * = %i, size2 = %i\n", */ - /* N, k, mat_size1, mat_size2); */ - /* LOG(setuplog.str); */ - /* } */ F64 *t = g_bspline_ctx.knotpoints; F64 Z = 1.0; @@ -97,7 +78,7 @@ function void set_up_first_matrices(Mat_F64 *H, Mat_F64 *H_l, Mat_F64 *B_inv) { for (U32 i = 0; i < H->size1; i++) { for (U32 j = 0; j < H->size2; j++) { U32 bspl_index_i = - i + 1; // The second Bspline has index 1 in our array etc. + i + 1; // The second Bspline has index 1 in our array etc. U32 bspl_index_j = j + 1; // This logic assumes 1-indexed bsplines @@ -167,7 +148,7 @@ function void set_up_first_matrices(Mat_F64 *H, Mat_F64 *H_l, Mat_F64 *B_inv) { H->size1, g_bspline_ctx.num_bsplines - 1) .str); scratch_release(scratch); - print_mat_F64(H); + //print_mat_F64(H); LOG("\n"); // print_mat_F64(B); } @@ -222,9 +203,9 @@ function void EntryPoint(void) { //- Set up grid and write to file. set_up_grid(g_base_arena); - temp_wavefunction_F64 = (F64 *)PushArray(g_base_arena, F64, g_grid.num_steps); - write_array_binary_F64(str8_lit(grid_file_path_bin), g_grid.points, + write_array_binary_F64(str8_lit(grid_file_path_bin), + g_grid.points, g_grid.num_steps); write_array_F64(str8_lit(grid_file_path), g_grid.points, g_grid.num_steps, "%13.6e\n"); @@ -242,14 +223,14 @@ function void EntryPoint(void) { U32 k = g_bspline_ctx.order; U32 mat_size1 = N - k - 2; U32 mat_size2 = mat_size1; - Mat_F64 H_base = mat_F64(mat_size1, mat_size2); - Mat_F64 H_l_base = mat_F64(mat_size1, mat_size2); - Mat_F64 H_l = mat_F64(mat_size1, mat_size2); - Mat_F64 H = mat_F64(mat_size1, mat_size2); + Mat_F64 H_base = mat_F64(g_base_arena, mat_size1, mat_size2); + Mat_F64 H_l_base = mat_F64(g_base_arena, mat_size1, mat_size2); + Mat_F64 H_l = mat_F64(g_base_arena, mat_size1, mat_size2); + Mat_F64 H = mat_F64(g_base_arena, mat_size1, mat_size2); // This will be the inverse of B, but to start with we construct B. - Mat_F64 B_inv = mat_F64(mat_size1, mat_size2); + Mat_F64 B_inv = mat_F64(g_base_arena, mat_size1, mat_size2); // A is the actual matrix for each eigenvalue problem. - Mat_F64 A = mat_F64(H.size1, H.size2); + Mat_F64 A = mat_F64(g_base_arena, H.size1, H.size2); set_up_first_matrices(&H_base, &H_l_base, &B_inv); // Our problem is Hc = EBc, but we want to solve B^-1Hc = Ec, @@ -257,10 +238,12 @@ function void EntryPoint(void) { // zgeev mat_invert_F64(&B_inv); + // This arena is used to push results from f. ex eigenvalue computations. + Arena *mkl_arena = m_make_arena(); // For each angular momentum for (U32 ang_mom_idx = 0; ang_mom_idx < NUM_ANGULAR_MOMENTA; ang_mom_idx++) { mat_F64_copy_to_dst(&H, &H_base); - F64 l = angular_momenta[ang_mom_idx]; + F64 l = g_angular_momenta[ang_mom_idx]; if (l > 1e-16) { F64 l_factor = l * (l + 1.0); U64 mat_size = H_l.size1 * H_l.size2; @@ -276,12 +259,10 @@ function void EntryPoint(void) { S32 n = A.size1; cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n, n, n, 1.0, B_inv.data, n, H.data, n, 0.0, A.data, n); - LOG("Matrix A: \n"); - print_mat_F64(&A); + //LOG("Matrix A: \n"); + //print_mat_F64(&A); } - // This arena is used to push results from f. ex eigenvalue computations. - Arena *mkl_arena = m_make_arena(); Eigensolution_F64 eigensolution = {0}; // Solve generalised eigenvalue problem { @@ -312,13 +293,12 @@ function void EntryPoint(void) { LOG("Failed to compute eigenvalues in dgeev\n"); exit(1); } - - eigensolution.right_eigenvectors = mat_F64_from_data(ldvr, size1, vr); - + + F64 *right_eigenvectors = vr; U64 *sorted_indices = PushArray(mkl_arena, U64, size1); sort_and_get_indices_F64(wr, sorted_indices, size1); sort_by_indices_F64(wi, sorted_indices, size1); - // print_eigenvalues( "Eigenvalues sorted: ", size1, wr, wi ); + print_eigenvalues((U32)l, size1, wr, wi ); U32 i = 0; F64 energy = -1000.0; @@ -334,9 +314,11 @@ function void EntryPoint(void) { // compute_wf_norm_F64(eigensolution.right_eigenvectors.matrix[energy_index], // size1, n, ang_mom_idx); + U64 eigvec_idx = mat_get_col_major_idx(0, energy_index, size1); + F64 *eigvecs = &right_eigenvectors[eigvec_idx]; write_array_F64( get_eigenvector_filename(g_filename_arena, n, ang_mom_idx), - eigensolution.right_eigenvectors.matrix[energy_index], size1, + eigvecs, size1, "%13.6e\n"); i += 1; @@ -348,5 +330,7 @@ function void EntryPoint(void) { free((void *)work); } + } + m_arena_release(mkl_arena); } diff --git a/src/os/win32/os_core_win32.c b/src/os/win32/os_core_win32.c index f295ceb..926bab5 100644 --- a/src/os/win32/os_core_win32.c +++ b/src/os/win32/os_core_win32.c @@ -174,28 +174,28 @@ OS_file_read(Arena *arena, OS_Handle file, U64 min, U64 max) { if(handle == INVALID_HANDLE_VALUE) { // TODO(anton): accumulate errors } else { - U64 bytes_to_read = AbsoluteValueU64(max - min); - U64 bytes_actually_read = 0; - result.str = PushArray(arena, U8, bytes_to_read); - result.size = 0; - U8 *ptr = result.str; - U8 *one_past_last = result.str + bytes_to_read; - - for(;;) { - U64 unread = (U64)(one_past_last - ptr); - DWORD to_read = (DWORD)(ClampTop(unread, U32Max)); - DWORD did_read = 0; - // TODO(anton): Understand WINAPI - if(!ReadFile(handle, ptr, to_read, &did_read, 0)) { - break; - } - ptr += did_read; - result.size += did_read; - if(ptr >= one_past_last) { - break; - } - + U64 bytes_to_read = AbsoluteValueU64(max - min); + U64 bytes_actually_read = 0; + result.str = PushArray(arena, U8, bytes_to_read); + result.size = 0; + U8 *ptr = result.str; + U8 *one_past_last = result.str + bytes_to_read; + + for(;;) { + U64 unread = (U64)(one_past_last - ptr); + DWORD to_read = (DWORD)(ClampTop(unread, U32Max)); + DWORD did_read = 0; + // TODO(anton): Understand WINAPI + if(!ReadFile(handle, ptr, to_read, &did_read, 0)) { + break; } + ptr += did_read; + result.size += did_read; + if(ptr >= one_past_last) { + break; + } + + } } return result; } @@ -224,24 +224,24 @@ OS_file_write(Arena *arena, OS_Handle file, U64 off, String8List data, OS_ErrorL } else for(String8Node *node = data.first; node != 0; node = node->next) { - U8 *ptr = node->string.str; - U8 *opl = ptr + node->string.size; - for(;;) + U8 *ptr = node->string.str; + U8 *opl = ptr + node->string.size; + for(;;) + { + U64 unwritten = (U64)(opl - ptr); + DWORD to_write = (DWORD)(ClampTop(unwritten, U32Max)); + DWORD did_write = 0; + // TODO(anton): understand winapi + if(!WriteFile(handle, ptr, to_write, &did_write, 0)) { - U64 unwritten = (U64)(opl - ptr); - DWORD to_write = (DWORD)(ClampTop(unwritten, U32Max)); - DWORD did_write = 0; - // TODO(anton): understand winapi - if(!WriteFile(handle, ptr, to_write, &did_write, 0)) - { - goto fail_out; - } - ptr += did_write; - if(ptr >= opl) - { - break; - } + goto fail_out; } + ptr += did_write; + if(ptr >= opl) + { + break; + } + } } fail_out:; }