.clangd

@@ -0,0 +1,24 @@
+CompileFlags:
+  Add:
+    - -xc
+    - -std=c11
+    - -Wno-unused-includes
+    - -IE:/dev/hf_again/src
+    - -IE:/lib/intel_mkl/mkl/2025.3/include
+    - -include
+    - E:/dev/hf_again/src/base/base_inc.h
+    - -include
+    - E:/dev/hf_again/src/os/os_inc.h
+    - -include
+    - E:/dev/hf_again/src/hf/bsplines_and_grid.h
+    - -include
+    - E:/dev/hf_again/src/hf/file_io.h
+    - -include
+    - E:/dev/hf_again/src/hf/hf_base.h
+    - -include
+    - E:/dev/hf_again/src/hf/hf_tests.c
+  Remove: [-xobjective-c++-header]
+  Compiler: clang
+
+Index:
+  Background: Build

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

build.bat

@@ -32,6 +32,16 @@ cl %CommonCompilerFlags% %Sources% /I"%mkl_root%\include" ^
 set LastError=%ERRORLEVEL%
 popd
 
+(
+  echo [
+  echo   {
+  echo     "directory": "E:/dev/hf_again/src",
+  echo     "command": "clang -xc -std=c11 -g -O0 -Wno-unused-includes -IE:/dev/hf_again/src -I%mkl_root%/include main.c",
+  echo     "file": "E:/dev/hf_again/src/main.c"
+  echo   }
+  echo ]
+) > compile_commands.json
+
 
 echo Build complete
 

compile_commands.json

@@ -0,0 +1,7 @@
+[
+  {
+    "directory": "E:/dev/hf_again/src",
+    "command": "clang -xc -std=c11 -g -O0 -Wno-unused-includes -IE:/dev/hf_again/src -IE:/lib/intel_mkl/mkl/2025.3/include main.c",
+    "file": "E:/dev/hf_again/src/main.c"
+  }
+]

src/base/base_inc.c

@@ -2,4 +2,3 @@
 #include "base_memory.c"
 #include "base_strings.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 *pairs, U64 size, U64 low, U64 high) {
-  F64 pivot = pairs[high].value;
-  U64 i = low - 1;
-
-  SortPair_F64 temp = {0};
-  for (U64 j = low; j < high; j++) {
-    if (pairs[j].value <= pivot) {
-      i += 1;
-      temp = pairs[i];
-      pairs[i] = pairs[j];
-      pairs[j] = temp;
-    }
-  }
-
-  U64 final_pivot_pos = i + 1;
-  temp = pairs[final_pivot_pos];
-  pairs[final_pivot_pos] = pairs[high];
-  pairs[high] = temp;
-
-  return final_pivot_pos;
-}
-
-function void qsort_F64(SortPair_F64 *pairs, U64 size, U64 low, U64 high) {
-  if (low < high) {
-    U64 pivot_index = qsort_partition_F64(pairs, size, low, high);
-    qsort_F64(pairs, size, low, pivot_index - 1);
-    qsort_F64(pairs, 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/bsplines_and_grid.c

@@ -1,27 +1,16 @@
-
-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;
 
-function void
-bspline_ctx_assign(BSplineCtx *ctx) {
-  g_bspline_ctx = ctx;
-}
-
-function void 
-grid_assign(Grid *grid) {
-  g_grid = grid;
-}
 
 function F64 
 bspline_recursion(F64 x, U32 k, U32 i)
 {
-  F64 *t = g_bspline_ctx->knotpoints;
+  F64 *t = g_bspline_ctx.knotpoints;
 
-  F64 tolerance = 1e-14;
   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 && x == g_grid.end) 
     {
       // TODO(anton):
       // This is like a hack to get the last bspline to be 1 at the last point.
@@ -29,7 +18,7 @@ bspline_recursion(F64 x, U32 k, U32 i)
       // to unity at the last point, actually. I have to check this.
       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;
     } else {
@@ -50,13 +39,15 @@ bspline_recursion(F64 x, U32 k, U32 i)
 
     return term1 + term2;
   }
+
+
 }
 
 
 function F64 
 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);
   return out;
 }
@@ -65,9 +56,9 @@ compute_bspline_F64(F64 x_coord, U32 index)
 function F64
 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 *t = g_bspline_ctx->knotpoints;
+  F64 *t = g_bspline_ctx.knotpoints;
   F64 term1_enum = bspline_recursion(x_coord, k-1, index);
   F64 term2_enum = bspline_recursion(x_coord, k-1, index+1);
   F64 term1_denom = t[index+k-1] - t[index];
@@ -81,17 +72,17 @@ compute_dBspline_F64(F64 x_coord, U32 index)
 function void
 set_up_grid(Arena *arena)
 {
-  g_grid->start = GRID_START_POINT;
-  g_grid->end = GRID_END_POINT;
-  g_grid->num_steps = GRID_NUM_STEPS;
+  g_grid.start = GRID_START_POINT;
+  g_grid.end = GRID_END_POINT;
+  g_grid.num_steps = 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;
-  g_grid->points[0] = g_grid->start;
-  g_grid->points[g_grid->num_steps-1] = g_grid->end;
-  for(U32 i = 1; i < g_grid->num_steps-1; i++)
+  g_grid.points = PushArray(arena, 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[g_grid.num_steps-1] = g_grid.end;
+  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;
   }
 
 }
@@ -106,13 +97,13 @@ get_bspline_index_size(U32 size1, U32 i, U32 j)
 function inline U32
 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
 get_bspline_grid_index(U32 i, U32 j)
 {
-  return g_grid->num_steps * j + i;
+  return g_grid.num_steps * j + i;
 }
 
 
