arenas for matrices and some cleanup

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

out/eigvec_n1_l0.dat

@@ -19,16 +19,16 @@
  2.445164e-08
  7.783642e-09
  2.403077e-09
- 7.695295e-10
+ 7.695297e-10
- 2.325414e-10
+ 2.325411e-10
- 7.588069e-11
+ 7.588094e-11
- 2.207375e-11
+ 2.207388e-11
- 7.527071e-12
+ 7.526951e-12
- 2.032901e-12
+ 2.032667e-12
- 7.619167e-13
+ 7.620890e-13
- 1.762993e-13
+ 1.763452e-13
- 8.066878e-14
+ 8.070080e-14
- 1.303573e-14
+ 1.321055e-14
- 9.445442e-15
+ 9.522466e-15
--2.760762e-16
+-3.487734e-16
- 1.847288e-15
+ 1.639149e-15

out/eigvec_n2_l2.dat

@@ -1,34 +1,34 @@
--2.548436e-03
+ 2.548436e-03
--1.768047e-02
+ 1.768047e-02
--6.801789e-02
+ 6.801789e-02
--1.434030e-01
+ 1.434030e-01
--2.235356e-01
+ 2.235356e-01
--2.907379e-01
+ 2.907379e-01
--3.357788e-01
+ 3.357788e-01
--3.567533e-01
+ 3.567533e-01
--3.563788e-01
+ 3.563788e-01
--3.395486e-01
+ 3.395486e-01
--3.116223e-01
+ 3.116223e-01
--2.774422e-01
+ 2.774422e-01
--2.408901e-01
+ 2.408901e-01
--2.047894e-01
+ 2.047894e-01
--1.709988e-01
+ 1.709988e-01
--1.405887e-01
+ 1.405887e-01
--1.140369e-01
+ 1.140369e-01
--9.140819e-02
+ 9.140819e-02
--7.250238e-02
+ 7.250238e-02
--5.696782e-02
+ 5.696782e-02
--4.438296e-02
+ 4.438296e-02
--3.431126e-02
+ 3.431126e-02
--2.633551e-02
+ 2.633551e-02
--2.007697e-02
+ 2.007697e-02
--1.520409e-02
+ 1.520409e-02
--1.143438e-02
+ 1.143438e-02
--8.531931e-03
+ 8.531931e-03
--6.302777e-03
+ 6.302777e-03
--4.589046e-03
+ 4.589046e-03
--3.262862e-03
+ 3.262862e-03
--2.220364e-03
+ 2.220364e-03
--1.376048e-03
+ 1.376048e-03
--5.929916e-04
+ 5.929916e-04
--1.933374e-04
+ 1.933374e-04

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

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_WRITE_STRING_LIST_TO_FILE
-  if (debug) {
+  String8List log_list = {0};
-    String8List log_list = {0};
+  str8_list_push(arena, &log_list, str8_lit("Wrote array to"));
-    str8_list_push(arena, &log_list, str8_lit("Wrote array to"));
+  str8_list_push(arena, &log_list, path);
-    str8_list_push(arena, &log_list, path);
+  StringJoin join = {0};
-    StringJoin join = {0};
+  join.sep = str8_lit(" ");
-    join.sep = str8_lit(" ");
+  join.post = str8_lit("\n");
-    join.post = str8_lit("\n");
+  String8 log_msg = str8_list_join(arena, log_list, &join);
-    String8 log_msg = str8_list_join(arena, log_list, &join);
+  LOG(log_msg.str);
-    LOG(log_msg.str);
+#endif
-  }
+
   OS_file_close(file_handle);
 }
 

src/hf/hf_base.c

@@ -2,78 +2,43 @@
 //~ Globals
 global GaussLegendre g_gauss_legendre = {0};
 
-
+function inline U32
-function Mat_F64 
+mat_get_col_major_idx(U32 i, U32 j, U32 size1) {
-mat_F64(U32 size1, U32 size2)
+  return j * size1 + i;
-{
-  Mat_F64 out = {0};
-  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];
-  }
-
-  return out;
 }
 
 function Mat_F64 
-mat_F64_from_data(U32 size1, U32 size2, F64 *data)
+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 *));
+  out.data = PushArray(arena, F64, size1*size2);
-  out.data = data;
+  MemoryZero(out.data, size1*size2*sizeof(F64));
-
-  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 index = mat_get_col_major_idx(i, j, mat->size1);
-  U32 col = j < mat->size2 ? j : mat->size2-1;
+  mat->data[index] = val;
-  mat->matrix[row][col] = 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 index = mat_get_col_major_idx(i, j, mat->size1);
-  U32 col = j < mat->size2 ? j : mat->size2-1;
+  return mat->data[index];
-  return mat->matrix[row][col];
 }
 
 
 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);
+  // We assume that the dst matrix has been initialised and has
-  MemoryCopy(dst->data, src->data, data_byte_size);
+  // 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;
+}
 
 
 

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 inline U32 mat_get_col_major_idx(U32 i, U32 j, U32 size1);
-function Mat_F64 mat_F64_from_data(U32 size1, U32 size2, F64 *data);
+function Mat_F64 mat_F64(Arena* arena, U32 size1, U32 size2);
-function void mat_F64_free(Mat_F64 *mat);
 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 */

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 g_angular_momenta[NUM_ANGULAR_MOMENTA] = {0.0, 1.0, 2.0};
-global F64 *temp_wavefunction_F64;
 
 //////
 //~
 
 /* 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_base = mat_F64(g_base_arena, mat_size1, mat_size2);
-  Mat_F64 H_l_base = mat_F64(mat_size1, mat_size2);
+  Mat_F64 H_l_base = mat_F64(g_base_arena, mat_size1, mat_size2);
-  Mat_F64 H_l = mat_F64(mat_size1, mat_size2);
+  Mat_F64 H_l = mat_F64(g_base_arena, mat_size1, mat_size2);
-  Mat_F64 H = mat_F64(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");
+      //LOG("Matrix A: \n");
-      print_mat_F64(&A);
+      //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
     {
@@ -313,12 +294,11 @@ function void EntryPoint(void) {
         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);
 }

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_to_read = AbsoluteValueU64(max - min);
-        U64 bytes_actually_read = 0;
+      U64 bytes_actually_read = 0;
-        result.str = PushArray(arena, U8, bytes_to_read);
+      result.str = PushArray(arena, U8, bytes_to_read);
-        result.size = 0;
+      result.size = 0;
-        U8 *ptr = result.str;
+      U8 *ptr = result.str;
-        U8 *one_past_last = result.str + bytes_to_read;
+      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;
-            }
       
+      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 *ptr = node->string.str;
-        U8 *opl = ptr + node->string.size;
+      U8 *opl = ptr + node->string.size;
-        for(;;)
+      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);
+            goto fail_out;
-            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;
-            }
         }
+        ptr += did_write;
+        if(ptr >= opl)
+        {
+            break;
+        }
+      }
     }
     fail_out:;
 }