TODO.txt

@@ -3,7 +3,7 @@
 //-
 /**
  *
- [DONE] Use arena for matrices instead of malloc maybe
+ [-] 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.695297e-10
+ 7.695295e-10
- 2.325411e-10
+ 2.325414e-10
- 7.588094e-11
+ 7.588069e-11
- 2.207388e-11
+ 2.207375e-11
- 7.526951e-12
+ 7.527071e-12
- 2.032667e-12
+ 2.032901e-12
- 7.620890e-13
+ 7.619167e-13
- 1.763452e-13
+ 1.762993e-13
- 8.070080e-14
+ 8.066878e-14
- 1.321055e-14
+ 1.303573e-14
- 9.522466e-15
+ 9.445442e-15
--3.487734e-16
+-2.760762e-16
- 1.639149e-15
+ 1.847288e-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_inc.c

@@ -1,5 +1,4 @@
 #include "base_math.c"
 #include "base_memory.c"
 #include "base_strings.c"
-#include "base_thread_context.c"
+#include "base_thread_context.c"
-#include "base_sort.c"

src/base/base_inc.h

@@ -7,6 +7,5 @@
 #include "base_memory.h"
 #include "base_strings.h"
 #include "base_thread_context.h"
-#include "base_sort.h"
 
 #endif //BASE_H

src/base/base_memory.c

@@ -147,7 +147,6 @@ 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/base/base_sort.c

@@ -1,60 +0,0 @@
-
-U64 qsort_partition_F64(SortPair_F64 *array, U64 size, U64 low, U64 high) {
-  F64 pivot = array[high].value;
-  U64 i = low - 1;
-
-  SortPair_F64 temp = {0};
-  for (U64 j = low; j < high; j++) {
-    if (array[j].value <= pivot) {
-      i += 1;
-      temp = array[i];
-      array[i] = array[j];
-      array[j] = temp;
-    }
-  }
-
-  U64 final_pivot_pos = i + 1;
-  temp = array[final_pivot_pos];
-  array[final_pivot_pos] = array[high];
-  array[high] = temp;
-
-  return final_pivot_pos;
-}
-
-function void qsort_F64(SortPair_F64 *array, U64 size, U64 low, U64 high) {
-  if (low < high) {
-    U64 pivot_index = qsort_partition_F64(array, size, low, high);
-    qsort_F64(array, size, low, pivot_index - 1);
-    qsort_F64(array, size, pivot_index + 1, high);
-  }
-}
-
-function void sort_and_get_indices_F64(F64 *array, U64 *indices, U64 size) {
-  SortPair_F64 *pairs = malloc(size * sizeof(SortPair_F64));
-  for (U64 i = 0; i < size; i++) {
-    pairs[i].value = array[i];
-    pairs[i].original_index = i;
-  }
-
-  qsort_F64(pairs, size, 0, size - 1);
-
-  for (U32 i = 0; i < size; i++) {
-    array[i] = pairs[i].value;
-    indices[i] = pairs[i].original_index;
-  }
-
-  free(pairs);
-}
-
-function void sort_by_indices_F64(F64 *array, U64 *indices, U64 size) {
-
-  F64 *temp = malloc(size * sizeof(F64));
-  for (U64 i = 0; i < size; i++) {
-    U64 original_index = indices[i];
-    temp[i] = array[original_index];
-  }
-  for (U64 i = 0; i < size; i++) {
-    array[i] = temp[i];
-  }
-  free(temp);
-}

src/base/base_sort.h

@@ -1,15 +0,0 @@
-#ifndef BASE_SORT_H
-#define BASE_SORT_H
-
-typedef struct SortPair_F64 SortPair_F64;
-struct SortPair_F64 {
-  F64 value;
-  U64 original_index;
-};
-
-function U64 qsort_partition_F64(SortPair_F64 *array, U64 size, U64 low, U64 high);
-function void qsort_F64(SortPair_F64 *array, U64 size, U64 low, U64 high);
-function void sort_and_get_indices_F64(F64 *array, U64 *indices, U64 size);
-function void sort_by_indices_F64(F64 *array, U64 *indices, U64 size);
-
-#endif //BASE_SORT_H

src/hf/file_io.c

@@ -1,5 +1,4 @@
-#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) {
@@ -14,7 +13,6 @@ 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"));
@@ -24,7 +22,6 @@ 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);
@@ -37,17 +34,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);
 
-#if DEBUG_WRITE_STRING_LIST_TO_FILE
+  U32 debug = 1;
-  String8List log_list = {0};
+  if (debug) {
-  str8_list_push(arena, &log_list, str8_lit("Wrote array to"));
+    String8List log_list = {0};
-  str8_list_push(arena, &log_list, path);
+    str8_list_push(arena, &log_list, str8_lit("Wrote array to"));
-  StringJoin join = {0};
+    str8_list_push(arena, &log_list, path);
-  join.sep = str8_lit(" ");
+    StringJoin join = {0};
-  join.post = str8_lit("\n");
+    join.sep = str8_lit(" ");
-  String8 log_msg = str8_list_join(arena, log_list, &join);
+    join.post = str8_lit("\n");
-  LOG(log_msg.str);
+    String8 log_msg = str8_list_join(arena, log_list, &join);
-#endif
+    LOG(log_msg.str);
-
+  }
   OS_file_close(file_handle);
 }
 

src/hf/hf_base.c

@@ -2,43 +2,78 @@
 //~ 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(Arena *arena, U32 size1, U32 size2)
+mat_F64(U32 size1, U32 size2)
 {
   Mat_F64 out = {0};
-  out.arena = arena;
   out.size1 = size1;
   out.size2 = size2;
-  out.data = PushArray(arena, F64, size1*size2);
+  out.matrix = (F64 **)malloc(size1 * sizeof(F64 *));
-  MemoryZero(out.data, size1*size2*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 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 index = mat_get_col_major_idx(i, j, mat->size1);
+  U32 row = i < mat->size1 ? i : mat->size1-1;
-  mat->data[index] = val;
+  U32 col = j < mat->size2 ? j : mat->size2-1;
+  mat->matrix[row][col] = val;
   return;
 }
 
 function inline F64
 mat_F64_get(Mat_F64 *mat, U32 i, U32 j)
 {
-  U32 index = mat_get_col_major_idx(i, j, mat->size1);
+  U32 row = i < mat->size1 ? i : mat->size1-1;
-  return mat->data[index];
+  U32 col = j < mat->size2 ? j : mat->size2-1;
+  return mat->matrix[row][col];
 }
 
 
 function Mat_F64
-mat_F64_copy(Arena *arena, Mat_F64 *src) 
+mat_F64_copy(Mat_F64 *src) 
 {
-  Mat_F64 out = mat_F64(arena, src->size1, src->size2);
+  Mat_F64 out = mat_F64(src->size1, src->size2);
   U64 data_byte_size = src->size1 * src->size2 * sizeof(F64);
   MemoryCopy(out.data, src->data, data_byte_size);
   return out;
@@ -47,11 +82,8 @@ mat_F64_copy(Arena *arena, Mat_F64 *src)
 function void
 mat_F64_copy_to_dst(Mat_F64 *dst, Mat_F64 *src) 
 {
-  // We assume that the dst matrix has been initialised and has
+  U64 data_byte_size = src->size1 * src->size2 * sizeof(F64);
-  // memory for available for the src data.
+  MemoryCopy(dst->data, src->data, data_byte_size);
-  // 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
@@ -193,10 +225,7 @@ 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,11 +29,12 @@ struct GaussLegendre {
 //~ Base math and utility functions
 
 // Mat_F64 functions
-function inline U32 mat_get_col_major_idx(U32 i, U32 j, U32 size1);
+function Mat_F64 mat_F64(U32 size1, U32 size2);
-function Mat_F64 mat_F64(Arena* arena, 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 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(Arena *arena, Mat_F64 *src);
+function Mat_F64 mat_F64_copy(Mat_F64 *src);
 function void print_mat_F64(Mat_F64 *mat);
 
 function void mat_invert_F64(Mat_F64 *mat);
@@ -45,6 +46,4 @@ 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,5 +1,4 @@
-//~
+// ---
-//----
 // Header includes
 #include "base/base_inc.h"
 #include "os/os_inc.h"
@@ -8,7 +7,7 @@
 #include "hf/file_io.h"
 #include "hf/hf_base.h"
 
-//----
+// ---
 // .C includes
 #include "base/base_inc.c"
 #include "os/os_entry_point.c"
@@ -21,9 +20,6 @@
 // TODO make this a separate module that can be compiled instead
 // #include "hf/tests.c"
 
-//////
-//~ 
-
 typedef struct Eigensolution_F64 Eigensolution_F64;
 struct Eigensolution_F64 {
   F64 *eigenvalues_re;
@@ -32,27 +28,103 @@ 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;
+};
+
+typedef struct SortPair_F64 SortPair_F64;
+struct SortPair_F64 {
+  F64 value;
+  U64 original_index;
+};
+
 global Arena *g_base_arena = 0;
 global Arena *g_filename_arena = 0;
 #define NUM_ANGULAR_MOMENTA 3
-global F64 g_angular_momenta[NUM_ANGULAR_MOMENTA] = {0.0, 1.0, 2.0};
+global F64 angular_momenta[NUM_ANGULAR_MOMENTA] = {0.0, 1.0, 2.0};
+global F64 *temp_wavefunction_F64;
 
-//////
+U64 qsort_partition_F64(SortPair_F64 *array, U64 size, U64 low, U64 high) {
-//~
+  F64 pivot = array[high].value;
+  U64 i = low - 1;
+
+  SortPair_F64 temp = {0};
+  for (U64 j = low; j < high; j++) {
+    if (array[j].value <= pivot) {
+      i += 1;
+      temp = array[i];
+      array[i] = array[j];
+      array[j] = temp;
+    }
+  }
+
+  U64 final_pivot_pos = i + 1;
+  temp = array[final_pivot_pos];
+  array[final_pivot_pos] = array[high];
+  array[high] = temp;
+
+  return final_pivot_pos;
+}
+
+function void qsort_F64(SortPair_F64 *array, U64 size, U64 low, U64 high) {
+  if (low < high) {
+    U64 pivot_index = qsort_partition_F64(array, size, low, high);
+    qsort_F64(array, size, low, pivot_index - 1);
+    qsort_F64(array, size, pivot_index + 1, high);
+  }
+}
+
+function void sort_and_get_indices_F64(F64 *array, U64 *indices, U64 size) {
+  SortPair_F64 *pairs = malloc(size * sizeof(SortPair_F64));
+  for (U64 i = 0; i < size; i++) {
+    pairs[i].value = array[i];
+    pairs[i].original_index = i;
+  }
+
+  qsort_F64(pairs, size, 0, size - 1);
+
+  for (U32 i = 0; i < size; i++) {
+    array[i] = pairs[i].value;
+    indices[i] = pairs[i].original_index;
+  }
+
+  free(pairs);
+}
+
+function void sort_by_indices_F64(F64 *array, U64 *indices, U64 size) {
+
+  F64 *temp = malloc(size * sizeof(F64));
+  for (U64 i = 0; i < size; i++) {
+    U64 original_index = indices[i];
+    temp[i] = array[original_index];
+  }
+  for (U64 i = 0; i < size; i++) {
+    array[i] = temp[i];
+  }
+  free(temp);
+}
 
 /* Auxiliary routine: printing a matrix */