@@ -121,36 +112,31 @@ set_up_bspline_context(Arena* arena)
 {
   // Create knotpoint sequence.
   U32 k = 4;
-  U32 bspl_N = BSPLINES_NUM;
-  g_bspline_ctx->order = k; 
-  g_bspline_ctx->num_knotpoints = bspl_N; 
-  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->arena = arena;
-  g_bspline_ctx->knotpoints = PushArray(arena, F64, g_bspline_ctx->num_knotpoints);
+  U32 bspl_N = 40;
+  g_bspline_ctx.order = k; 
+  g_bspline_ctx.num_knotpoints = bspl_N; 
+  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.arena = arena;
+  g_bspline_ctx.knotpoints = PushArray(arena, F64, g_bspline_ctx.num_knotpoints);
   // 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
-  U32 first_phys_index = k-1;
-  U32 phys_point_last_index = g_bspline_ctx->num_phys_points + k-1;
-  for(U32 i = 0; i <= first_phys_index; i++) 
+  U32 phys_point_last_index = g_bspline_ctx.num_phys_points + k-1;
+  for(U32 i = 0; i < k; i++) 
   {
-    g_bspline_ctx->knotpoints[i] = g_grid->start;
+    g_bspline_ctx.knotpoints[i] = g_grid.start;
   }
-  g_bspline_ctx->first_phys_index = first_phys_index;
   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. including last physical point
-  F64 last_physical = g_grid->end;
-  for(U32 i = phys_point_last_index; i < g_bspline_ctx->num_knotpoints; i++)
+  // Set the last points
+  F64 last_physical = g_grid.end;
+  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;
   }
-  g_bspline_ctx->last_phys_index = phys_point_last_index;
 }
 
 function void 
@@ -201,11 +187,11 @@ function void
 set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena)
 {
 
-  U64 num_bsplines = g_bspline_ctx->num_bsplines;
-  U64 k = g_bspline_ctx->order;
-  F64 *t = g_bspline_ctx->knotpoints;
-  U32 num_grid_points = g_grid->num_steps;
-  U32 num_knotpoints = g_bspline_ctx->num_knotpoints;
+  U64 num_bsplines = g_bspline_ctx.num_bsplines;
+  U64 k = g_bspline_ctx.order;
+  F64 *t = g_bspline_ctx.knotpoints;
+  U32 num_grid_points = g_grid.num_steps;
+  U32 num_knotpoints = g_bspline_ctx.num_knotpoints;
 
   // For sanity check we make the first 4 bsplines by hand.
   {
@@ -221,7 +207,7 @@ set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena)
     F64 *dBspl9 = PushArray(arena, F64, num_grid_points);
     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);
       bspl1[i] = compute_bspline_F64(x, 1);
       bspl2[i] = compute_bspline_F64(x, 2);
@@ -234,35 +220,35 @@ set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena)
       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\\bspline1.dat"), bspl1, num_grid_points, "%13.6e\n");
-    write_array_F64(str8_lit("E:\\dev\\hf_again\\out\\bspline2.dat"), bspl2, num_grid_points, "%13.6e\n");
-    write_array_F64(str8_lit("E:\\dev\\hf_again\\out\\bspline3.dat"), bspl3, num_grid_points, "%13.6e\n");
-    write_array_F64(str8_lit("E:\\dev\\hf_again\\out\\bspline9.dat"), bspl9, num_grid_points, "%13.6e\n");
-    write_array_F64(str8_lit("E:\\dev\\hf_again\\out\\dBspline0.dat"), dBspl0, num_grid_points, "%13.6e\n");
-    write_array_F64(str8_lit("E:\\dev\\hf_again\\out\\dBspline1.dat"), dBspl1, num_grid_points, "%13.6e\n");
-    write_array_F64(str8_lit("E:\\dev\\hf_again\\out\\dBspline2.dat"), dBspl2, num_grid_points, "%13.6e\n");
-    write_array_F64(str8_lit("E:\\dev\\hf_again\\out\\dBspline3.dat"), dBspl3, num_grid_points, "%13.6e\n");
-    write_array_F64(str8_lit("E:\\dev\\hf_again\\out\\dBspline9.dat"), dBspl9, num_grid_points, "%13.6e\n");
+    write_array_F64(str8_lit("D:\\dev\\hf_again\\out\\bspline0.dat"), bspl0, num_grid_points, "%13.6e\n");
+    write_array_F64(str8_lit("D:\\dev\\hf_again\\out\\bspline1.dat"), bspl1, num_grid_points, "%13.6e\n");
+    write_array_F64(str8_lit("D:\\dev\\hf_again\\out\\bspline2.dat"), bspl2, num_grid_points, "%13.6e\n");
+    write_array_F64(str8_lit("D:\\dev\\hf_again\\out\\bspline3.dat"), bspl3, num_grid_points, "%13.6e\n");
+    write_array_F64(str8_lit("D:\\dev\\hf_again\\out\\bspline9.dat"), bspl9, num_grid_points, "%13.6e\n");
+    write_array_F64(str8_lit("D:\\dev\\hf_again\\out\\dBspline0.dat"), dBspl0, num_grid_points, "%13.6e\n");
+    write_array_F64(str8_lit("D:\\dev\\hf_again\\out\\dBspline1.dat"), dBspl1, num_grid_points, "%13.6e\n");
+    write_array_F64(str8_lit("D:\\dev\\hf_again\\out\\dBspline2.dat"), dBspl2, num_grid_points, "%13.6e\n");
+    write_array_F64(str8_lit("D:\\dev\\hf_again\\out\\dBspline3.dat"), dBspl3, num_grid_points, "%13.6e\n");
+    write_array_F64(str8_lit("D:\\dev\\hf_again\\out\\dBspline9.dat"), dBspl9, num_grid_points, "%13.6e\n");
   }
 
   {
     ArenaTemp scratch = scratch_get(0,0);
 
-  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.bsplines_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 i = 0; i < num_grid_points; i++)
     {
       U32 index = get_bspline_grid_index(i, j);
-      F64 bspline_value = compute_bspline_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->dBsplines_grid[index] = dBspline_value;
+      F64 bspline_value = compute_bspline_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.dBsplines_grid[index] = dBspline_value;
       if(g_debug_bspline_matrix && j == 0) 
       {
         String8 out = str8_pushf(scratch.arena, "%i %i \t %13.6e\n", i, index, bspline_value); 
@@ -274,35 +260,35 @@ set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena)
     scratch_release(scratch);
   }
 
