pulled out intel mkl zgeev test

src/main.c

@@ -21,7 +21,7 @@
 #include "base/base_inc.c"
 #include "os/os_inc.c"
 #include "os/os_entry_point.c"
-#include "test_mkl.c"
+//#include "test_mkl.c"
 
 #define grid_file_path_bin "D:\\dev\\eigsol_gpu\\out\\grid.bin"
 #define grid_file_path "D:\\dev\\eigsol_gpu\\out\\grid.dat"
@@ -58,6 +58,23 @@ struct BSplineCtx
   F64 *bsplines;
 };
 
+typedef struct Orbital Orbital;
+struct Orbital 
+{
+  U32 n;
+  U32 l;
+  U32 j;
+};
+
+
+typedef struct Atom Atom;
+struct Atom 
+{
+  U32 N;
+  Orbital *orbitals;
+};
+
+
 //~ Globals
 global Grid g_grid = {0};
 global BSplineCtx g_bspline_ctx = {0};
@@ -75,7 +92,7 @@ write_array_binary_F64(String8 path_to_file, F64 *values, U32 array_size)
     temp.size = sizeof(F64)*array_size;
     str8_list_push(scratch.arena, &list, temp);
     OS_file_write(scratch.arena, file_handle, 0, list, 0);
-    
+
     String8List log_list = {0};
     str8_list_push(scratch.arena, &log_list, str8_lit("Wrote binary array data to"));
     str8_list_push(scratch.arena, &log_list, path_to_file);
@@ -112,19 +129,19 @@ write_string_list_to_file(Arena *arena, String8 path, String8List *list)
 }
 
 
-  function void
+function void
 write_array_F64(String8 path_to_file, F64 *values, U32 array_size, char* fmt) 
 {
-    ArenaTemp scratch = scratch_get(0, 0);
+  ArenaTemp scratch = scratch_get(0, 0);
-    String8List list = {0};
+  String8List list = {0};
-    for(U32 i = 0; i < array_size; i++) 
+  for(U32 i = 0; i < array_size; i++) 
-    {
+  {
-      str8_list_pushf(scratch.arena, &list, fmt, values[i]);
+    str8_list_pushf(scratch.arena, &list, fmt, values[i]);
-    }
+  }
-    write_string_list_to_file(scratch.arena, path_to_file, &list);
+  write_string_list_to_file(scratch.arena, path_to_file, &list);
 }
 
-  function F64 
+function F64 
 bspline_recursion(F64 x, U32 k, U32 i)
 {
   F64 *t = g_bspline_ctx.knotpoints;
@@ -165,7 +182,7 @@ bspline_recursion(F64 x, U32 k, U32 i)
 }
 
 
-  function F64 
+function F64 
 get_bspline_F64(F64 x_coord, U32 index)
 {
   U32 k = g_bspline_ctx.order;
@@ -173,7 +190,7 @@ get_bspline_F64(F64 x_coord, U32 index)
   return out;
 }
 
-  function void
+function void
 set_up_grid(Arena *arena)
 {
   g_grid.start = 0.0;
@@ -192,7 +209,7 @@ set_up_grid(Arena *arena)
 }
 
 
-  function void
+function void
 set_up_bspline_context(Arena* arena)
 {
   // Create knotpoint sequence.
@@ -225,7 +242,7 @@ set_up_bspline_context(Arena* arena)
 }
 
 
-  function void
+function void
 write_bsplines_to_matrix_F64(Arena *arena)
 {
 
@@ -304,13 +321,8 @@ write_bsplines_to_matrix_F64(Arena *arena)
 
 }
 
