diff --git a/src/main.c b/src/main.c index ddd7d32..2e87633 100644 --- a/src/main.c +++ b/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); - String8List list = {0}; - for(U32 i = 0; i < array_size; i++) - { - str8_list_pushf(scratch.arena, &list, fmt, values[i]); - } - write_string_list_to_file(scratch.arena, path_to_file, &list); + ArenaTemp scratch = scratch_get(0, 0); + String8List list = {0}; + for(U32 i = 0; i < array_size; i++) + { + str8_list_pushf(scratch.arena, &list, fmt, values[i]); + } + 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 EntryPoint(void) +function void bspline_things() { - 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(); + +} diff --git a/src/test_mkl.c b/src/test_mkl.c index 8698c54..41a9bbb 100644 --- a/src/test_mkl.c +++ b/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 -#define TEST_MKL_N 4 +#ifndef TEST_ZGEEV_N +#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]; - MKL_Complex16 w[TEST_MKL_N], vl[TEST_MKL_N*TEST_MKL_N], vr[TEST_MKL_N*TEST_MKL_N]; - MKL_Complex16 a[TEST_MKL_N*TEST_MKL_N] = { + double rwork[2*TEST_ZGEEV_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_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 LDA N +#define TEST_DSYEVD_N 5 +#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 a[LDA*N] = { + double w[TEST_DSYEVD_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, diff --git a/timeBuild.ctm b/timeBuild.ctm index 93129fc..f03b3da 100644 Binary files a/timeBuild.ctm and b/timeBuild.ctm differ