-  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.bsplines = 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 j = 0; j < num_bsplines; j++)
     {
       U32 index = get_bspline_index(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.bsplines[index] = compute_bspline_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_bsplines, 
       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_bsplines, 
       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_bsplines, 
       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_bsplines, 
       str8_lit(dBspline_knots_array_file_path));

src/hf/bsplines_and_grid.h

@@ -3,12 +3,8 @@
 
 // coordinates in Bohr radii
 #define GRID_START_POINT 0.0 
-#define GRID_END_POINT 300.0
-
-#define GRID_NUM_STEPS 100 // Just for plotting/writing, 
-                           // else we stick to knotpts and quadrature
-
-#define BSPLINES_NUM 2000
+#define GRID_END_POINT 40.0
+#define GRID_NUM_STEPS 300
 
 typedef struct Grid Grid;
 struct Grid 
@@ -28,8 +24,6 @@ struct BSplineCtx
   U32 num_knotpoints;
   U32 num_bsplines;
   U32 num_phys_points;
-  U32 first_phys_index;
-  U32 last_phys_index;
   F64 *bsplines_grid; // In grid points
   F64 *dBsplines_grid; // First deriv in grid points
   F64 *bsplines; // In knotpoints
@@ -47,6 +41,5 @@ function inline U32 get_bspline_grid_index(U32 i, U32 j);
 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 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 */

src/hf/file_io.c

@@ -1,5 +1,4 @@
-#define DEBUG_LOG_WRITE_STRING_LIST_TO_FILE 1
-#define DEBUG_LOG_WRITE_ARRAY_BINARY 1
+
 
 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
-  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
-
+  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);
+  }
   OS_file_close(file_handle);
 }
 

src/hf/file_io.h

@@ -1,13 +1,13 @@
 #ifndef FILE_IO_H
 #define FILE_IO_H
 
-#define grid_file_path_bin "E:\\dev\\hf_again\\out\\grid.bin"
-#define grid_file_path "E:\\dev\\hf_again\\out\\grid.dat"
-#define knotpoints_file_path "E:\\dev\\hf_again\\out\\knotpoints.dat"
-#define bspline_grid_array_file_path "E:\\dev\\hf_again\\out\\bsplines_grid.dat"
-#define dBspline_grid_array_file_path "E:\\dev\\hf_again\\out\\dBsplines_grid.dat"
-#define bspline_knots_array_file_path "E:\\dev\\hf_again\\out\\bsplines_knots.dat"
-#define dBspline_knots_array_file_path "E:\\dev\\hf_again\\out\\dBsplines_knots.dat"
+#define grid_file_path_bin "D:\\dev\\hf_again\\out\\grid.bin"
+#define grid_file_path "D:\\dev\\hf_again\\out\\grid.dat"
+#define knotpoints_file_path "D:\\dev\\hf_again\\out\\knotpoints.dat"
+#define bspline_grid_array_file_path "D:\\dev\\hf_again\\out\\bsplines_grid.dat"
+#define dBspline_grid_array_file_path "D:\\dev\\hf_again\\out\\dBsplines_grid.dat"
+#define bspline_knots_array_file_path "D:\\dev\\hf_again\\out\\bsplines_knots.dat"
+#define dBspline_knots_array_file_path "D:\\dev\\hf_again\\out\\dBsplines_knots.dat"
 
 
 function void write_string_list_to_file(Arena *arena, String8 path, String8List *list);

src/hf/hf_base.c

@@ -1,43 +1,79 @@
-//~ Globals
-GaussLegendre g_gauss_legendre = {0};
 