-
+function void bspline_things() 
-function
-void EntryPoint(void) 
 {
-  OS_InitReceipt os_receipt = OS_init();
-  OS_InitGfxReceipt os_gfx_receipt = OS_gfx_init(os_receipt);
-
   Arena *arena = m_make_arena();
 
   //- Set up grid and write to file.
@@ -327,9 +339,101 @@ void EntryPoint(void)
   //- Then we generate the BSplines and save them off for reference and debugging.
   write_bsplines_to_matrix_F64(arena);
 
-  //- Test MKL
+}
-  test_mkl_zgeev();
+
-  test_mkl_dsyevd();
+
+
+/* Auxiliary routine: printing a matrix */
+function void
+print_matrix_Z64( char* desc, int m, int n, Z64* a, int lda ) {
+  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( i = 0; i < m; i++ ) {
+    for( j = 0; j < n; j++ ) {
+      String8 outstr = str8_pushf(scratch.arena, " (%6.2f,%6.2f)", a[i+j*lda].re, a[i+j*lda].im );
+      LOG(outstr.str);
+      //printf(" (%6.2f,%6.2f)", a[i+j*lda].real, a[i+j*lda].imag );
+    }
+    LOG(newline.str);
+    //printf("\n");
+  }
+}
+
+
+
+function void EntryPoint(void) 
+{
+  OS_InitReceipt os_receipt = OS_init();
+  OS_InitGfxReceipt os_gfx_receipt = OS_gfx_init(os_receipt);
+
+  Arena *arena = m_make_arena();
+
+  LOG("\n\n---- Calling Intel MKL zgeev test (Using Z64 instead of MKL_Complex16 etc) ---- \n\n");
+
+  {
+    U32 N = 4;
+    S32 n = N, lda = N, ldvl = N, ldvr = N, info, lwork;
+    Z64 wkopt;
+    Z64 *work;
+
+    F64 *rwork = PushArray(arena, F64, 2*N);
+    Z64 *w = PushArray(arena, Z64, N);
+    Z64 *vl = PushArray(arena, Z64, N*N);
+    Z64 *vr = PushArray(arena, Z64, N*N);
+    Z64 *a = PushArray(arena, Z64, N*N);
+    a[0]  = (Z64){-3.84,  2.25};
+    a[1]  = (Z64){-0.66,  0.83}; 
+    a[2]  = (Z64){-3.99, -4.73};
+    a[3]  = (Z64){ 7.74,  4.18};
+    a[4]  = (Z64){-8.94, -4.75};
+    a[5]  = (Z64){-4.40, -3.82};
+    a[6]  = (Z64){-5.88, -6.60};
+    a[7]  = (Z64){ 3.66, -7.53};
+    a[8]  = (Z64){ 8.95, -6.53};
+    a[9]  = (Z64){-3.50, -4.26};
+    a[10] = (Z64){-3.36, -0.40};
+    a[11] = (Z64){ 2.58,  3.60};
+    a[12] = (Z64){-9.87,  4.82};
+    a[13] = (Z64){-3.15,  7.36};
+    a[14] = (Z64){-0.75,  5.23};
+    a[15] = (Z64){ 4.59,  5.41};
+
+    /* Executable statements */
+    LOG( " ZGEEV Example Program Results\n" );
+    /* Query and allocate the optimal workspace */
+    lwork = -1;
+    zgeev( "Vectors", "Vectors", &n, a, &lda, w, vl, &ldvl, vr, &ldvr,
+        &wkopt, &lwork, rwork, &info );
+    lwork = (S32)wkopt.re;
+    work = (Z64*)malloc( lwork*sizeof(Z64) );
+    /* Solve eigenproblem */
+    zgeev( "Vectors", "Vectors", &n, a, &lda, w, vl, &ldvl, vr, &ldvr,
+        work, &lwork, rwork, &info );
+    /* Check for convergence */
+    if( info > 0 ) {
+      LOG( "The algorithm failed to compute eigenvalues.\n" );
+      exit( 1 );
+    }
+    /* Print eigenvalues */
+    print_matrix_Z64( "Eigenvalues", 1, n, w, 1 );
+    /* Print left eigenvectors */
+    print_matrix_Z64( "Left eigenvectors", n, n, vl, ldvl );
+    /* Print right eigenvectors */
+    print_matrix_Z64( "Right eigenvectors", n, n, vr, ldvr );
+    /* Free workspace */
+    free( (void*)work );
+    //exit( 0 );
+  } /* End of ZGEEV Example */
+
+
+
+
+
+
 
   LOG("\n\n---- Skipping CUDA part for now ---- \n\n");
   //cuda_entry_point();
@@ -341,4 +445,11 @@ void EntryPoint(void)
 }
 
 
+  function void 
+testing_MKL()
+{
 
+  test_mkl_zgeev();
+  test_mkl_dsyevd();
+
+}

src/test_mkl.c

@@ -70,23 +70,23 @@
 function void 
 print_matrix_cmplx16( char* desc, int m, int n, MKL_Complex16* a, int lda );
 
-#ifndef TEST_MKL_N
+#ifndef TEST_ZGEEV_N
-#define TEST_MKL_N 4
+#define TEST_ZGEEV_N 4
 #endif
 
 /* Main program */
 static void
 test_mkl_zgeev(void) {
   /* Locals */
-  int N = TEST_MKL_N;
+  int N = TEST_ZGEEV_N;
   MKL_INT n = N, lda = N, ldvl = N, ldvr = N, info, lwork;
   MKL_Complex16 wkopt;
   MKL_Complex16* work;
   /* Local arrays */
   /* rwork dimension should be at least 2*n */
-  double rwork[2*TEST_MKL_N];
+  double rwork[2*TEST_ZGEEV_N];
-  MKL_Complex16 w[TEST_MKL_N], vl[TEST_MKL_N*TEST_MKL_N], vr[TEST_MKL_N*TEST_MKL_N];
+  MKL_Complex16 w[TEST_ZGEEV_N], vl[TEST_ZGEEV_N*TEST_ZGEEV_N], vr[TEST_ZGEEV_N*TEST_ZGEEV_N];
-  MKL_Complex16 a[TEST_MKL_N*TEST_MKL_N] = {
+  MKL_Complex16 a[TEST_ZGEEV_N*TEST_ZGEEV_N] = {
     {-3.84,  2.25}, {-0.66,  0.83}, {-3.99, -4.73}, { 7.74,  4.18},
     {-8.94, -4.75}, {-4.40, -3.82}, {-5.88, -6.60}, { 3.66, -7.53},
     { 8.95, -6.53}, {-3.50, -4.26}, {-3.36, -0.40}, { 2.58,  3.60},
@@ -202,21 +202,21 @@ print_matrix_cmplx16( char* desc, int m, int n, MKL_Complex16* a, int lda ) {
 extern void print_matrix( char* desc, int m, int n, double* a, int lda );
 
 /* Parameters */
-#define N 5
+#define TEST_DSYEVD_N 5
-#define LDA N
+#define LDA TEST_DSYEVD_N
 
 /* Main program */
 function void 
 test_mkl_dsyevd(void) {
   /* Locals */
-  MKL_INT n = N, lda = LDA, info, lwork, liwork;
+  MKL_INT n = TEST_DSYEVD_N, lda = LDA, info, lwork, liwork;
   MKL_INT iwkopt;
   MKL_INT* iwork;
   double wkopt;
   double* work;
   /* Local arrays */
-  double w[N];
+  double w[TEST_DSYEVD_N];
-  double a[LDA*N] = {
+  double a[LDA*TEST_DSYEVD_N] = {
     6.39,  0.00,  0.00,  0.00,  0.00,
     0.13,  8.37,  0.00,  0.00,  0.00,
     -8.23, -4.46, -9.58,  0.00,  0.00,