-function void print_eigenvalues(S32 l, S32 n, F64 *wr, F64 *wi) {
+function void print_eigenvalues(char *desc, int n, F64 *wr, F64 *wi) {
   ArenaTemp scratch = scratch_get(0, 0);
-  S32 i, j;
+  int i, j;
   String8 newline = str8_lit("\n");
-  String8 header = str8_pushf(scratch.arena, "\n ------- \n"
+  String8 header = str8_pushf(scratch.arena, "\n %s\n", desc);
-                              "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) Hartree, %4.5f eV \n", 
+        str8_pushf(scratch.arena, " (%4.5f, %4.5f)\n", wr[j], wi[j]);
-                   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 );
   }
@@ -62,13 +134,20 @@ function void print_eigenvalues(S32 l, S32 n, F64 *wr, F64 *wi) {
   scratch_release(scratch);
 }
 
-function void 
+function void set_up_first_matrices(Mat_F64 *H, Mat_F64 *H_l, Mat_F64 *B_inv) {
-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;
@@ -78,7 +157,7 @@ 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
@@ -148,7 +227,7 @@ 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);
   }
@@ -188,9 +267,6 @@ function void compute_wf_norm_F64(F64 *coeffs, U64 coeff_size, U64 n, U64 l) {
   scratch_release(scratch);
 }
 