-function inline U32
-mat_ij_to_index(U32 i, U32 j, U32 size1) {
-  return j * size1 + i;
+//~ Globals
+global GaussLegendre g_gauss_legendre = {0};
+
+
+function Mat_F64 
+mat_F64(U32 size1, U32 size2)
+{
+  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(Arena *arena, U32 size1, U32 size2)
+mat_F64_from_data(U32 size1, U32 size2, F64 *data)
 {
   Mat_F64 out = {0};
-  out.arena = arena;
   out.size1 = size1;
   out.size2 = size2;
-  out.data = PushArray(arena, F64, size1*size2);
-  MemoryZero(out.data, size1*size2*sizeof(F64));
+  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_ij_to_index(i, j, mat->size1);
-  mat->data[index] = val;
+  U32 row = i < mat->size1 ? i : mat->size1-1;
+  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_ij_to_index(i, j, mat->size1);
-  return mat->data[index];
+  U32 row = i < mat->size1 ? i : mat->size1-1;
+  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;
@@ -46,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
-  // 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));
+  U64 data_byte_size = src->size1 * src->size2 * sizeof(F64);
+  MemoryCopy(dst->data, src->data, data_byte_size);
 }
 
 function void
@@ -192,473 +225,10 @@ mat_invert_F64(Mat_F64 *mat)
 }
 
 
-function F64
-hartree2eV(F64 energy_hartree) {
-  return 27.2114*energy_hartree;
-}
-
-function Problem problem_create(void) {
-  Problem out = {0};
-  out.arena = m_make_arena();
-  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 */
-function void print_eigenvalues(S32 l, S32 n, F64 *wr, F64 *wi) {
-  ArenaTemp scratch = scratch_get(0, 0);
-  S32 i, j;
-  String8 newline = str8_lit("\n");
-  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) 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 );
-  }
-  LOG(newline.str);
-  // printf("\n");
-
-  scratch_release(scratch);
-}
-
-function void 
-set_up_base_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
-  // 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);
-    BSplineCtx *bspl_ctx = &problem->bspline_ctx;
-    F64 *t = bspl_ctx->knotpoints;
-    F64 Z = (F64)problem->atom.Z;
-    U32 k = bspl_ctx->order;
-
-    // Skipping first bspline
-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.
-        U32 bspl_index_j = j + 1;
-
-        // This logic assumes 1-indexed bsplines
-        F64 abs_index_diff =
-            fabs((F64)(bspl_index_i + 1) - (F64)(bspl_index_j + 1));
-        if (!(abs_index_diff > ((F64)k - 1.0))) {
-          // We do Gaussian quadrature between each knot point,
-          // so we need to figure out where to start.
-          // We start integration in the first shared knotpoint, which is the
-          // one of the highest index.
-          U32 start_knotpoint_index =
-              bspl_index_i < bspl_index_j ? bspl_index_j : bspl_index_i;
-          // And we integrate over the next k knotpoints.
-          U32 end_knotpoint_index =
-              bspl_index_i < bspl_index_j ? bspl_index_i + k : bspl_index_j + k;
-
-          F64 term1 = 0.0;
-          F64 term2 = 0.0;
-          F64 term3 = 0.0;
-          F64 Bmat_term = 0.0;
-
-          for (U32 knotpoint_idx = start_knotpoint_index;
-               knotpoint_idx < end_knotpoint_index; knotpoint_idx++) {
-            F64 a = t[knotpoint_idx];
-            F64 b = t[knotpoint_idx + 1];
-            F64 prefac = 0.5 * (b - a);
-
-            // Only integrate non-zero intervals
-            if (prefac > 1e-16) {
-              for (U32 gq_i = 0; gq_i < g_gauss_legendre.order; gq_i++) {
-                F64 w = g_gauss_legendre.weights[gq_i];
-                F64 z = g_gauss_legendre.abscissae[gq_i];
-                F64 r = (z * prefac) + ((a + b) * 0.5);
-                F64 term_prefac = (prefac * w);
-                F64 dB_i = compute_dBspline_F64(r, bspl_index_i);
-                F64 dB_j = compute_dBspline_F64(r, bspl_index_j);
-                F64 B_i = compute_bspline_F64(r, bspl_index_i);
-                F64 B_j = compute_bspline_F64(r, bspl_index_j);
-                term1 += term_prefac * dB_i * dB_j;
-                term2 += term_prefac * B_i * B_j / (r * r);
-                term3 += term_prefac * B_i * B_j / r;
-                Bmat_term += term_prefac * B_i * B_j;
-              }
-            }
-          }
-
-          F64 H_term_sum = 0.5 * term1 + (-Z) * term3;
-          F64 H_l_term = 0.5 * term2;
-          /* String8 debug = str8_pushf(scratch.arena,
-           * "(i=%i,j=%i,t_i=%4.4f,t_i=%4.4f,term1=%.4e,term2=%.4e,term3=%.4e,term_sum=%.4e)
-           * \n", */
-          /*     bspl_index_i, bspl_index_j,
-           * t[bspl_index_i+k-1],t[bspl_index_j+k-1],term1,term2,term3,term_sum);
-           */
-          /* LOG(debug.str); */
-          mat_F64_set(H, i, j, H_term_sum);
-          mat_F64_set(H_l, i, j, H_l_term);
-          mat_F64_set(B_inv, i, j, Bmat_term);
-          // mat_F64_set(&H, i, j, abs_index_diff);
-        }
-        // mat_F64_set(&H, i, j, abs_index_diff);
-      }
-      // LOG("\n");
-    }
-
-    LOG(str8_pushf(scratch.arena, "H.size1=N-k-2=%i, last bspline index=%i \n",
-                   H->size1, bspl_ctx->num_bsplines - 1)
-            .str);
-    scratch_release(scratch);
-    //print_mat_F64(H);
-    LOG("\n");
-    // print_mat_F64(B);
-  }
-}
-
-function void 
-compute_wf_norm_F64(F64 *coeffs, U64 coeff_size, U64 n, U64 l) {
-  ArenaTemp scratch = scratch_get(0, 0);
-
-  // Gauss legendre integration
-  //
-  F64 norm = 0.0;
-  for (U64 i = 0; i < g_grid->num_steps - 1; i++) {
-    F64 a = g_grid->points[i];
-    F64 b = g_grid->points[i + 1];
-    F64 prefac = 0.5 * (b - a);
-
-    // Only integrate non-zero intervals
-    if (prefac > 1e-16) {
-      for (U32 gq_i = 0; gq_i < g_gauss_legendre.order; gq_i++) {
-        F64 w = g_gauss_legendre.weights[gq_i];
-        F64 z = g_gauss_legendre.abscissae[gq_i];
-        F64 r = (z * prefac) + ((a + b) * 0.5);
-        F64 term_prefac = (prefac * w);
-        F64 wf_at_r = 0.0;
-
-        for (U64 j = 0; j < coeff_size; j++) {
-          wf_at_r += coeffs[j] * compute_bspline_F64(r, j + 1);
-        }
-
-        norm += term_prefac * wf_at_r * wf_at_r;
-      }
-    }
-  }
-  String8 out =
-      str8_pushf(scratch.arena, "n:%i, l:%i norm: %.2f \n", n, l, norm);
-  LOG(out.str);
-  scratch_release(scratch);
-}
-
-function F64 naive_kernel(F64 r1, F64 r2, F64 k, F64 alpha, F64 beta) {
-  if(r1 < r2) {
-    return pow(r1, 2.0*k+1.0)*exp(-alpha*r1)*exp(-beta*r2);
-  } else {
-    return exp(-alpha*r1)*pow(r2, 2.0*k+1.0)*exp(-beta*r2);
-  }
-}
-
-function F64 factorial(S32 k) {
-  if(k == 0) {
-    return 1.0;
-  } else if(k == 1) {
-    return 1.0;
-  } else {
-    F64 kf = (F64)k;
-    S32 next = k - 1;
-    return kf * factorial(next);
-  }
-}
-
-/////////////////////////////////////////////////////////////////////////////////////////
-//~ Testing HF integrals.
-function void 
-hf_coulomb_integrals_test(void) 
-{
-  
-  Problem problem = problem_create();
-  LOG(str8_pushf(problem.arena, "Created Problem-struct for %s \n", problem.atom.name.str).str);
-
-  set_up_gauss_legendre_points(problem.arena);
-
-  //- Set up grid and write to file.
-  grid_assign(&problem.grid);
-  set_up_grid(problem.arena);
-
-  write_array_binary_F64(str8_lit(grid_file_path_bin), 
-                         problem.grid.points,
-                         problem.grid.num_steps);
-
-  write_array_F64(str8_lit(grid_file_path), 
-                  problem.grid.points,
-                  problem.grid.num_steps,
-                  "%13.6e\n");
-
-  //- The BSpline context is the knotpoints and the BSpline order etc.
-  bspline_ctx_assign(&problem.bspline_ctx);
-  set_up_bspline_context(problem.arena);
-  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
-  // debugging.
-  set_up_bsplines_at_points_and_write_matrix_F64(problem.arena);
-
-  U32 N = problem.bspline_ctx.num_knotpoints;
-  U32 k = problem.bspline_ctx.order;
-  U32 mat_size1 = N - k - 2;
-  U32 mat_size2 = mat_size1;
-  
-  LOG(" *** Running tests for coulomb operator integrals! *** \n ");
-  Mat_F64 test_mat = mat_F64(problem.arena, mat_size1, mat_size2);
-  Mat_F64 double_int_mat = mat_F64(problem.arena, mat_size1, mat_size2);
-
-  // We want to compute integrals over two r-coordinates r1 and r2, of the form
-  //   int1(int2( F(r1) (r^k_smaller / r^(k+1)_greater) G(r2) )
-  // We will test with analytic functions
-  //   F(r) = r^(k+1) exp(-alpha*r)
-  //   G(r) = r^(k+1) exp(-beta*r)
-  // Because then whenever r1 < r2 the kernel will be
-  //   r1^(2k+1)exp(-alpha*r1)*exp(-beta*r2)
-  // and whenever r2 > r1 we have
-  //   exp(-alpha*r1)*r2^(2k+1)exp(-beta*r2)
-  // That is we cancel the denominator.
-  // In particular we can get an analytic result as a function of k, alpha and beta:
-  //   (2k+1)! / ((alpha*beta)(alpha+beta)^(2k+1))
-  // The reason is so that we can
-  //  1) Make sure that we are integrating over two r-coordinates properly.
-  //  2) Compare a "fast" method to just brute force double integral loop to make sur
-  //     it is correct.
-  //  3) Expand the method to the general spline case to compute the Coulomb matrix 
-  //     element with confidence that it is correct.
-  {
-    BSplineCtx *bspl_ctx = &problem.bspline_ctx;
-    U32 bspl_order = bspl_ctx->order;
-    F64 *t1 = bspl_ctx->knotpoints;
-    F64 k_mpole = 1.0;
-
-    F64 alpha = 2.0;
-    F64 beta = 3.0;
-    // The set of r coordinates will be all knotpoints and the quadrature points in
-    // between them
-    F64 result = 0.0;
-    U32 phys_point_first_index = bspl_ctx->first_phys_index;;
-    U32 phys_point_last_index = bspl_ctx->num_phys_points + k-1;
-    for(S32 r1_i = phys_point_first_index; r1_i < phys_point_last_index-1; r1_i += 1) 
-    {
-      F64 a1 = t1[r1_i];
-      F64 b1 = t1[r1_i + 1];
-      F64 c1 = (a1+b1);
-      F64 h1 = (b1-a1);
-      for(S32 r2_i = phys_point_first_index; r2_i < phys_point_last_index-1; r2_i += 1) 
-      {
-        F64 a2 = t1[r2_i];
-        F64 b2 = t1[r2_i + 1];
-        F64 c2 = (a2+b2);
-        F64 h2 = (b2-a2);
-        for(S32 q1 = 0; q1 < g_gauss_legendre.order; q1 += 1) {
-          F64 w1 = 0.5 * h1 * g_gauss_legendre.weights[q1];
-          F64 z1 = g_gauss_legendre.abscissae[q1];
-          F64 r1 = (z1 * (0.5 * h1)) + (c1 * 0.5);
-          for(S32 q2 = 0; q2 < g_gauss_legendre.order; q2 += 1) {
-            F64 w2 = 0.5 * h2 * g_gauss_legendre.weights[q2];
-            F64 z2 = g_gauss_legendre.abscissae[q2];
-            F64 r2 = (z2 * (0.5 * h2)) + (c2 * 0.5);
-            F64 kernel = naive_kernel(r1, r2, k_mpole, alpha, beta);
-            result += w1 * w2 * kernel;
-          }
-        }
-      }
-    }
-
-    ArenaTemp scratch = scratch_get(0,0);
-    // Analytic 
-    S32 k_mpole_int = 1;
-    F64 fact = factorial(2*k_mpole_int + 1);
-    F64 reference = fact / ((alpha*beta)*pow((alpha+beta),2.0*k_mpole+1.0));
-    LOG(str8_pushf(scratch.arena, "Reference: %.12f \n", reference).str);
-    LOG(str8_pushf(scratch.arena, "Computed: %.12f \n", result).str);
-    
-    scratch_release(scratch);
-  }
-
-}
-
-/////////////////////////////////////////////////////////////////////////////////////////
-//~ 
-// Main problem function, called from entry point.
-function void
-hf_eigenvalues_main(void) 
-{
-  
-  Problem problem = problem_create();
-  LOG(str8_pushf(problem.arena, "Created Problem-struct for %s \n", problem.atom.name.str).str);
-
-  set_up_gauss_legendre_points(problem.arena);
-
-  //- Set up grid and write to file.
-  grid_assign(&problem.grid);
-  set_up_grid(problem.arena);
-
-  write_array_binary_F64(str8_lit(grid_file_path_bin), 
-                         problem.grid.points,
-                         problem.grid.num_steps);
-
-  write_array_F64(str8_lit(grid_file_path), 
-                  problem.grid.points,
-                  problem.grid.num_steps,
-                  "%13.6e\n");
-
-  //- The BSpline context is the knotpoints and the BSpline order etc.
-  bspline_ctx_assign(&problem.bspline_ctx);
-  set_up_bspline_context(problem.arena);
-  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
-  // debugging.
-  set_up_bsplines_at_points_and_write_matrix_F64(problem.arena);
-
-  U32 N = problem.bspline_ctx.num_knotpoints;
-  U32 k = problem.bspline_ctx.order;
-  U32 mat_size1 = N - k - 2;
-  U32 mat_size2 = mat_size1;
-  problem.H_base = mat_F64(problem.arena, mat_size1, mat_size2);
-  problem.H_l_base = mat_F64(problem.arena, mat_size1, mat_size2);
-  problem.H_l = mat_F64(problem.arena, mat_size1, mat_size2);
-  problem.H = mat_F64(problem.arena, mat_size1, mat_size2);
-  // This will be the inverse of B, but to start with we construct B.
-  problem.B_inv = mat_F64(problem.arena, mat_size1, mat_size2);
-  // A is the actual matrix for each eigenvalue problem.
-  problem.A = mat_F64(problem.arena, problem.H.size1, problem.H.size2);
-  set_up_base_matrices(&problem, 
-                       &problem.H_base,
-                       &problem.H_l_base,
-                       &problem.B_inv);
-
-  // 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
-  // zgeev
-  mat_invert_F64(&problem.B_inv);
-
-  // This arena is used to push results from f. ex eigenvalue computations.
-  // For each angular momentum
-  for (U32 ang_mom_idx = 0; ang_mom_idx < MAX_NUM_ANGULAR_MOMENTA; ang_mom_idx++) {
-    ArenaTemp scratch = scratch_get(0, 0);
-    mat_F64_copy_to_dst(&problem.H, &problem.H_base);
-    F64 l = problem.angular_momentum_l[ang_mom_idx];
-    Eigensolution_F64 *eigsol = &problem.eigsols[ang_mom_idx];
-    eigsol->l = (U32)l;
-
-    if (l > 1e-16) {
-      F64 l_factor = l * (l + 1.0);
-      U64 mat_size = problem.H_l.size1 * problem.H_l.size2;
-      mat_F64_copy_to_dst(&problem.H_l, &problem.H_l_base);
-      // Multiply l(l+1)
-      cblas_dscal(mat_size, l_factor, problem.H_l.data, 1);
-      // Add H = H_base + H_l
-      cblas_daxpy(mat_size, 1.0, problem.H_l.data, 1, problem.H.data, 1);
-    }
-    //LOG("Matrix H: \n");
-    //print_mat_F64(&problem.H);
-
-    // Multiply to get A = B^-1 H
-    {
-      S32 n = problem.A.size1;
-      cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n, n, n, 1.0,
-                  problem.B_inv.data, n, problem.H.data, n, 0.0, problem.A.data, n);
-      //LOG("Matrix A: \n");
-      //print_mat_F64(&problem.A);
-    }
-
-    // Solve generalised eigenvalue problem
-    {
-      S32 size1 = problem.A.size1;
-      S32 lda = size1;
-      S32 ldvl = size1;
-      S32 ldvr = size1;
-      S32 info;
-      S32 lwork;
-
-      F64 wkopt;
-      F64 *work;
-
-      eigsol->eigvals_re= PushArray(problem.arena, F64, size1);
-      F64 *wr = eigsol->eigvals_re;
-      eigsol->eigvals_im = PushArray(problem.arena, F64, size1);
-      F64 *wi = eigsol->eigvals_im;
-      eigsol->left_eigvecs = mat_F64(problem.arena, ldvl, size1);
-      F64 *vl = eigsol->left_eigvecs.data;
-      eigsol->right_eigvecs = mat_F64(problem.arena, size1, ldvr);
-      F64 *vr = eigsol->right_eigvecs.data;
-
-      lwork = -1;
-      F64 *a = problem.A.data;
-      dgeev("Vectors", "Vectors", &size1, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr,
-            &wkopt, &lwork, &info);
-      lwork = (S32)wkopt;
-      //work = (F64 *)malloc(lwork * sizeof(F64));
-      work = PushArray(scratch.arena, F64, lwork);
-      dgeev("Vectors", "Vectors", &size1, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr,
-            work, &lwork, &info);
-      if (info > 0) {
-        LOG("Failed to compute eigenvalues in dgeev\n");
-        exit(1);
-      }
-      
-      // Sort real and imaginary eigenvalues by real part
-      U64 *sorted_indices = PushArray(scratch.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 );
-
-      U32 i = 0;
-      F64 energy = -1.0*(F64)(U32Max); // Just have an unreasonable energy bound
-      U32 counter = 0;
-      while (energy < 0.0) {
-        energy = wr[i];
-        U64 energy_index = sorted_indices[i];
-        U64 n = 1 + i;
-        if (ang_mom_idx > 0) {
-          n = 2 + i;
-        }
-
-        // compute_wf_norm_F64(eigensolution.right_eigenvectors.matrix[energy_index],
-        // size1, n, ang_mom_idx);
-
-        U64 eigvec_idx = mat_ij_to_index(0, energy_index, size1);
-        F64 *eigvecs = &eigsol->right_eigvecs.data[eigvec_idx];
-        write_array_F64(
-            get_eigenvector_filename(scratch.arena, n, ang_mom_idx),
-            eigvecs, size1,
-            "%13.6e\n");
-
-        i += 1;
-        counter += 1;
-        if (counter > 10) {
-          break;
-        }
-      }
-
-      scratch_release(scratch);
-    }
-
-  }
-
-}
+
+
+
+
 
 
 

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;
 };
 
