Compare commits
No commits in common. "d670d018b714dd5989b260c2f62bfe2659aadef3" and "fc2c897727799dad746db03c9a97c00823fecba5" have entirely different histories.
d670d018b7
...
fc2c897727
40
build.bat
40
build.bat
@ -1,17 +1,29 @@
|
||||
@echo off
|
||||
set CommonCompilerFlags=/nologo /Zi /FC /Od
|
||||
|
||||
set mkl_root=E:/lib/intel_mkl/mkl/2025.3
|
||||
set mkl_core=%mkl_root%/lib/mkl_core.lib
|
||||
set mkl_intel_lp64=%mkl_root%/lib/mkl_intel_lp64.lib
|
||||
set mkl_intel_thread=%mkl_root%/lib/mkl_intel_thread.lib
|
||||
set MKLCOMPILER=E:/lib/intel_mkl/compiler/2025.3
|
||||
set libiomp5md=%MKLCOMPILER%/lib/libiomp5md.lib
|
||||
set libiompdll_path=E:\lib\intel_mkl\compiler\2025.3\bin
|
||||
ctime -begin timeBuild.ctm
|
||||
|
||||
@rem /WX /W4 /wd4201 /wd4100 /wd4189 /wd4244 /wd4127 /wd4456
|
||||
@rem set CommonCompilerFlags="/nologo /Zi /FC"
|
||||
set CommonCompilerFlags=/nologo /Zi /FC /Od
|
||||
@rem /WX /W4 /wd4201 /wd4100 /wd4189 /wd4244 /wd4127 /wd4456
|
||||
@rem
|
||||
|
||||
|
||||
set mkl_root=D:/lib/oneAPI_mkl/mkl/2021.3.0
|
||||
set mkl_core=%mkl_root%/lib/intel64/mkl_core.lib
|
||||
set mkl_intel_lp64=%mkl_root%/lib/intel64/mkl_intel_lp64.lib
|
||||
set mkl_intel_thread=%mkl_root%/lib/intel64/mkl_intel_thread.lib
|
||||
set MKLCOMPILER=D:/lib/oneAPI_mkl/compiler/2021.3.0/windows/compiler
|
||||
set libiomp5md=%MKLCOMPILER%/lib/intel64_win/libiomp5md.lib
|
||||
|
||||
set libiompdll_path=D:\lib\oneAPI_mkl\compiler\2021.3.0\windows\redist\intel64_win\compiler
|
||||
set libiompdll_name=libiomp5md.dll
|
||||
set libiompdll=%libiompdll_path%\%libiompdll_name%
|
||||
|
||||
set Sources=../src/main.c
|
||||
|
||||
|
||||
|
||||
IF NOT EXIST .\out mkdir .\out
|
||||
IF NOT EXIST .\build mkdir .\build
|
||||
pushd .\build
|
||||
@ -24,19 +36,15 @@ if not exist "%libiompdll_name%" (
|
||||
) else (
|
||||
echo Copied openmp dll: %libiompdll_name%
|
||||
)
|
||||
|
||||
)
|
||||
|
||||
cl %CommonCompilerFlags% %Sources% /I"%mkl_root%\include" ^
|
||||
/link %mkl_core% %mkl_intel_lp64% %mkl_intel_thread% %libiomp5md%
|
||||
cl %CommonCompilerFlags% %Sources% /I"%mkl_root%\include" /link %mkl_core% %mkl_intel_lp64% %mkl_intel_thread% %libiomp5md%
|
||||
|
||||
set LastError=%ERRORLEVEL%
|
||||
popd
|
||||
|
||||
|
||||
echo Build complete
|
||||
|
||||
set PATH=%PATH%;%mkl_root%\bin
|
||||
set PATH=%PATH%;%MKLCOMPILER%\bin
|
||||
set LastError=%ERRORLEVEL%
|
||||
ctime -end timeBuild.ctm %LastError%
|
||||
IF NOT %LastError%==0 GOTO :end
|
||||
|
||||
:end
|
||||
|
||||
File diff suppressed because one or more lines are too long
@ -1,34 +1,34 @@
|
||||
4.079148e-01
|
||||
6.284958e-01
|
||||
5.853385e-01
|
||||
2.830382e-01
|
||||
1.166302e-01
|
||||
4.437820e-02
|
||||
1.608212e-02
|
||||
5.645430e-03
|
||||
1.935000e-03
|
||||
6.523220e-04
|
||||
2.166921e-04
|
||||
7.130839e-05
|
||||
2.321027e-05
|
||||
7.519194e-06
|
||||
2.411514e-06
|
||||
7.732122e-07
|
||||
2.449929e-07
|
||||
7.808282e-08
|
||||
2.445164e-08
|
||||
7.783642e-09
|
||||
2.403077e-09
|
||||
7.695295e-10
|
||||
2.325414e-10
|
||||
7.588069e-11
|
||||
2.207375e-11
|
||||
7.527071e-12
|
||||
2.032901e-12
|
||||
7.619167e-13
|
||||
1.762993e-13
|
||||
8.066878e-14
|
||||
1.303573e-14
|
||||
9.445442e-15
|
||||
-2.760762e-16
|
||||
1.847288e-15
|
||||
-4.079148e-01
|
||||
-6.284958e-01
|
||||
-5.853385e-01
|
||||
-2.830382e-01
|
||||
-1.166302e-01
|
||||
-4.437820e-02
|
||||
-1.608212e-02
|
||||
-5.645430e-03
|
||||
-1.935000e-03
|
||||
-6.523220e-04
|
||||
-2.166921e-04
|
||||
-7.130839e-05
|
||||
-2.321027e-05
|
||||
-7.519194e-06
|
||||
-2.411514e-06
|
||||
-7.732122e-07
|
||||
-2.449929e-07
|
||||
-7.808282e-08
|
||||
-2.445164e-08
|
||||
-7.783642e-09
|
||||
-2.403077e-09
|
||||
-7.695295e-10
|
||||
-2.325412e-10
|
||||
-7.588092e-11
|
||||
-2.207390e-11
|
||||
-7.526818e-12
|
||||
-2.032911e-12
|
||||
-7.621029e-13
|
||||
-1.765319e-13
|
||||
-8.079409e-14
|
||||
-1.311978e-14
|
||||
-9.257067e-15
|
||||
1.842392e-16
|
||||
-1.701083e-15
|
||||
|
||||
@ -1,34 +1,34 @@
|
||||
1.067985e-01
|
||||
9.343962e-02
|
||||
-1.122171e-01
|
||||
-3.647821e-01
|
||||
-4.769945e-01
|
||||
-4.732031e-01
|
||||
-4.077962e-01
|
||||
-3.223604e-01
|
||||
-2.403128e-01
|
||||
-1.716945e-01
|
||||
-1.187758e-01
|
||||
-8.011225e-02
|
||||
-5.294166e-02
|
||||
-3.440233e-02
|
||||
-2.204206e-02
|
||||
-1.395428e-02
|
||||
-8.743348e-03
|
||||
-5.429335e-03
|
||||
-3.344953e-03
|
||||
-2.046453e-03
|
||||
-1.244260e-03
|
||||
-7.523093e-04
|
||||
-4.525787e-04
|
||||
-2.710226e-04
|
||||
-1.616214e-04
|
||||
-9.600751e-05
|
||||
-5.681871e-05
|
||||
-3.349503e-05
|
||||
-1.964742e-05
|
||||
-1.142544e-05
|
||||
-6.511164e-06
|
||||
-3.501003e-06
|
||||
-1.380734e-06
|
||||
-4.297736e-07
|
||||
-1.067985e-01
|
||||
-9.343962e-02
|
||||
1.122171e-01
|
||||
3.647821e-01
|
||||
4.769945e-01
|
||||
4.732031e-01
|
||||
4.077962e-01
|
||||
3.223604e-01
|
||||
2.403128e-01
|
||||
1.716945e-01
|
||||
1.187758e-01
|
||||
8.011225e-02
|
||||
5.294166e-02
|
||||
3.440233e-02
|
||||
2.204206e-02
|
||||
1.395428e-02
|
||||
8.743348e-03
|
||||
5.429335e-03
|
||||
3.344953e-03
|
||||
2.046453e-03
|
||||
1.244260e-03
|
||||
7.523093e-04
|
||||
4.525787e-04
|
||||
2.710226e-04
|
||||
1.616214e-04
|
||||
9.600751e-05
|
||||
5.681871e-05
|
||||
3.349503e-05
|
||||
1.964742e-05
|
||||
1.142544e-05
|
||||
6.511164e-06
|
||||
3.501003e-06
|
||||
1.380734e-06
|
||||
4.297736e-07
|
||||
|
||||
@ -1,34 +1,34 @@
|
||||
-2.548436e-03
|
||||
-1.768047e-02
|
||||
-6.801789e-02
|
||||
-1.434030e-01
|
||||
-2.235356e-01
|
||||
-2.907379e-01
|
||||
-3.357788e-01
|
||||
-3.567533e-01
|
||||
-3.563788e-01
|
||||
-3.395486e-01
|
||||
-3.116223e-01
|
||||
-2.774422e-01
|
||||
-2.408901e-01
|
||||
-2.047894e-01
|
||||
-1.709988e-01
|
||||
-1.405887e-01
|
||||
-1.140369e-01
|
||||
-9.140819e-02
|
||||
-7.250238e-02
|
||||
-5.696782e-02
|
||||
-4.438296e-02
|
||||
-3.431126e-02
|
||||
-2.633551e-02
|
||||
-2.007697e-02
|
||||
-1.520409e-02
|
||||
-1.143438e-02
|
||||
-8.531931e-03
|
||||
-6.302777e-03
|
||||
-4.589046e-03
|
||||
-3.262862e-03
|
||||
-2.220364e-03
|
||||
-1.376048e-03
|
||||
-5.929916e-04
|
||||
-1.933374e-04
|
||||
2.548436e-03
|
||||
1.768047e-02
|
||||
6.801789e-02
|
||||
1.434030e-01
|
||||
2.235356e-01
|
||||
2.907379e-01
|
||||
3.357788e-01
|
||||
3.567533e-01
|
||||
3.563788e-01
|
||||
3.395486e-01
|
||||
3.116223e-01
|
||||
2.774422e-01
|
||||
2.408901e-01
|
||||
2.047894e-01
|
||||
1.709988e-01
|
||||
1.405887e-01
|
||||
1.140369e-01
|
||||
9.140819e-02
|
||||
7.250238e-02
|
||||
5.696782e-02
|
||||
4.438296e-02
|
||||
3.431126e-02
|
||||
2.633551e-02
|
||||
2.007697e-02
|
||||
1.520409e-02
|
||||
1.143438e-02
|
||||
8.531931e-03
|
||||
6.302777e-03
|
||||
4.589046e-03
|
||||
3.262862e-03
|
||||
2.220364e-03
|
||||
1.376048e-03
|
||||
5.929916e-04
|
||||
1.933374e-04
|
||||
|
||||
34
out/eigvec_n6_l2.dat
Normal file
34
out/eigvec_n6_l2.dat
Normal file
@ -0,0 +1,34 @@
|
||||
-5.741527e-03
|
||||
-3.569258e-02
|
||||
-1.169176e-01
|
||||
-1.897552e-01
|
||||
-1.957073e-01
|
||||
-1.207026e-01
|
||||
4.698912e-03
|
||||
1.300312e-01
|
||||
2.099986e-01
|
||||
2.200159e-01
|
||||
1.605827e-01
|
||||
5.287137e-02
|
||||
-7.055343e-02
|
||||
-1.760989e-01
|
||||
-2.377474e-01
|
||||
-2.422956e-01
|
||||
-1.907935e-01
|
||||
-9.665084e-02
|
||||
1.856075e-02
|
||||
1.303638e-01
|
||||
2.165004e-01
|
||||
2.609945e-01
|
||||
2.565772e-01
|
||||
2.052819e-01
|
||||
1.173563e-01
|
||||
8.888416e-03
|
||||
-1.013070e-01
|
||||
-1.949579e-01
|
||||
-2.572314e-01
|
||||
-2.788515e-01
|
||||
-2.571990e-01
|
||||
-1.963299e-01
|
||||
-9.506447e-02
|
||||
-3.283758e-02
|
||||
34
out/eigvec_n7_l0.dat
Normal file
34
out/eigvec_n7_l0.dat
Normal file
@ -0,0 +1,34 @@
|
||||
4.415584e-02
|
||||
-2.080198e-03
|
||||
-1.208549e-01
|
||||
-1.232666e-01
|
||||
1.381534e-02
|
||||
1.544869e-01
|
||||
1.992137e-01
|
||||
1.305172e-01
|
||||
-6.669893e-03
|
||||
-1.442392e-01
|
||||
-2.249528e-01
|
||||
-2.214819e-01
|
||||
-1.397799e-01
|
||||
-1.097338e-02
|
||||
1.224366e-01
|
||||
2.207053e-01
|
||||
2.574974e-01
|
||||
2.251569e-01
|
||||
1.344371e-01
|
||||
9.825503e-03
|
||||
-1.176775e-01
|
||||
-2.183055e-01
|
||||
-2.700506e-01
|
||||
-2.627792e-01
|
||||
-1.994307e-01
|
||||
-9.443627e-02
|
||||
2.999191e-02
|
||||
1.487476e-01
|
||||
2.388476e-01
|
||||
2.836426e-01
|
||||
2.755399e-01
|
||||
2.168929e-01
|
||||
1.066895e-01
|
||||
3.713090e-02
|
||||
34
out/eigvec_n7_l1.dat
Normal file
34
out/eigvec_n7_l1.dat
Normal file
@ -0,0 +1,34 @@
|
||||
3.344934e-02
|
||||
1.015861e-01
|
||||
1.437945e-01
|
||||
3.591032e-02
|
||||
-1.135295e-01
|
||||
-1.940374e-01
|
||||
-1.636417e-01
|
||||
-4.787872e-02
|
||||
9.288932e-02
|
||||
1.980292e-01
|
||||
2.296506e-01
|
||||
1.809186e-01
|
||||
7.223784e-02
|
||||
-5.977311e-02
|
||||
-1.754378e-01
|
||||
-2.433252e-01
|
||||
-2.473296e-01
|
||||
-1.887683e-01
|
||||
-8.403244e-02
|
||||
4.084245e-02
|
||||
1.569740e-01
|
||||
2.391540e-01
|
||||
2.708134e-01
|
||||
2.466788e-01
|
||||
1.729256e-01
|
||||
6.512299e-02
|
||||
-5.537940e-02
|
||||
-1.657922e-01
|
||||
-2.461079e-01
|
||||
-2.824849e-01
|
||||
-2.693014e-01
|
||||
-2.096634e-01
|
||||
-1.025569e-01
|
||||
-3.559765e-02
|
||||
@ -1,8 +1,9 @@
|
||||
#ifndef BASE_TYPES_H
|
||||
#define BASE_TYPES_H
|
||||
|
||||
#include <math.h>
|
||||
#include <stdint.h>
|
||||
#include <stdlib.h>
|
||||
#include <stdio.h>
|
||||
|
||||
/////////////////////////
|
||||
//~ Macros
|
||||
@ -12,17 +13,17 @@
|
||||
|
||||
// TODO(anton): Understand this, yoinked from rjf's layer.
|
||||
#if LANG_CPP
|
||||
#define no_name_mangle extern "C"
|
||||
# define no_name_mangle extern "C"
|
||||
#else
|
||||
#define no_name_mangle
|
||||
# define no_name_mangle
|
||||
#endif
|
||||
|
||||
// TODO(anton): OS_WINDOWS dll import/export macros
|
||||
|
||||
/////////////////////////
|
||||
//- Keywords
|
||||
// Static is stupid and means different things depending on context in C and
|
||||
// C++. These defines increases readability.
|
||||
// Static is stupid and means different things depending on context in C and C++.
|
||||
// These defines increases readability.
|
||||
#define function static // Function internal to compilation unit.
|
||||
#define local_persist static
|
||||
#define global static
|
||||
@ -30,123 +31,111 @@
|
||||
|
||||
// TODO(anton): Understand and add good comment on this.
|
||||
#if LANG_CPP
|
||||
#define root_global no_name_mangle
|
||||
#define root_function function
|
||||
# define root_global no_name_mangle
|
||||
# define root_function function
|
||||
#else
|
||||
#define root_global extern
|
||||
#define root_function function
|
||||
# define root_global extern
|
||||
# define root_function function
|
||||
#endif
|
||||
|
||||
#define inline_function inline static
|
||||
|
||||
#if OS_WINDOWS
|
||||
#pragma section(".roglob", read)
|
||||
#define read_only __declspec(allocate(".roglob"))
|
||||
# pragma section(".roglob", read)
|
||||
# define read_only __declspec(allocate(".roglob"))
|
||||
#else
|
||||
#define read_only
|
||||
# define read_only
|
||||
#endif
|
||||
|
||||
#if COMPILER_MSVC
|
||||
#define per_thread __declspec(thread)
|
||||
# define per_thread __declspec(thread)
|
||||
#else
|
||||
#error Thread keyword not abstracted on compiler.
|
||||
# error Thread keyword not abstracted on compiler.
|
||||
#endif
|
||||
|
||||
/////////////////////////
|
||||
//- Memory operations
|
||||
// It's nice to put these in macros, so we can swap out the functionality from
|
||||
// standard library, eventually.
|
||||
// It's nice to put these in macros, so we can swap out the functionality from standard library, eventually.
|
||||
#define MemoryCopy memcpy
|
||||
#define MemoryMove memmove
|
||||
#define MemorySet memset
|
||||
#define MemorySet memset
|
||||
|
||||
// NOTE(anton): This gives a 4127 compiler warning for the sizeof conditional.
|
||||
// This should be ignored
|
||||
#define MemoryCopyStruct(dst, src) \
|
||||
do { \
|
||||
Assert(sizeof(*(dst)) == sizeof(*(src))); \
|
||||
MemoryCopy((dst), (src), sizeof(*(dst))); \
|
||||
} while (0)
|
||||
// NOTE(anton): This gives a 4127 compiler warning for the sizeof conditional. This should be ignored
|
||||
#define MemoryCopyStruct(dst, src) do { Assert(sizeof(*(dst)) == sizeof(*(src))); MemoryCopy((dst), (src), sizeof(*(dst))); } while(0)
|
||||
|
||||
#define MemoryZero(ptr, size) MemorySet((ptr), 0, (size))
|
||||
#define MemoryZeroStruct(ptr) MemoryZero((ptr), sizeof(*(ptr)))
|
||||
#define MemoryZeroArray(arr) MemoryZero((arr), sizeof(arr))
|
||||
#define MemoryZeroArray(arr) MemoryZero((arr), sizeof(arr))
|
||||
|
||||
/////////////////////////
|
||||
//- Integer/pointer/array/type manipulations
|
||||
|
||||
#define ArrayCount(a) (sizeof(a) / sizeof((a)[0]))
|
||||
#define IntFromPtr(p) (U64)(((U8 *)p) - 0)
|
||||
#define PtrFromInt(i) (void *)(((U8 *)0) + i)
|
||||
#define Member(type, member_name) ((type *)0)->member_name
|
||||
#define IntFromPtr(p) (U64)(((U8*)p) - 0)
|
||||
#define PtrFromInt(i) (void*)(((U8*)0) + i)
|
||||
#define Member(type, member_name) ((type *)0)->member_name
|
||||
// TODO(anton): Understand why this becomes offset actually
|
||||
#define OffsetOf(type, member_name) IntFromPtr(&Member(type, member_name))
|
||||
// TODO(anton): Understand this
|
||||
#define BaseFromMember(type, member_name, ptr) \
|
||||
(type *)((U8 *)(ptr) - OffsetOf(type, member_name))
|
||||
#define BaseFromMember(type, member_name, ptr) (type *)((U8 *)(ptr) - OffsetOf(type, member_name))
|
||||
|
||||
#define Bytes(n) (n)
|
||||
#define Bytes(n) (n)
|
||||
#define Kilobytes(n) (n << 10) // 2^10 == 1024 etc
|
||||
#define Megabytes(n) (n << 20)
|
||||
#define Gigabytes(n) (((U64)n) << 30)
|
||||
#define Terabytes(n) (((U64)n) << 40)
|
||||
|
||||
#define Thousand(n) ((n) * 1000)
|
||||
#define Million(n) ((n) * 1000000)
|
||||
#define Billion(n) ((n) * 1000000000LL)
|
||||
#define Thousand(n) ((n)*1000)
|
||||
#define Million(n) ((n)*1000000)
|
||||
#define Billion(n) ((n)*1000000000LL)
|
||||
|
||||
#define AbsoluteValueU64(x) (U64) llabs((U64)(x))
|
||||
#define AbsoluteValueU64(x) (U64)llabs((U64)(x))
|
||||
|
||||
/////////////////////////
|
||||
//- Linked list helpers
|
||||
|
||||
#define CheckNull(p) ((p) == 0)
|
||||
#define SetNull(p) ((p) = 0)
|
||||
#define CheckNull(p) ((p)==0)
|
||||
#define SetNull(p) ((p)=0)
|
||||
// Link list helper macros that are a bit involved
|
||||
|
||||
// Suffixes N,P,Z means that we have (N)ext, (P)rev arguments and/or a (Z)ero
|
||||
// check and/or set argument f, l, n are "first", "last", "node" I think? DLL
|
||||
// Doubly Linked List: Each node has a prev and next pointer. Operations: Push
|
||||
// back, Push front, remove
|
||||
#define DLLInsert_NPZ(f, l, p, n, next, prev, zchk, zset) \
|
||||
(zchk(f) ? (((f) = (l) = (n)), zset((n)->next), zset((n)->prev)) \
|
||||
: zchk(p) ? (zset((n)->prev), (n)->next = (f), \
|
||||
(zchk(f) ? (0) : ((f)->prev = (n))), (f) = (n)) \
|
||||
: ((zchk((p)->next) ? (0) : (((p)->next->prev) = (n))), \
|
||||
(n)->next = (p)->next, (n)->prev = (p), (p)->next = (n), \
|
||||
((p) == (l) ? (l) = (n) : (0))))
|
||||
// Suffixes N,P,Z means that we have (N)ext, (P)rev arguments and/or a (Z)ero check and/or set argument
|
||||
// f, l, n are "first", "last", "node" I think?
|
||||
// DLL
|
||||
// Doubly Linked List: Each node has a prev and next pointer. Operations: Push back, Push front, remove
|
||||
#define DLLInsert_NPZ(f,l,p,n,next,prev,zchk,zset) \
|
||||
(zchk(f) ? (((f) = (l) = (n)), zset((n)->next), zset((n)->prev)) :\
|
||||
zchk(p) ? (zset((n)->prev), (n)->next = (f), (zchk(f) ? (0) : ((f)->prev = (n))), (f) = (n)) :\
|
||||
((zchk((p)->next) ? (0) : (((p)->next->prev) = (n))), (n)->next = (p)->next, (n)->prev = (p), (p)->next = (n),\
|
||||
((p) == (l) ? (l) = (n) : (0))))
|
||||
|
||||
#define DLLPushBack_NPZ(f, l, n, next, prev, zchk, zset) \
|
||||
DLLInsert_NPZ(f, l, l, n, next, prev, zchk, zset)
|
||||
#define DLLPushBack_NPZ(f,l,n,next,prev,zchk,zset) DLLInsert_NPZ(f,l,l,n,next,prev,zchk,zset)
|
||||
|
||||
#define DLLPushBack_NP(f, l, n, next, prev, zchk) \
|
||||
(zchk(f) ? ((f) = (l) = (n), (n)->next = (n)->prev = 0) \
|
||||
: ((n)->prev = (l), (l)->next = (n), (l) = (n), (n)->next = 0))
|
||||
#define DLLPushBack_NP(f, l, n, next, prev, zchk) \
|
||||
(zchk(f) ? ((f)=(l)=(n),(n)->next=(n)->prev=0) : ((n)->prev=(l),(l)->next=(n),(l)=(n),(n)->next=0))
|
||||
|
||||
// If f == n we put f to f->next, and f->prev = 0.
|
||||
// Else if l == n, we put l=l->prev, l->next = 0.
|
||||
// If l != n and f != n we set n->next->prev to n->prev, and n->prev->next to
|
||||
// n->next
|
||||
// If l != n and f != n we set n->next->prev to n->prev, and n->prev->next to n->next
|
||||
|
||||
#define DLLRemove_NP(f, l, n, next, prev) \
|
||||
(((f) == (n) ? ((f) = (f)->next, (f)->prev = 0) \
|
||||
: (l) == (n) \
|
||||
? ((l) = (l)->prev, (l)->next = 0) \
|
||||
: ((n)->next->prev = (n)->prev, (n)->prev->next = (n)->next)))
|
||||
#define DLLRemove_NP(f, l, n, next, prev) (((f) == (n) ? \
|
||||
((f)=(f)->next, (f)->prev=0) : \
|
||||
(l) == (n) ? \
|
||||
((l)=(l)->prev, (l)->next=0) : \
|
||||
((n)->next->prev=(n)->prev, \
|
||||
(n)->prev->next=(n)->next) ))
|
||||
|
||||
#define DLLRemove_NPZ(f, l, n, next, prev, zchk, zset) \
|
||||
(((f) == (n)) ? ((f) = (f)->next, (zchk(f) ? (zset(l)) : zset((f)->prev))) \
|
||||
: ((l) == (n)) ? ((l) = (l)->prev, (zchk(l) ? (zset(f)) : zset((l)->next))) \
|
||||
: ((zchk((n)->next) ? (0) : ((n)->next->prev = (n)->prev)), \
|
||||
(zchk((n)->prev) ? (0) : ((n)->prev->next = (n)->next))))
|
||||
#define DLLRemove_NPZ(f,l,n,next,prev,zchk,zset) (((f)==(n))?\
|
||||
((f)=(f)->next, (zchk(f) ? (zset(l)) : zset((f)->prev))):\
|
||||
((l)==(n))?\
|
||||
((l)=(l)->prev, (zchk(l) ? (zset(f)) : zset((l)->next))):\
|
||||
((zchk((n)->next) ? (0) : ((n)->next->prev=(n)->prev)),\
|
||||
(zchk((n)->prev) ? (0) : ((n)->prev->next=(n)->next))))
|
||||
|
||||
#define DLLPushBack(f, l, n) \
|
||||
DLLPushBack_NPZ(f, l, n, next, prev, CheckNull, SetNull)
|
||||
#define DLLPushBack(f, l, n) DLLPushBack_NPZ(f, l, n, next, prev, CheckNull, SetNull)
|
||||
// For front push I can just switch prev/next!
|
||||
#define DLLPushFront(f, l, n) \
|
||||
DLLPushBack_NPZ(l, f, n, prev, next, CheckNull, SetNull)
|
||||
#define DLLRemove(f, l, n) \
|
||||
DLLRemove_NPZ(f, l, n, next, prev, CheckNull, SetNull)
|
||||
#define DLLPushFront(f, l, n) DLLPushBack_NPZ(l, f, n, prev, next, CheckNull, SetNull)
|
||||
#define DLLRemove(f, l, n) DLLRemove_NPZ(f, l, n, next, prev, CheckNull, SetNull)
|
||||
|
||||
|
||||
// SLL, queue or stack.
|
||||
// These are from rjf's layer.
|
||||
@ -155,78 +144,63 @@
|
||||
// Queue
|
||||
// Queue has only a next pointer. But we can push from front also.
|
||||
// zchk = zero check, zset = zero set
|
||||
#define QueuePush_NZ(f, l, n, next, zchk, zset) \
|
||||
(zchk(f) ? (((f) = (l) = (n)), zset((n)->next)) \
|
||||
: ((l)->next = (n), (l) = (n), zset((n)->next)))
|
||||
#define QueuePush_NZ(f, l, n, next, zchk, zset) (zchk(f)?\
|
||||
(((f)=(l)=(n)), zset((n)->next)):\
|
||||
((l)->next=(n),(l)=(n),zset((n)->next)))
|
||||
|
||||
#define QueuePushFront_NZ(f, l, n, next, zchk, zset) \
|
||||
(zchk(f) ? ((f) = (l) = (n)), zset((n)->next) : ((n)->next = (f)), \
|
||||
((f) = (n)))
|
||||
#define QueuePushFront_NZ(f, l, n, next, zchk, zset) ( zchk(f) ? \
|
||||
((f)=(l)=(n)), zset((n)->next) : \
|
||||
((n)->next = (f)), ((f) = (n)) )
|
||||
|
||||
#define QueuePop_NZ(f, l, next, zchk, zset) \
|
||||
((f) == (l) ? (zset(f), zset(l)) : ((f) = (f)->next))
|
||||
#define QueuePop_NZ(f, l, next, zchk, zset) ( (f)==(l) ? \
|
||||
(zset(f), zset(l)) : ((f)=(f)->next))
|
||||
|
||||
#define QueuePush(f, l, n) QueuePush_NZ(f, l, n, next, CheckNull, SetNull)
|
||||
#define QueuePushFront(f, l, n) \
|
||||
QueuePushFront_NZ(f, l, n, next, CheckNull, SetNull)
|
||||
#define QueuePushFront(f, l, n) QueuePushFront_NZ(f, l, n, next, CheckNull, SetNull)
|
||||
#define QueuePop(f, l) QueuePop_NZ(f, l, next, CheckNull, SetNull)
|
||||
|
||||
////////////////
|
||||
// Stack
|
||||
#define StackPush_N(f, n, next) \
|
||||
((n)->next = (f), (f) = (n)) // Take the first element and set it to n->next,
|
||||
// and set the first element to the node n.
|
||||
#define StackPop_NZ(f, next, zchk) \
|
||||
(zchk(f) ? 0 : ((f) = (f)->next)) // If first element is not zero we say that
|
||||
// the first element is f->next, ie we pop f
|
||||
// and put f->next on top.
|
||||
#define StackPush_N(f, n, next) ((n)->next=(f), (f)=(n)) // Take the first element and set it to n->next, and set the first element to the node n.
|
||||
#define StackPop_NZ(f, next, zchk) (zchk(f) ? 0 : ((f)=(f)->next)) // If first element is not zero we say that the first element is f->next, ie we pop f and put f->next on top.
|
||||
|
||||
#define StackPush(f, n) StackPush_N(f, n, next)
|
||||
#define StackPop(f) StackPop_NZ(f, next, CheckNull)
|
||||
|
||||
/////////////////////////
|
||||
//- Clamp/min/max
|
||||
#define Min(a, b) (((a) < (b)) ? (a) : (b))
|
||||
#define Max(a, b) (((a) > (b)) ? (a) : (b))
|
||||
#define ClampTop(x, a) \
|
||||
Min(x, a) // "Top" since we are cutting off anything above Min(x,a)
|
||||
#define ClampBot(a, x) \
|
||||
Max(a, x) // "Bot" since we're cutting off anything below Max(a,x)
|
||||
#define Min(a, b) (((a)<(b)) ? (a) : (b))
|
||||
#define Max(a, b) (((a)>(b)) ? (a) : (b))
|
||||
#define ClampTop(x, a) Min(x,a) // "Top" since we are cutting off anything above Min(x,a)
|
||||
#define ClampBot(a, x) Max(a,x) // "Bot" since we're cutting off anything below Max(a,x)
|
||||
// If a > x we get a, else we see if b < x and then get b if true, else x.
|
||||
// TODO(anton): Is this actually what we want from a Clamp?
|
||||
#define Clamp(a, x, b) (((a) > (x)) ? (a) : ((b) < (x)) ? (b) : (x))
|
||||
#define Clamp(a, x, b) (((a)>(x))?(a):((b)<(x))?(b):(x))
|
||||
|
||||
//- loop
|
||||
#define DeferLoop(start, end) \
|
||||
for (int _i_ = ((start), 0); _i_ == 0; _i_ += 1, (end))
|
||||
#define DeferLoopChecked(begin, end) \
|
||||
for (int _i_ = 2 * !(begin); (_i_ == 2 ? ((end), 0) : !_i_); _i_ += 1, (end))
|
||||
#define DeferLoop(start, end) for(int _i_ = ((start), 0); _i_ == 0; _i_ += 1, (end))
|
||||
#define DeferLoopChecked(begin, end) for(int _i_ = 2 * !(begin); (_i_ == 2 ? ((end), 0) : !_i_); _i_ += 1, (end))
|
||||
|
||||
#define EachEnumVal(type, it) type it = (type)0; it < type##_COUNT; it = (type)(it+1)
|
||||
#define EachNonZeroEnumVal(type, it) type it = (type)1; it < type##_COUNT; it = (type)(it+1)
|
||||
|
||||
#define EachEnumVal(type, it) \
|
||||
type it = (type)0; \
|
||||
it < type##_COUNT; \
|
||||
it = (type)(it + 1)
|
||||
#define EachNonZeroEnumVal(type, it) \
|
||||
type it = (type)1; \
|
||||
it < type##_COUNT; \
|
||||
it = (type)(it + 1)
|
||||
|
||||
/////////////////////////
|
||||
//~ Base types
|
||||
typedef int8_t S8;
|
||||
typedef int16_t S16;
|
||||
typedef int32_t S32;
|
||||
typedef int64_t S64;
|
||||
typedef uint8_t U8;
|
||||
typedef int8_t S8;
|
||||
typedef int16_t S16;
|
||||
typedef int32_t S32;
|
||||
typedef int64_t S64;
|
||||
typedef uint8_t U8;
|
||||
typedef uint16_t U16;
|
||||
typedef uint32_t U32;
|
||||
typedef uint64_t U64;
|
||||
typedef S8 B8;
|
||||
typedef S16 B16;
|
||||
typedef S32 B32;
|
||||
typedef S64 B64;
|
||||
typedef float F32;
|
||||
typedef double F64;
|
||||
typedef S8 B8;
|
||||
typedef S16 B16;
|
||||
typedef S32 B32;
|
||||
typedef S64 B64;
|
||||
typedef float F32;
|
||||
typedef double F64;
|
||||
typedef void VoidFunction(void);
|
||||
|
||||
/////////////////////////
|
||||
@ -240,6 +214,7 @@ read_only global U64 U64Max = 0xFFFFFFFFFFFFFFFF;
|
||||
// TODO(anton): Rest of the limits, unsigned and signed integer values
|
||||
read_only global U32 SignF32 = 0x80000000;
|
||||
|
||||
|
||||
//- compiler, shut up! helpers
|
||||
#define unused_variable(name) (void)name
|
||||
|
||||
@ -247,32 +222,41 @@ read_only global U32 SignF32 = 0x80000000;
|
||||
//~ Base enums
|
||||
|
||||
// Describing a 2-coordinate system
|
||||
typedef enum Axis2 { Axis2_Invalid = -1, Axis2_X, Axis2_Y, Axis2_COUNT } Axis2;
|
||||
typedef enum Axis2
|
||||
{
|
||||
Axis2_Invalid = -1,
|
||||
Axis2_X,
|
||||
Axis2_Y,
|
||||
Axis2_COUNT
|
||||
}
|
||||
Axis2;
|
||||
#define Axis2_flip(a) ((Axis2)(!(a)))
|
||||
|
||||
// Corners of a rectangle.
|
||||
// 00 ----- 10
|
||||
// | |
|
||||
// 01 ----- 11
|
||||
typedef enum Corner {
|
||||
typedef enum Corner
|
||||
{
|
||||
Corner_Invalid = -1,
|
||||
Corner_00,
|
||||
Corner_01,
|
||||
Corner_10,
|
||||
Corner_11,
|
||||
Corner_COUNT
|
||||
} Corner;
|
||||
}
|
||||
Corner;
|
||||
|
||||
////////////////////////////////
|
||||
//~ Member Offset Helper
|
||||
|
||||
typedef struct MemberOffset MemberOffset;
|
||||
struct MemberOffset {
|
||||
struct MemberOffset
|
||||
{
|
||||
U64 v;
|
||||
};
|
||||
|
||||
#define MemberOff(S, member) \
|
||||
(MemberOffset) { OffsetOf(S, member) }
|
||||
#define MemberOff(S, member) (MemberOffset){OffsetOf(S, member)}
|
||||
#define MemberOffLit(S, member) {OffsetOf(S, member)}
|
||||
#define MemberFromOff(ptr, type, memoff) (*(type *)((U8 *)ptr + memoff.v))
|
||||
|
||||
@ -280,36 +264,26 @@ struct MemberOffset {
|
||||
//~ Assertions
|
||||
|
||||
#if OS_WINDOWS
|
||||
#define break_debugger() __debugbreak()
|
||||
# define break_debugger() __debugbreak()
|
||||
#else
|
||||
#error not implemented
|
||||
# error not implemented
|
||||
#endif
|
||||
|
||||
#undef Assert
|
||||
#define Assert(b) \
|
||||
do { \
|
||||
if (!(b)) { \
|
||||
break_debugger(); \
|
||||
} \
|
||||
} while (0)
|
||||
#define Assert(b) do { if(!(b)) { break_debugger(); } } while(0)
|
||||
|
||||
#if !defined(LOG_NOT_IMPLEMENTED)
|
||||
#define LOG_NOT_IMPLEMENTED \
|
||||
printf("\nFATAL ERROR: Not implemented yet.\n"); \
|
||||
Assert(false); \
|
||||
exit(1);
|
||||
# define LOG_NOT_IMPLEMENTED printf("\nFATAL ERROR: Not implemented yet.\n"); Assert(false); exit(1);
|
||||
#endif
|
||||
|
||||
/////////////////////////
|
||||
//~ Bit patterns
|
||||
#define AlignUpToPow2(bytes_to_align, alignment_bytes) \
|
||||
(((bytes_to_align) + (alignment_bytes - 1)) & ~(alignment_bytes - 1))
|
||||
#define AlignUpToPow2(bytes_to_align, alignment_bytes) (((bytes_to_align) + (alignment_bytes - 1)) & ~(alignment_bytes - 1))
|
||||
|
||||
inline_function F32 absolute_value_F32(F32 f) {
|
||||
union {
|
||||
U32 u;
|
||||
F32 f;
|
||||
} x;
|
||||
inline_function F32
|
||||
absolute_value_F32(F32 f)
|
||||
{
|
||||
union { U32 u; F32 f; } x;
|
||||
x.f = f;
|
||||
x.u = x.u & ~SignF32;
|
||||
return x.f;
|
||||
@ -317,17 +291,4 @@ inline_function F32 absolute_value_F32(F32 f) {
|
||||
|
||||
// TODO(anton): Understand rjf's bit patterns
|
||||
|
||||
///////////////
|
||||
/// LOgging
|
||||
|
||||
#define ENABLE_LOGGING 1
|
||||
#if ENABLE_LOGGING
|
||||
#define LOG(msg) \
|
||||
{ \
|
||||
LPCSTR lpmsg = (LPCSTR)msg; \
|
||||
OutputDebugString(lpmsg); \
|
||||
}
|
||||
#else
|
||||
#define LOG(msg)
|
||||
#endif
|
||||
#endif // BASE_TYPES_H
|
||||
#endif //BASE_TYPES_H
|
||||
|
||||
@ -2,15 +2,15 @@ global BSplineCtx g_bspline_ctx = {0};
|
||||
global Grid g_grid = {0};
|
||||
global U32 g_debug_bspline_matrix = 0;
|
||||
|
||||
|
||||
function F64
|
||||
bspline_recursion(F64 x, U32 k, U32 i)
|
||||
{
|
||||
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.
|
||||
@ -39,6 +39,8 @@ bspline_recursion(F64 x, U32 k, U32 i)
|
||||
|
||||
return term1 + term2;
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
@ -220,16 +222,16 @@ set_up_bsplines_at_points_and_write_matrix_F64(Arena *arena)
|
||||
|
||||
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");
|
||||
}
|
||||
|
||||
{
|
||||
|
||||
@ -1,21 +1,21 @@
|
||||
|
||||
|
||||
function void write_array_binary_F64(String8 path_to_file, F64 *values,
|
||||
U32 array_size) {
|
||||
OS_Handle file_handle =
|
||||
OS_file_open(OS_AccessFlag_Write | OS_AccessFlag_CreateNew, path_to_file);
|
||||
function void
|
||||
write_array_binary_F64(String8 path_to_file, F64 *values, U32 array_size)
|
||||
{
|
||||
OS_Handle file_handle = OS_file_open(OS_AccessFlag_Write | OS_AccessFlag_CreateNew,
|
||||
path_to_file);
|
||||
{
|
||||
ArenaTemp scratch = scratch_get(0, 0);
|
||||
String8List list = {0};
|
||||
String8 temp = {0};
|
||||
temp.str = (U8 *)values;
|
||||
temp.size = sizeof(F64) * array_size;
|
||||
temp.str = (U8*)values;
|
||||
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, str8_lit("Wrote binary array data to"));
|
||||
str8_list_push(scratch.arena, &log_list, path_to_file);
|
||||
StringJoin join = {0};
|
||||
join.sep = str8_lit(" ");
|
||||
@ -27,15 +27,17 @@ function void write_array_binary_F64(String8 path_to_file, F64 *values,
|
||||
OS_file_close(file_handle);
|
||||
}
|
||||
|
||||
function void write_string_list_to_file(Arena *arena, String8 path,
|
||||
String8List *list) {
|
||||
function void
|
||||
write_string_list_to_file(Arena *arena, String8 path, String8List *list)
|
||||
{
|
||||
|
||||
OS_Handle file_handle =
|
||||
OS_file_open(OS_AccessFlag_Write | OS_AccessFlag_CreateNew, path);
|
||||
OS_Handle file_handle = OS_file_open(OS_AccessFlag_Write | OS_AccessFlag_CreateNew,
|
||||
path);
|
||||
OS_file_write(arena, file_handle, 0, *list, 0);
|
||||
|
||||
U32 debug = 1;
|
||||
if (debug) {
|
||||
if(debug)
|
||||
{
|
||||
String8List log_list = {0};
|
||||
str8_list_push(arena, &log_list, str8_lit("Wrote array to"));
|
||||
str8_list_push(arena, &log_list, path);
|
||||
@ -48,20 +50,24 @@ function void write_string_list_to_file(Arena *arena, String8 path,
|
||||
OS_file_close(file_handle);
|
||||
}
|
||||
|
||||
function void write_array_F64(String8 path_to_file, F64 *values, U32 array_size,
|
||||
char *fmt) {
|
||||
|
||||
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++) {
|
||||
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);
|
||||
scratch_release(scratch);
|
||||
}
|
||||
|
||||
function String8 get_eigenvector_filename(Arena *arena, U32 energy_idx,
|
||||
U32 angular_momentum_idx) {
|
||||
String8 out = str8_pushf(arena, "E:\\dev\\hf_again\\out\\eigvec_n%i_l%i.dat",
|
||||
energy_idx, angular_momentum_idx);
|
||||
|
||||
function String8
|
||||
get_eigenvector_filename(Arena *arena, U32 energy_idx, U32 angular_momentum_idx)
|
||||
{
|
||||
String8 out = str8_pushf(arena, "D:\\dev\\hf_again\\out\\eigvec_n%i_l%i.dat", energy_idx, angular_momentum_idx);
|
||||
return out;
|
||||
}
|
||||
|
||||
@ -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);
|
||||
|
||||
@ -1,18 +1,18 @@
|
||||
#ifndef HF_BASE_H
|
||||
#define HF_BASE_H
|
||||
|
||||
#include <mkl_cblas.h>
|
||||
#include <mkl_lapack.h>
|
||||
|
||||
// Complex number with double precision
|
||||
typedef struct Z64 Z64;
|
||||
struct Z64 {
|
||||
struct Z64
|
||||
{
|
||||
F64 re;
|
||||
F64 im;
|
||||
};
|
||||
|
||||
typedef struct Mat_F64 Mat_F64;
|
||||
struct Mat_F64 {
|
||||
struct Mat_F64
|
||||
{
|
||||
U32 size1;
|
||||
U32 size2;
|
||||
F64 **matrix;
|
||||
@ -20,7 +20,8 @@ struct Mat_F64 {
|
||||
};
|
||||
|
||||
typedef struct GaussLegendre GaussLegendre;
|
||||
struct GaussLegendre {
|
||||
struct GaussLegendre
|
||||
{
|
||||
U32 order;
|
||||
F64 *weights;
|
||||
F64 *abscissae;
|
||||
@ -40,10 +41,10 @@ function void print_mat_F64(Mat_F64 *mat);
|
||||
function void mat_invert_F64(Mat_F64 *mat);
|
||||
|
||||
// Gauss-Legendre
|
||||
function void set_up_gauss_legendre_points(Arena *arena);
|
||||
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 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 );
|
||||
|
||||
#endif /* HF_BASE_H */
|
||||
|
||||
179
src/hf/tests.c
179
src/hf/tests.c
@ -2,72 +2,83 @@
|
||||
|
||||
typedef F64 (*func_F64)(F64);
|
||||
|
||||
function F64 cos2(F64 x) { return cos(x) * cos(x); }
|
||||
function F64
|
||||
cos2(F64 x)
|
||||
{
|
||||
return cos(x)*cos(x);
|
||||
}
|
||||
|
||||
function F64 gq_integration(F64 a, F64 b, func_F64 func) {
|
||||
function F64
|
||||
gq_integration(F64 a, F64 b, func_F64 func)
|
||||
{
|
||||
|
||||
F64 prefac = 0.5 * (b - a);
|
||||
F64 prefac = 0.5*(b-a);
|
||||
F64 gq_sum = 0.0;
|
||||
|
||||
for (U32 i = 0; i < g_gauss_legendre.order; i++) {
|
||||
for(U32 i = 0; i < g_gauss_legendre.order; i++)
|
||||
{
|
||||
F64 w = g_gauss_legendre.weights[i];
|
||||
F64 z = g_gauss_legendre.abscissae[i];
|
||||
F64 shift = (z * prefac) + ((a + b) * 0.5);
|
||||
F64 shift = (z*prefac)+((a+b)*0.5);
|
||||
F64 funct_value_at_shift = func(shift);
|
||||
gq_sum = gq_sum + prefac * w * funct_value_at_shift;
|
||||
gq_sum = gq_sum + prefac*w*funct_value_at_shift;
|
||||
}
|
||||
|
||||
return gq_sum;
|
||||
|
||||
}
|
||||
|
||||
function void test_gauss_legendre() {
|
||||
function void
|
||||
test_gauss_legendre()
|
||||
{
|
||||
|
||||
// Test GL qudrature by integrating cos^2 from 0 to 2pi, it should equal pi.
|
||||
{
|
||||
ArenaTemp scratch = scratch_get(0, 0);
|
||||
ArenaTemp scratch = scratch_get(0,0);
|
||||
|
||||
F64 pi = 4.0 * atan(1.0);
|
||||
String8 pi_out = str8_pushf(
|
||||
scratch.arena, "Pi from standard library 4.0*atan(1.0) = %.16f \n", pi);
|
||||
F64 pi = 4.0*atan(1.0);
|
||||
String8 pi_out = str8_pushf(scratch.arena, "Pi from standard library 4.0*atan(1.0) = %.16f \n", pi);
|
||||
LOG(pi_out.str);
|
||||
|
||||
F64 a = 0.0;
|
||||
F64 b = 2.0 * pi;
|
||||
F64 b = 2.0*pi;
|
||||
F64 gq_sum_single_interval = gq_integration(a, b, cos2);
|
||||
|
||||
String8 out = str8_pushf(
|
||||
scratch.arena,
|
||||
"Integration result for a single interval: %.16f, Reference: %.16f \n",
|
||||
gq_sum_single_interval, pi);
|
||||
String8 out = str8_pushf(scratch.arena, "Integration result for a single interval: %.16f, Reference: %.16f \n", gq_sum_single_interval, pi);
|
||||
LOG(out.str);
|
||||
|
||||
// Test with several smaller intervals instead.
|
||||
F64 aa[10] = {0.0};
|
||||
F64 bb[10] = {0.0};
|
||||
F64 delta = b / 10.0;
|
||||
for (U32 i = 0; i < 10; i++) {
|
||||
aa[i] = i * delta;
|
||||
bb[i] = (i + 1) * delta;
|
||||
String8 intervals =
|
||||
str8_pushf(scratch.arena, "%i, a=%f, b=%f \n", i, aa[i], bb[i]);
|
||||
F64 delta = b/10.0;
|
||||
for(U32 i = 0; i < 10; i++)
|
||||
{
|
||||
aa[i] = i*delta;
|
||||
bb[i] = (i+1)*delta;
|
||||
String8 intervals = str8_pushf(scratch.arena, "%i, a=%f, b=%f \n", i, aa[i], bb[i]);
|
||||
LOG(intervals.str);
|
||||
}
|
||||
|
||||
F64 gq_sum_several_intervals = 0.0;
|
||||
for (U32 i = 0; i < 10; i++) {
|
||||
for(U32 i = 0; i < 10; i++)
|
||||
{
|
||||
gq_sum_several_intervals += gq_integration(aa[i], bb[i], cos2);
|
||||
|
||||
}
|
||||
|
||||
out = str8_pushf(
|
||||
scratch.arena,
|
||||
"Integration result for ten intervals: %.16f, Reference: %.16f \n",
|
||||
gq_sum_several_intervals, pi);
|
||||
out = str8_pushf(scratch.arena, "Integration result for ten intervals: %.16f, Reference: %.16f \n", gq_sum_several_intervals, pi);
|
||||
LOG(out.str);
|
||||
|
||||
|
||||
scratch_release(scratch);
|
||||
}
|
||||
}
|
||||
|
||||
function void mkl_things(void) {
|
||||
|
||||
|
||||
|
||||
function void mkl_things(void)
|
||||
{
|
||||
|
||||
OS_InitReceipt os_receipt = OS_init();
|
||||
OS_InitGfxReceipt os_gfx_receipt = OS_gfx_init(os_receipt);
|
||||
@ -75,90 +86,102 @@ function void mkl_things(void) {
|
||||
Arena *arena = m_make_arena();
|
||||
|
||||
U32 N = 4;
|
||||
Z64 *main_A = PushArray(arena, Z64, N * N);
|
||||
main_A[0] = (Z64){-3.84, 2.25};
|
||||
main_A[1] = (Z64){-0.66, 0.83};
|
||||
main_A[2] = (Z64){-3.99, -4.73};
|
||||
main_A[3] = (Z64){7.74, 4.18};
|
||||
main_A[4] = (Z64){-8.94, -4.75};
|
||||
main_A[5] = (Z64){-4.40, -3.82};
|
||||
main_A[6] = (Z64){-5.88, -6.60};
|
||||
main_A[7] = (Z64){3.66, -7.53};
|
||||
main_A[8] = (Z64){8.95, -6.53};
|
||||
main_A[9] = (Z64){-3.50, -4.26};
|
||||
Z64 *main_A = PushArray(arena, Z64, N*N);
|
||||
main_A[0] = (Z64){-3.84, 2.25};
|
||||
main_A[1] = (Z64){-0.66, 0.83};
|
||||
main_A[2] = (Z64){-3.99, -4.73};
|
||||
main_A[3] = (Z64){ 7.74, 4.18};
|
||||
main_A[4] = (Z64){-8.94, -4.75};
|
||||
main_A[5] = (Z64){-4.40, -3.82};
|
||||
main_A[6] = (Z64){-5.88, -6.60};
|
||||
main_A[7] = (Z64){ 3.66, -7.53};
|
||||
main_A[8] = (Z64){ 8.95, -6.53};
|
||||
main_A[9] = (Z64){-3.50, -4.26};
|
||||
main_A[10] = (Z64){-3.36, -0.40};
|
||||
main_A[11] = (Z64){2.58, 3.60};
|
||||
main_A[12] = (Z64){-9.87, 4.82};
|
||||
main_A[13] = (Z64){-3.15, 7.36};
|
||||
main_A[14] = (Z64){-0.75, 5.23};
|
||||
main_A[15] = (Z64){4.59, 5.41};
|
||||
main_A[11] = (Z64){ 2.58, 3.60};
|
||||
main_A[12] = (Z64){-9.87, 4.82};
|
||||
main_A[13] = (Z64){-3.15, 7.36};
|
||||
main_A[14] = (Z64){-0.75, 5.23};
|
||||
main_A[15] = (Z64){ 4.59, 5.41};
|
||||
|
||||
LOG("\n\n---- Calling Intel MKL zgeev test (Using Z64 instead of "
|
||||
"MKL_Complex16 etc) ---- \n\n");
|
||||
LOG("\n\n---- Calling Intel MKL zgeev test (Using Z64 instead of MKL_Complex16 etc) ---- \n\n");
|
||||
|
||||
{
|
||||
S32 n = N, lda = N, ldvl = N, ldvr = N, info, lwork;
|
||||
Z64 wkopt;
|
||||
Z64 *work;
|
||||
|
||||
F64 *rwork = PushArray(arena, F64, 2 * N);
|
||||
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);
|
||||
for (U32 j = 0; j < N; j++) {
|
||||
for (U32 i = 0; i < N; i++) {
|
||||
U32 index = i * N + j;
|
||||
Z64 *vl = PushArray(arena, Z64, N*N);
|
||||
Z64 *vr = PushArray(arena, Z64, N*N);
|
||||
Z64 *a = PushArray(arena, Z64, N*N);
|
||||
for(U32 j = 0; j < N; j++)
|
||||
{
|
||||
for(U32 i = 0; i < N; i++)
|
||||
{
|
||||
U32 index = i*N+j;
|
||||
|
||||
a[index] = main_A[index];
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
/* Executable statements */
|
||||
LOG(" ZGEEV Example Program Results\n");
|
||||
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);
|
||||
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));
|
||||
work = (Z64*)malloc( lwork*sizeof(Z64) );
|
||||
/* Solve eigenproblem */
|
||||
zgeev("Vectors", "Vectors", &n, a, &lda, w, vl, &ldvl, vr, &ldvr, work,
|
||||
&lwork, rwork, &info);
|
||||
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);
|
||||
if( info > 0 ) {
|
||||
LOG( "The algorithm failed to compute eigenvalues.\n" );
|
||||
exit( 1 );
|
||||
}
|
||||
/* Print eigenvalues */
|
||||
print_matrix_Z64("Eigenvalues", 1, n, w, 1);
|
||||
print_matrix_Z64( "Eigenvalues", 1, n, w, 1 );
|
||||
/* Print left eigenvectors */
|
||||
print_matrix_Z64("Left eigenvectors", n, n, vl, ldvl);
|
||||
print_matrix_Z64( "Left eigenvectors", n, n, vl, ldvl );
|
||||
/* Print right eigenvectors */
|
||||
print_matrix_Z64("Right eigenvectors", n, n, vr, ldvr);
|
||||
print_matrix_Z64( "Right eigenvectors", n, n, vr, ldvr );
|
||||
/* Free workspace */
|
||||
free((void *)work);
|
||||
free( (void*)work );
|
||||
} /* End of ZGEEV Example */
|
||||
|
||||
|
||||
|
||||
LOG("\n\n--- End of EntryPoint, exiting program. \n\n");
|
||||
}
|
||||
|
||||
function void test_matrix() {
|
||||
|
||||
function void
|
||||
test_matrix()
|
||||
{
|
||||
|
||||
Mat_F64 test_mat = mat_F64(5, 5);
|
||||
|
||||
{
|
||||
ArenaTemp scratch = scratch_get(0, 0);
|
||||
ArenaTemp scratch = scratch_get(0,0);
|
||||
|
||||
for (U32 i = 0; i < test_mat.size1; i++) {
|
||||
for (U32 j = 0; j < test_mat.size2; j++) {
|
||||
for(U32 i = 0; i < test_mat.size1; i++)
|
||||
{
|
||||
for(U32 j = 0; j < test_mat.size2; j++)
|
||||
{
|
||||
F64 val = i * test_mat.size2 + j;
|
||||
mat_F64_set(&test_mat, i, j, val);
|
||||
}
|
||||
}
|
||||
|
||||
for (U32 i = 0; i < test_mat.size1; i++) {
|
||||
for (U32 j = 0; j < test_mat.size2; j++) {
|
||||
for(U32 i = 0; i < test_mat.size1; i++)
|
||||
{
|
||||
for(U32 j = 0; j < test_mat.size2; j++)
|
||||
{
|
||||
F64 val = mat_F64_get(&test_mat, i, j);
|
||||
String8 out_str = str8_pushf(scratch.arena, " %2.2f", val);
|
||||
LOG(out_str.str);
|
||||
@ -170,10 +193,14 @@ function void test_matrix() {
|
||||
}
|
||||
}
|
||||
|
||||
function void testing_MKL() {
|
||||
|
||||
function void
|
||||
testing_MKL()
|
||||
{
|
||||
|
||||
test_mkl_zgeev();
|
||||
test_mkl_dsyevd();
|
||||
|
||||
}
|
||||
|
||||
/* function void */
|
||||
@ -196,14 +223,12 @@ function void testing_MKL() {
|
||||
/* printf( " DGEEV Example Program Results\n" ); */
|
||||
/* /1* Query and allocate the optimal workspace *1/ */
|
||||
/* lwork = -1; */
|
||||
/* dgeev( "Vectors", "Vectors", &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr,
|
||||
*/
|
||||
/* dgeev( "Vectors", "Vectors", &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr, */
|
||||
/* &wkopt, &lwork, &info ); */
|
||||
/* lwork = (MKL_INT)wkopt; */
|
||||
/* work = (double*)malloc( lwork*sizeof(double) ); */
|
||||
/* /1* Solve eigenproblem *1/ */
|
||||
/* dgeev( "Vectors", "Vectors", &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr,
|
||||
*/
|
||||
/* dgeev( "Vectors", "Vectors", &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr, */
|
||||
/* work, &lwork, &info ); */
|
||||
/* /1* Check for convergence *1/ */
|
||||
/* if( info > 0 ) { */
|
||||
|
||||
332
src/main.c
332
src/main.c
@ -1,27 +1,40 @@
|
||||
|
||||
#define ENABLE_LOGGING 1
|
||||
#if ENABLE_LOGGING
|
||||
#define LOG(msg) { OutputDebugString(msg); }
|
||||
#else
|
||||
#define LOG(msg)
|
||||
#endif
|
||||
|
||||
// ---
|
||||
// Header includes
|
||||
#include <stdlib.h>
|
||||
#include <math.h>
|
||||
#include <mkl_cblas.h>
|
||||
#include "base/base_inc.h"
|
||||
#include "os/os_inc.h"
|
||||
|
||||
#include "hf/bsplines_and_grid.h"
|
||||
#include "hf/file_io.h"
|
||||
#include "hf/hf_base.h"
|
||||
#include "hf/bsplines_and_grid.h"
|
||||
|
||||
// ---
|
||||
// .C includes
|
||||
#include "base/base_inc.c"
|
||||
#include "os/os_entry_point.c"
|
||||
#include "os/os_inc.c"
|
||||
#include "os/os_entry_point.c"
|
||||
|
||||
#include "hf/bsplines_and_grid.c"
|
||||
#include "hf/file_io.c"
|
||||
#include "hf/hf_base.c"
|
||||
#include "hf/bsplines_and_grid.c"
|
||||
|
||||
#include "hf/tests.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 {
|
||||
struct Eigensolution_F64
|
||||
{
|
||||
F64 *eigenvalues_re;
|
||||
F64 *eigenvalues_im;
|
||||
Mat_F64 right_eigenvectors;
|
||||
@ -29,38 +42,47 @@ struct Eigensolution_F64 {
|
||||
};
|
||||
|
||||
typedef struct Orbital Orbital;
|
||||
struct Orbital {
|
||||
struct Orbital
|
||||
{
|
||||
U32 n;
|
||||
U32 l;
|
||||
U32 j;
|
||||
Eigensolution_F64 eigensolution;
|
||||
};
|
||||
|
||||
|
||||
typedef struct Atom Atom;
|
||||
struct Atom {
|
||||
struct Atom
|
||||
{
|
||||
U32 N;
|
||||
Orbital *orbitals;
|
||||
};
|
||||
|
||||
typedef struct SortPair_F64 SortPair_F64;
|
||||
struct SortPair_F64 {
|
||||
struct SortPair_F64
|
||||
{
|
||||
F64 value;
|
||||
U64 original_index;
|
||||
};
|
||||
|
||||
global Arena *g_base_arena = 0;
|
||||
global Arena *g_filename_arena = 0;
|
||||
|
||||
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) {
|
||||
|
||||
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) {
|
||||
for (U64 j = low; j < high; j++)
|
||||
{
|
||||
if(array[j].value <= pivot)
|
||||
{
|
||||
i += 1;
|
||||
temp = array[i];
|
||||
array[i] = array[j];
|
||||
@ -76,198 +98,221 @@ U64 qsort_partition_F64(SortPair_F64 *array, U64 size, U64 low, U64 high) {
|
||||
return final_pivot_pos;
|
||||
}
|
||||
|
||||
function void qsort_F64(SortPair_F64 *array, U64 size, U64 low, U64 high) {
|
||||
if (low < high) {
|
||||
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) {
|
||||
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++) {
|
||||
for(U64 i = 0; i < size; i++)
|
||||
{
|
||||
pairs[i].value = array[i];
|
||||
pairs[i].original_index = i;
|
||||
}
|
||||
|
||||
qsort_F64(pairs, size, 0, size - 1);
|
||||
qsort_F64(pairs, size, 0, size-1);
|
||||
|
||||
for (U32 i = 0; i < size; i++) {
|
||||
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) {
|
||||
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++) {
|
||||
for(U64 i = 0; i < size; i++)
|
||||
{
|
||||
U64 original_index = indices[i];
|
||||
temp[i] = array[original_index];
|
||||
}
|
||||
for (U64 i = 0; i < size; i++) {
|
||||
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);
|
||||
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);
|
||||
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]);
|
||||
//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 );
|
||||
//printf(" (%6.2f,%6.2f)", a[i+j*lda].real, a[i+j*lda].imag );
|
||||
}
|
||||
LOG(newline.str);
|
||||
// printf("\n");
|
||||
//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
|
||||
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); */
|
||||
/* } */
|
||||
|
||||
/* { */
|
||||
/* 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;
|
||||
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))) {
|
||||
F64 abs_index_diff = abs((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;
|
||||
// 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;
|
||||
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++) {
|
||||
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);
|
||||
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++) {
|
||||
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 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 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);
|
||||
*/
|
||||
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);
|
||||
}
|
||||
// mat_F64_set(&H, i, j, abs_index_diff);
|
||||
//mat_F64_set(&H, i, j, abs_index_diff);
|
||||
}
|
||||
// LOG("\n");
|
||||
//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);
|
||||
//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);
|
||||
//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);
|
||||
|
||||
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++) {
|
||||
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);
|
||||
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++) {
|
||||
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 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);
|
||||
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;
|
||||
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);
|
||||
String8 out = str8_pushf(scratch.arena, "n:%i, l:%i norm: %.2f \n", n, l, norm);
|
||||
LOG(out.str);
|
||||
scratch_release(scratch);
|
||||
}
|
||||
|
||||
function void EntryPoint(void) {
|
||||
|
||||
function void
|
||||
EntryPoint(void)
|
||||
{
|
||||
|
||||
OS_InitReceipt os_receipt = OS_init();
|
||||
OS_InitGfxReceipt os_gfx_receipt = OS_gfx_init(os_receipt);
|
||||
@ -281,18 +326,16 @@ function void EntryPoint(void) {
|
||||
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");
|
||||
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");
|
||||
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.
|
||||
//- 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;
|
||||
@ -310,16 +353,16 @@ function void EntryPoint(void) {
|
||||
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
|
||||
// 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++) {
|
||||
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);
|
||||
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)
|
||||
@ -331,12 +374,13 @@ function void EntryPoint(void) {
|
||||
// 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);
|
||||
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};
|
||||
@ -352,22 +396,23 @@ function void EntryPoint(void) {
|
||||
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);
|
||||
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);
|
||||
&wkopt, &lwork, &info);
|
||||
lwork = (S32)wkopt;
|
||||
work = (F64 *)malloc(lwork * sizeof(F64));
|
||||
work = (F64 *)malloc(lwork * sizeof(F64) );
|
||||
dgeev("Vectors", "Vectors", &size1, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr,
|
||||
work, &lwork, &info);
|
||||
if (info > 0) {
|
||||
work, &lwork, &info);
|
||||
if(info > 0)
|
||||
{
|
||||
LOG("Failed to compute eigenvalues in dgeev\n");
|
||||
exit(1);
|
||||
exit( 1 );
|
||||
}
|
||||
|
||||
eigensolution.right_eigenvectors = mat_F64_from_data(ldvr, size1, vr);
|
||||
@ -375,35 +420,40 @@ function void EntryPoint(void) {
|
||||
U64 *sorted_indices = PushArray(mkl_arena, U64, size1);
|
||||
sort_and_get_indices_F64(wr, sorted_indices, size1);
|
||||
sort_by_indices_F64(wi, sorted_indices, size1);
|
||||
// print_eigenvalues( "Eigenvalues sorted: ", size1, wr, wi );
|
||||
//print_eigenvalues( "Eigenvalues sorted: ", size1, wr, wi );
|
||||
|
||||
|
||||
U32 i = 0;
|
||||
F64 energy = -1000.0;
|
||||
U32 counter = 0;
|
||||
while (energy < 0.0) {
|
||||
while(energy < 0.0)
|
||||
{
|
||||
energy = wr[i];
|
||||
U64 energy_index = sorted_indices[i];
|
||||
U64 n = 1 + i;
|
||||
if (ang_mom_idx > 0) {
|
||||
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");
|
||||
//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) {
|
||||
if(counter > 10)
|
||||
{
|
||||
break;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
free((void *)work);
|
||||
free( (void*)work );
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user