-/////////////////
-//~
-// Main entry point
 function void EntryPoint(void) {
 
   OS_InitReceipt os_receipt = OS_init();
@@ -203,9 +279,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), 
+  write_array_binary_F64(str8_lit(grid_file_path_bin), g_grid.points,
-                         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");
@@ -223,14 +299,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(g_base_arena, mat_size1, mat_size2);
+  Mat_F64 H_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_base = mat_F64(mat_size1, mat_size2);
-  Mat_F64 H_l = mat_F64(g_base_arena, mat_size1, mat_size2);
+  Mat_F64 H_l = mat_F64(mat_size1, mat_size2);
-  Mat_F64 H = mat_F64(g_base_arena, mat_size1, mat_size2);
+  Mat_F64 H = mat_F64(mat_size1, mat_size2);
   // 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(mat_size1, mat_size2);
   // 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(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,
@@ -238,12 +314,10 @@ 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 = g_angular_momenta[ang_mom_idx];
+    F64 l = 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;
@@ -259,10 +333,12 @@ 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
     {
@@ -293,12 +369,13 @@ function void EntryPoint(void) {
         LOG("Failed to compute eigenvalues in dgeev\n");
         exit(1);
       }
-      
+
-      F64 *right_eigenvectors = vr;
+      eigensolution.right_eigenvectors = mat_F64_from_data(ldvr, size1, 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((U32)l, size1, wr, wi );
+      // print_eigenvalues( "Eigenvalues sorted: ", size1, wr, wi );
 
       U32 i = 0;
       F64 energy = -1000.0;
@@ -314,11 +391,9 @@ 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),
-            eigvecs, size1,
+            eigensolution.right_eigenvectors.matrix[energy_index], size1,
             "%13.6e\n");
 
         i += 1;
@@ -330,7 +405,5 @@ 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(;;) {
+        for(;;) {
-        U64 unread = (U64)(one_past_last - ptr);
+            U64 unread = (U64)(one_past_last - ptr);
-        DWORD to_read = (DWORD)(ClampTop(unread, U32Max));
+            DWORD to_read = (DWORD)(ClampTop(unread, U32Max));
-        DWORD did_read = 0;
+            DWORD did_read = 0;
-        // TODO(anton): Understand WINAPI
+            // TODO(anton): Understand WINAPI
-        if(!ReadFile(handle, ptr, to_read, &did_read, 0)) {
+            if(!ReadFile(handle, ptr, to_read, &did_read, 0)) {
-            break;
+                break;
+            }
+            ptr += did_read;
+            result.size += did_read;
+            if(ptr >= one_past_last) {
+                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))
         {
-            goto fail_out;
+            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;
+            }
         }
-        ptr += did_write;
-        if(ptr >= opl)
-        {
-            break;
-        }
-      }
     }
     fail_out:;
 }