@@ -26,85 +26,17 @@ struct GaussLegendre {
   F64 *abscissae;
 };
 
-// We have one set of eigenvalue problem solutions for each angular momentum 
-// quantum number
-typedef struct Eigensolution_F64 Eigensolution_F64;
-struct Eigensolution_F64 {
-  U32 l;
-  F64 *eigvals_re;
-  F64 *eigvals_im;
-  Mat_F64 right_eigvecs;
-  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
-};
-
-typedef struct Atom Atom;
-struct Atom {
-  String8 name;
-  U32 Z;
-  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;
-  // Problem matrices. Since we can perform the "mass" integral once we call that
-  // "base", and similarly we can compute the integral for the l-dependent part once,
-  // and call that H_l_base. Then we can add whatever numerical factors we want on the 
-  // base matrices, and finally combine them.
-  Mat_F64 H_base;
-  Mat_F64 H;
-  Mat_F64 H_l_base;
-  Mat_F64 H_l;
-  // 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
-  // zgeev for Ac = Ec
-  Mat_F64 B_inv;
-  Mat_F64 A;
-};
-
 //~ Base math and utility functions
 
 // Mat_F64 functions
-function inline U32 mat_ij_to_index(U32 i, U32 j, U32 size1);
-function Mat_F64 mat_F64(Arena* arena, U32 size1, U32 size2);
+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 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);
 
 // Gauss-Legendre
@@ -113,19 +45,5 @@ function void set_up_gauss_legendre_points(Arena *arena);
 // Random utility
 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);
-function void print_eigenvalues(S32 l, S32 n, F64 *wr, F64 *wi);
-function void compute_wf_norm_F64(F64 *coeffs, U64 coeff_size, U64 n, U64 l);
-
-// Problem
-function Problem problem_create();
-function void set_up_base_matrices(Problem *problem,
-                                   Mat_F64 *H,
-                                   Mat_F64 *H_l,
-                                   Mat_F64 *B_inv);
-function void hf_main();
-
-
-
 
 #endif /* HF_BASE_H */

src/main.c

@@ -1,44 +1,409 @@
-//~
-//----
+// ---
 // Header includes
 #include "base/base_inc.h"
 #include "os/os_inc.h"
 
-#include "hf/hf_bsplines_and_grid.h"
-#include "hf/hf_file_io.h"
+#include "hf/bsplines_and_grid.h"
+#include "hf/file_io.h"
 #include "hf/hf_base.h"
 
-//----
+// ---
 // .C includes
 #include "base/base_inc.c"
 #include "os/os_entry_point.c"
 #include "os/os_inc.c"
 
-#include "hf/hf_bsplines_and_grid.c"
-#include "hf/hf_file_io.c"
+#include "hf/bsplines_and_grid.c"
+#include "hf/file_io.c"
 #include "hf/hf_base.c"
 
 // 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;
+  F64 *eigenvalues_im;
+  Mat_F64 right_eigenvectors;
+  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 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(char *desc, int n, F64 *wr, F64 *wi) {
+  ArenaTemp scratch = scratch_get(0, 0);
+  int i, j;
+  String8 newline = str8_lit("\n");
+  String8 header = str8_pushf(scratch.arena, "\n %s\n", desc);
+  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]);
+    LOG(outstr.str);
+    // printf(" (%6.2f,%6.2f)", a[i+j*lda].real, a[i+j*lda].imag );
+  }
+  LOG(newline.str);
+  // printf("\n");
+
+  scratch_release(scratch);
+}
+
+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;
+    U32 k = g_bspline_ctx.order;
+
+    // Skipping first bspline
+    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.
+        U32 bspl_index_j = j + 1;
+
+        // This logic assumes 1-indexed bsplines
+        F64 abs_index_diff =
+            fabs((F64)(bspl_index_i + 1) - (F64)(bspl_index_j + 1));
+        if (!(abs_index_diff > ((F64)k - 1.0))) {
+          // We do Gaussian quadrature between each knot point,
+          // so we need to figure out where to start.
+          // We start integration in the first shared knotpoint, which is the
+          // one of the highest index.
+          U32 start_knotpoint_index =
+              bspl_index_i < bspl_index_j ? bspl_index_j : bspl_index_i;
+          // And we integrate over the next k knotpoints.
+          U32 end_knotpoint_index =
+              bspl_index_i < bspl_index_j ? bspl_index_i + k : bspl_index_j + k;
+
+          F64 term1 = 0.0;
+          F64 term2 = 0.0;
+          F64 term3 = 0.0;
+          F64 Bmat_term = 0.0;
+
+          for (U32 knotpoint_idx = start_knotpoint_index;
+               knotpoint_idx < end_knotpoint_index; knotpoint_idx++) {
+            F64 a = t[knotpoint_idx];
+            F64 b = t[knotpoint_idx + 1];
+            F64 prefac = 0.5 * (b - a);
+
+            // Only integrate non-zero intervals
+            if (prefac > 1e-16) {
+              for (U32 gq_i = 0; gq_i < g_gauss_legendre.order; gq_i++) {
+                F64 w = g_gauss_legendre.weights[gq_i];
+                F64 z = g_gauss_legendre.abscissae[gq_i];
+                F64 r = (z * prefac) + ((a + b) * 0.5);
+                F64 term_prefac = (prefac * w);
+                F64 dB_i = compute_dBspline_F64(r, bspl_index_i);
+                F64 dB_j = compute_dBspline_F64(r, bspl_index_j);
+                F64 B_i = compute_bspline_F64(r, bspl_index_i);
+                F64 B_j = compute_bspline_F64(r, bspl_index_j);
+                term1 += term_prefac * dB_i * dB_j;
+                term2 += term_prefac * B_i * B_j / (r * r);
+                term3 += term_prefac * B_i * B_j / r;
+                Bmat_term += term_prefac * B_i * B_j;
+              }
+            }
+          }
+
+          F64 H_term_sum = 0.5 * term1 + (-Z) * term3;
+          F64 H_l_term = 0.5 * term2;
+          /* String8 debug = str8_pushf(scratch.arena,
+           * "(i=%i,j=%i,t_i=%4.4f,t_i=%4.4f,term1=%.4e,term2=%.4e,term3=%.4e,term_sum=%.4e)
+           * \n", */
+          /*     bspl_index_i, bspl_index_j,
+           * t[bspl_index_i+k-1],t[bspl_index_j+k-1],term1,term2,term3,term_sum);
+           */
+          /* LOG(debug.str); */
+          mat_F64_set(H, i, j, H_term_sum);
+          mat_F64_set(H_l, i, j, H_l_term);
+          mat_F64_set(B_inv, i, j, Bmat_term);
+          // mat_F64_set(&H, i, j, abs_index_diff);
+        }
+        // mat_F64_set(&H, i, j, abs_index_diff);
+      }
+      // LOG("\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)
+            .str);
+    scratch_release(scratch);
+    print_mat_F64(H);
+    LOG("\n");
+    // print_mat_F64(B);
+  }
+}
+
+function void compute_wf_norm_F64(F64 *coeffs, U64 coeff_size, U64 n, U64 l) {
+  ArenaTemp scratch = scratch_get(0, 0);
+
+  // Gauss legendre integration
+  //
+  F64 norm = 0.0;
+  for (U64 i = 0; i < g_grid.num_steps - 1; i++) {
+    F64 a = g_grid.points[i];
+    F64 b = g_grid.points[i + 1];
+    F64 prefac = 0.5 * (b - a);
+
+    // Only integrate non-zero intervals
+    if (prefac > 1e-16) {
+      for (U32 gq_i = 0; gq_i < g_gauss_legendre.order; gq_i++) {
+        F64 w = g_gauss_legendre.weights[gq_i];
+        F64 z = g_gauss_legendre.abscissae[gq_i];
+        F64 r = (z * prefac) + ((a + b) * 0.5);
+        F64 term_prefac = (prefac * w);
+        F64 wf_at_r = 0.0;
+
+        for (U64 j = 0; j < coeff_size; j++) {
+          wf_at_r += coeffs[j] * compute_bspline_F64(r, j + 1);
+        }
+
+        norm += term_prefac * wf_at_r * wf_at_r;
+      }
+    }
+  }
+  String8 out =
+      str8_pushf(scratch.arena, "n:%i, l:%i norm: %.2f \n", n, l, norm);
+  LOG(out.str);
+  scratch_release(scratch);
+}
 
-/////////////////
-//~
-// Main entry point
 function void EntryPoint(void) {
 
-  // Init subsystems
   OS_InitReceipt os_receipt = OS_init();
   OS_InitGfxReceipt os_gfx_receipt = OS_gfx_init(os_receipt);
 
-  // Main program
-  //hf_main();
-  hf_coulomb_integrals_test();
+  g_filename_arena = m_make_arena_reserve(Megabytes(2));
+  g_base_arena = m_make_arena();
+
+  set_up_gauss_legendre_points(g_base_arena);
+
+  //- 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,
+                         g_grid.num_steps);
+  write_array_F64(str8_lit(grid_file_path), g_grid.points, g_grid.num_steps,
+                  "%13.6e\n");
+
+  //- The BSpline context is the knotpoints and the BSpline order etc.
+  set_up_bspline_context(g_base_arena);
+  write_array_F64(str8_lit(knotpoints_file_path), g_bspline_ctx.knotpoints,
+                  g_bspline_ctx.num_knotpoints, "%13.6e\n");
+
+  //- Then we generate the BSplines and save them off for reference and
+  // debugging.
+  set_up_bsplines_at_points_and_write_matrix_F64(g_base_arena);
+
+  U32 N = g_bspline_ctx.num_knotpoints;
+  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);
+  // This will be the inverse of B, but to start with we construct B.
+  Mat_F64 B_inv = mat_F64(mat_size1, mat_size2);
+  // A is the actual matrix for each eigenvalue problem.
+  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,
+  // so we invert the B matrix and compute the product A = B^-1H before calling
+  // zgeev
+  mat_invert_F64(&B_inv);
+
+  // 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];
+    if (l > 1e-16) {
+      F64 l_factor = l * (l + 1.0);
+      U64 mat_size = H_l.size1 * H_l.size2;
+      mat_F64_copy_to_dst(&H_l, &H_l_base);
+      // Multiply l(l+1)
+      cblas_dscal(mat_size, l_factor, H_l.data, 1);
+      // Add H = H_base + H_l
+      cblas_daxpy(mat_size, 1.0, H_l.data, 1, H.data, 1);
+    }
+
+    // Multiply to get A = B^-1 H
+    {
+      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);
+    }
+
+    // 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
+    {
+      S32 size1 = A.size1;
+      S32 lda = size1;
+      S32 ldvl = size1;
+      S32 ldvr = size1;
+      S32 info;
+      S32 lwork;
+
+      F64 wkopt;
+      F64 *work;
+
+      F64 *wr = PushArray(mkl_arena, F64, size1);
+      F64 *wi = PushArray(mkl_arena, F64, size1);
+      F64 *vl = PushArray(mkl_arena, F64, ldvl * size1);
+      F64 *vr = PushArray(mkl_arena, F64, ldvr * size1);
+
+      lwork = -1;
+      F64 *a = A.data;
+      dgeev("Vectors", "Vectors", &size1, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr,
+            &wkopt, &lwork, &info);
+      lwork = (S32)wkopt;
+      work = (F64 *)malloc(lwork * sizeof(F64));
+      dgeev("Vectors", "Vectors", &size1, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr,
+            work, &lwork, &info);
+      if (info > 0) {
+        LOG("Failed to compute eigenvalues in dgeev\n");
+        exit(1);
+      }
+
+      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( "Eigenvalues sorted: ", size1, wr, wi );
+
+      U32 i = 0;
+      F64 energy = -1000.0;
+      U32 counter = 0;
+      while (energy < 0.0) {
+        energy = wr[i];
+        U64 energy_index = sorted_indices[i];
+        U64 n = 1 + i;
+        if (ang_mom_idx > 0) {
+          n = 2 + i;
+        }
+
+        // compute_wf_norm_F64(eigensolution.right_eigenvectors.matrix[energy_index],
+        // size1, n, ang_mom_idx);
+
+        write_array_F64(
+            get_eigenvector_filename(g_filename_arena, n, ang_mom_idx),
+            eigensolution.right_eigenvectors.matrix[energy_index], size1,
+            "%13.6e\n");
+
+        i += 1;
+        counter += 1;
+        if (counter > 10) {
+          break;
+        }
+      }
+
+      free((void *)work);
+    }
+  }
 }

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;
+        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;
-        }
+        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(;;)
-      {
-        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))
+        U8 *ptr = node->string.str;
+        U8 *opl = ptr + node->string.size;
+        for(;;)
         {
-            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:;
 }