anton/hf_again
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Compare commits

base: anton:d670d018b714dd5989b260c2f62bfe2659aadef3
Branches Tags
anton:main
anton:clang_test
..
compare: anton:fc2c897727799dad746db03c9a97c00823fecba5
Branches Tags
anton:main
anton:clang_test

No commits in common. "d670d018b714dd5989b260c2f62bfe2659aadef3" and "fc2c897727799dad746db03c9a97c00823fecba5" have entirely different histories.
d670d018b7 ... fc2c897727

15 changed files with 721 additions and 568 deletions
Show Stats Expand all files Collapse all files

40
build.bat
View File

@ -1,17 +1,29 @@
@echo off @echo off
set CommonCompilerFlags=/nologo /Zi /FC /Od
set mkl_root=E:/lib/intel_mkl/mkl/2025.3 ctime -begin timeBuild.ctm
set mkl_core=%mkl_root%/lib/mkl_core.lib
set mkl_intel_lp64=%mkl_root%/lib/mkl_intel_lp64.lib @rem /WX /W4 /wd4201 /wd4100 /wd4189 /wd4244 /wd4127 /wd4456
set mkl_intel_thread=%mkl_root%/lib/mkl_intel_thread.lib @rem set CommonCompilerFlags="/nologo /Zi /FC"
set MKLCOMPILER=E:/lib/intel_mkl/compiler/2025.3 set CommonCompilerFlags=/nologo /Zi /FC /Od
set libiomp5md=%MKLCOMPILER%/lib/libiomp5md.lib @rem /WX /W4 /wd4201 /wd4100 /wd4189 /wd4244 /wd4127 /wd4456
set libiompdll_path=E:\lib\intel_mkl\compiler\2025.3\bin @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_name=libiomp5md.dll
set libiompdll=%libiompdll_path%\%libiompdll_name% set libiompdll=%libiompdll_path%\%libiompdll_name%
set Sources=../src/main.c set Sources=../src/main.c
IF NOT EXIST .\out mkdir .\out IF NOT EXIST .\out mkdir .\out
IF NOT EXIST .\build mkdir .\build IF NOT EXIST .\build mkdir .\build
pushd .\build pushd .\build
@ -24,19 +36,15 @@ if not exist "%libiompdll_name%" (
) else ( ) else (
echo Copied openmp dll: %libiompdll_name% echo Copied openmp dll: %libiompdll_name%
) )
) )
cl %CommonCompilerFlags% %Sources% /I"%mkl_root%\include" ^ cl %CommonCompilerFlags% %Sources% /I"%mkl_root%\include" /link %mkl_core% %mkl_intel_lp64% %mkl_intel_thread% %libiomp5md%
/link %mkl_core% %mkl_intel_lp64% %mkl_intel_thread% %libiomp5md%
set LastError=%ERRORLEVEL% set LastError=%ERRORLEVEL%
popd popd
ctime -end timeBuild.ctm %LastError%
echo Build complete
set PATH=%PATH%;%mkl_root%\bin
set PATH=%PATH%;%MKLCOMPILER%\bin
set LastError=%ERRORLEVEL%
IF NOT %LastError%==0 GOTO :end IF NOT %LastError%==0 GOTO :end
:end :end

20
notebook/hf_again.ipynb
View File

File diff suppressed because one or more lines are too long

68
out/eigvec_n1_l0.dat
View File

@ -1,34 +1,34 @@
4.079148e-01 -4.079148e-01
6.284958e-01 -6.284958e-01
5.853385e-01 -5.853385e-01
2.830382e-01 -2.830382e-01
1.166302e-01 -1.166302e-01
4.437820e-02 -4.437820e-02
1.608212e-02 -1.608212e-02
5.645430e-03 -5.645430e-03
1.935000e-03 -1.935000e-03
6.523220e-04 -6.523220e-04
2.166921e-04 -2.166921e-04
7.130839e-05 -7.130839e-05
2.321027e-05 -2.321027e-05
7.519194e-06 -7.519194e-06
2.411514e-06 -2.411514e-06
7.732122e-07 -7.732122e-07
2.449929e-07 -2.449929e-07
7.808282e-08 -7.808282e-08
2.445164e-08 -2.445164e-08
7.783642e-09 -7.783642e-09
2.403077e-09 -2.403077e-09
7.695295e-10 -7.695295e-10
2.325414e-10 -2.325412e-10
7.588069e-11 -7.588092e-11
2.207375e-11 -2.207390e-11
7.527071e-12 -7.526818e-12
2.032901e-12 -2.032911e-12
7.619167e-13 -7.621029e-13
1.762993e-13 -1.765319e-13
8.066878e-14 -8.079409e-14
1.303573e-14 -1.311978e-14
9.445442e-15 -9.257067e-15
-2.760762e-16 1.842392e-16
1.847288e-15 -1.701083e-15

68
out/eigvec_n2_l0.dat
View File

@ -1,34 +1,34 @@
1.067985e-01 -1.067985e-01
9.343962e-02 -9.343962e-02
-1.122171e-01 1.122171e-01
-3.647821e-01 3.647821e-01
-4.769945e-01 4.769945e-01
-4.732031e-01 4.732031e-01
-4.077962e-01 4.077962e-01
-3.223604e-01 3.223604e-01
-2.403128e-01 2.403128e-01
-1.716945e-01 1.716945e-01
-1.187758e-01 1.187758e-01
-8.011225e-02 8.011225e-02
-5.294166e-02 5.294166e-02
-3.440233e-02 3.440233e-02
-2.204206e-02 2.204206e-02
-1.395428e-02 1.395428e-02
-8.743348e-03 8.743348e-03
-5.429335e-03 5.429335e-03
-3.344953e-03 3.344953e-03
-2.046453e-03 2.046453e-03
-1.244260e-03 1.244260e-03
-7.523093e-04 7.523093e-04
-4.525787e-04 4.525787e-04
-2.710226e-04 2.710226e-04
-1.616214e-04 1.616214e-04
-9.600751e-05 9.600751e-05
-5.681871e-05 5.681871e-05
-3.349503e-05 3.349503e-05
-1.964742e-05 1.964742e-05
-1.142544e-05 1.142544e-05
-6.511164e-06 6.511164e-06
-3.501003e-06 3.501003e-06
-1.380734e-06 1.380734e-06
-4.297736e-07 4.297736e-07

68
out/eigvec_n2_l2.dat
View File

@ -1,34 +1,34 @@
-2.548436e-03 2.548436e-03
-1.768047e-02 1.768047e-02
-6.801789e-02 6.801789e-02
-1.434030e-01 1.434030e-01
-2.235356e-01 2.235356e-01
-2.907379e-01 2.907379e-01
-3.357788e-01 3.357788e-01
-3.567533e-01 3.567533e-01
-3.563788e-01 3.563788e-01
-3.395486e-01 3.395486e-01
-3.116223e-01 3.116223e-01
-2.774422e-01 2.774422e-01
-2.408901e-01 2.408901e-01
-2.047894e-01 2.047894e-01
-1.709988e-01 1.709988e-01
-1.405887e-01 1.405887e-01
-1.140369e-01 1.140369e-01
-9.140819e-02 9.140819e-02
-7.250238e-02 7.250238e-02
-5.696782e-02 5.696782e-02
-4.438296e-02 4.438296e-02
-3.431126e-02 3.431126e-02
-2.633551e-02 2.633551e-02
-2.007697e-02 2.007697e-02
-1.520409e-02 1.520409e-02
-1.143438e-02 1.143438e-02
-8.531931e-03 8.531931e-03
-6.302777e-03 6.302777e-03
-4.589046e-03 4.589046e-03
-3.262862e-03 3.262862e-03
-2.220364e-03 2.220364e-03
-1.376048e-03 1.376048e-03
-5.929916e-04 5.929916e-04
-1.933374e-04 1.933374e-04

34
out/eigvec_n6_l2.dat Normal file
View 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
View 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
View 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

187
src/base/base_core.h
View File

@ -1,8 +1,9 @@
#ifndef BASE_TYPES_H #ifndef BASE_TYPES_H
#define BASE_TYPES_H #define BASE_TYPES_H
#include <math.h>
#include <stdint.h> #include <stdint.h>
#include <stdlib.h>
#include <stdio.h>
///////////////////////// /////////////////////////
//~ Macros //~ Macros
@ -21,8 +22,8 @@
///////////////////////// /////////////////////////
//- Keywords //- Keywords
// Static is stupid and means different things depending on context in C and // Static is stupid and means different things depending on context in C and C++.
// C++. These defines increases readability. // These defines increases readability.
#define function static // Function internal to compilation unit. #define function static // Function internal to compilation unit.
#define local_persist static #define local_persist static
#define global static #define global static
@ -54,19 +55,13 @@
///////////////////////// /////////////////////////
//- Memory operations //- Memory operations
// It's nice to put these in macros, so we can swap out the functionality from // It's nice to put these in macros, so we can swap out the functionality from standard library, eventually.
// standard library, eventually.
#define MemoryCopy memcpy #define MemoryCopy memcpy
#define MemoryMove memmove #define MemoryMove memmove
#define MemorySet memset #define MemorySet memset
// NOTE(anton): This gives a 4127 compiler warning for the sizeof conditional. // NOTE(anton): This gives a 4127 compiler warning for the sizeof conditional. This should be ignored
// This should be ignored #define MemoryCopyStruct(dst, src) do { Assert(sizeof(*(dst)) == sizeof(*(src))); MemoryCopy((dst), (src), sizeof(*(dst))); } while(0)
#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 MemoryZero(ptr, size) MemorySet((ptr), 0, (size))
#define MemoryZeroStruct(ptr) MemoryZero((ptr), sizeof(*(ptr))) #define MemoryZeroStruct(ptr) MemoryZero((ptr), sizeof(*(ptr)))
@ -82,8 +77,7 @@
// TODO(anton): Understand why this becomes offset actually // TODO(anton): Understand why this becomes offset actually
#define OffsetOf(type, member_name) IntFromPtr(&Member(type, member_name)) #define OffsetOf(type, member_name) IntFromPtr(&Member(type, member_name))
// TODO(anton): Understand this // TODO(anton): Understand this
#define BaseFromMember(type, member_name, ptr) \ #define BaseFromMember(type, member_name, ptr) (type *)((U8 *)(ptr) - OffsetOf(type, member_name))
(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 Kilobytes(n) (n << 10) // 2^10 == 1024 etc
@ -104,49 +98,44 @@
#define SetNull(p) ((p)=0) #define SetNull(p) ((p)=0)
// Link list helper macros that are a bit involved // 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 // Suffixes N,P,Z means that we have (N)ext, (P)rev arguments and/or a (Z)ero check and/or set argument
// check and/or set argument f, l, n are "first", "last", "node" I think? DLL // f, l, n are "first", "last", "node" I think?
// Doubly Linked List: Each node has a prev and next pointer. Operations: Push // DLL
// back, Push front, remove // 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) \ #define DLLInsert_NPZ(f,l,p,n,next,prev,zchk,zset) \
(zchk(f) ? (((f) = (l) = (n)), zset((n)->next), zset((n)->prev)) \ (zchk(f) ? (((f) = (l) = (n)), zset((n)->next), zset((n)->prev)) :\
: zchk(p) ? (zset((n)->prev), (n)->next = (f), \ zchk(p) ? (zset((n)->prev), (n)->next = (f), (zchk(f) ? (0) : ((f)->prev = (n))), (f) = (n)) :\
(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),\
: ((zchk((p)->next) ? (0) : (((p)->next->prev) = (n))), \
(n)->next = (p)->next, (n)->prev = (p), (p)->next = (n), \
((p) == (l) ? (l) = (n) : (0)))) ((p) == (l) ? (l) = (n) : (0))))
#define DLLPushBack_NPZ(f, 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)
DLLInsert_NPZ(f, l, l, n, next, prev, zchk, zset)
#define DLLPushBack_NP(f, l, n, next, prev, zchk) \ #define DLLPushBack_NP(f, l, n, next, prev, zchk) \
(zchk(f) ? ((f) = (l) = (n), (n)->next = (n)->prev = 0) \ (zchk(f) ? ((f)=(l)=(n),(n)->next=(n)->prev=0) : ((n)->prev=(l),(l)->next=(n),(l)=(n),(n)->next=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. // 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. // 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 // If l != n and f != n we set n->next->prev to n->prev, and n->prev->next to n->next
// n->next
#define DLLRemove_NP(f, l, n, next, prev) \ #define DLLRemove_NP(f, l, n, next, prev) (((f) == (n) ? \
(((f) == (n) ? ((f) = (f)->next, (f)->prev = 0) \ ((f)=(f)->next, (f)->prev=0) : \
: (l) == (n) \ (l) == (n) ? \
? ((l) = (l)->prev, (l)->next = 0) \ ((l)=(l)->prev, (l)->next=0) : \
: ((n)->next->prev = (n)->prev, (n)->prev->next = (n)->next))) ((n)->next->prev=(n)->prev, \
(n)->prev->next=(n)->next) ))
#define DLLRemove_NPZ(f, l, n, next, prev, zchk, zset) \ #define DLLRemove_NPZ(f,l,n,next,prev,zchk,zset) (((f)==(n))?\
(((f) == (n)) ? ((f) = (f)->next, (zchk(f) ? (zset(l)) : zset((f)->prev))) \ ((f)=(f)->next, (zchk(f) ? (zset(l)) : zset((f)->prev))):\
: ((l) == (n)) ? ((l) = (l)->prev, (zchk(l) ? (zset(f)) : zset((l)->next))) \ ((l)==(n))?\
: ((zchk((n)->next) ? (0) : ((n)->next->prev = (n)->prev)), \ ((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)))) (zchk((n)->prev) ? (0) : ((n)->prev->next=(n)->next))))
#define DLLPushBack(f, l, n) \ #define DLLPushBack(f, l, n) DLLPushBack_NPZ(f, l, n, next, prev, CheckNull, SetNull)
DLLPushBack_NPZ(f, l, n, next, prev, CheckNull, SetNull)
// For front push I can just switch prev/next! // For front push I can just switch prev/next!
#define DLLPushFront(f, l, n) \ #define DLLPushFront(f, l, n) DLLPushBack_NPZ(l, f, n, prev, next, CheckNull, SetNull)
DLLPushBack_NPZ(l, f, n, prev, next, CheckNull, SetNull) #define DLLRemove(f, l, n) DLLRemove_NPZ(f, l, n, next, prev, CheckNull, SetNull)
#define DLLRemove(f, l, n) \
DLLRemove_NPZ(f, l, n, next, prev, CheckNull, SetNull)
// SLL, queue or stack. // SLL, queue or stack.
// These are from rjf's layer. // These are from rjf's layer.
@ -155,31 +144,25 @@
// Queue // Queue
// Queue has only a next pointer. But we can push from front also. // Queue has only a next pointer. But we can push from front also.
// zchk = zero check, zset = zero set // zchk = zero check, zset = zero set
#define QueuePush_NZ(f, l, n, next, zchk, zset) \ #define QueuePush_NZ(f, l, n, next, zchk, zset) (zchk(f)?\
(zchk(f) ? (((f) = (l) = (n)), zset((n)->next)) \ (((f)=(l)=(n)), zset((n)->next)):\
: ((l)->next = (n), (l) = (n), zset((n)->next))) ((l)->next=(n),(l)=(n),zset((n)->next)))
#define QueuePushFront_NZ(f, l, n, next, zchk, zset) \ #define QueuePushFront_NZ(f, l, n, next, zchk, zset) ( zchk(f) ? \
(zchk(f) ? ((f) = (l) = (n)), zset((n)->next) : ((n)->next = (f)), \ ((f)=(l)=(n)), zset((n)->next) : \
((f) = (n))) ((n)->next = (f)), ((f) = (n)) )
#define QueuePop_NZ(f, l, next, zchk, zset) \ #define QueuePop_NZ(f, l, next, zchk, zset) ( (f)==(l) ? \
((f) == (l) ? (zset(f), zset(l)) : ((f) = (f)->next)) (zset(f), zset(l)) : ((f)=(f)->next))
#define QueuePush(f, l, n) QueuePush_NZ(f, l, n, next, CheckNull, SetNull) #define QueuePush(f, l, n) QueuePush_NZ(f, l, n, next, CheckNull, SetNull)
#define QueuePushFront(f, l, n) \ #define QueuePushFront(f, l, n) QueuePushFront_NZ(f, l, n, next, CheckNull, SetNull)
QueuePushFront_NZ(f, l, n, next, CheckNull, SetNull)
#define QueuePop(f, l) QueuePop_NZ(f, l, next, CheckNull, SetNull) #define QueuePop(f, l) QueuePop_NZ(f, l, next, CheckNull, SetNull)
//////////////// ////////////////
// Stack // Stack
#define StackPush_N(f, n, next) \ #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.
((n)->next = (f), (f) = (n)) // Take the first element and set it to n->next, #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.
// 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 StackPush(f, n) StackPush_N(f, n, next)
#define StackPop(f) StackPop_NZ(f, next, CheckNull) #define StackPop(f) StackPop_NZ(f, next, CheckNull)
@ -188,28 +171,19 @@
//- Clamp/min/max //- Clamp/min/max
#define Min(a, b) (((a)<(b)) ? (a) : (b)) #define Min(a, b) (((a)<(b)) ? (a) : (b))
#define Max(a, b) (((a)>(b)) ? (a) : (b)) #define Max(a, b) (((a)>(b)) ? (a) : (b))
#define ClampTop(x, a) \ #define ClampTop(x, a) Min(x,a) // "Top" since we are cutting off anything above Min(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 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. // 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? // 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 //- loop
#define DeferLoop(start, end) \ #define DeferLoop(start, end) for(int _i_ = ((start), 0); _i_ == 0; _i_ += 1, (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 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 //~ Base types
@ -240,6 +214,7 @@ read_only global U64 U64Max = 0xFFFFFFFFFFFFFFFF;
// TODO(anton): Rest of the limits, unsigned and signed integer values // TODO(anton): Rest of the limits, unsigned and signed integer values
read_only global U32 SignF32 = 0x80000000; read_only global U32 SignF32 = 0x80000000;
//- compiler, shut up! helpers //- compiler, shut up! helpers
#define unused_variable(name) (void)name #define unused_variable(name) (void)name
@ -247,32 +222,41 @@ read_only global U32 SignF32 = 0x80000000;
//~ Base enums //~ Base enums
// Describing a 2-coordinate system // 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))) #define Axis2_flip(a) ((Axis2)(!(a)))
// Corners of a rectangle. // Corners of a rectangle.
// 00 ----- 10 // 00 ----- 10
// | | // | |
// 01 ----- 11 // 01 ----- 11
typedef enum Corner { typedef enum Corner
{
Corner_Invalid = -1, Corner_Invalid = -1,
Corner_00, Corner_00,
Corner_01, Corner_01,
Corner_10, Corner_10,
Corner_11, Corner_11,
Corner_COUNT Corner_COUNT
} Corner; }
Corner;
//////////////////////////////// ////////////////////////////////
//~ Member Offset Helper //~ Member Offset Helper
typedef struct MemberOffset MemberOffset; typedef struct MemberOffset MemberOffset;
struct MemberOffset { struct MemberOffset
{
U64 v; U64 v;
}; };
#define MemberOff(S, member) \ #define MemberOff(S, member) (MemberOffset){OffsetOf(S, member)}
(MemberOffset) { OffsetOf(S, member) }
#define MemberOffLit(S, member) {OffsetOf(S, member)} #define MemberOffLit(S, member) {OffsetOf(S, member)}
#define MemberFromOff(ptr, type, memoff) (*(type *)((U8 *)ptr + memoff.v)) #define MemberFromOff(ptr, type, memoff) (*(type *)((U8 *)ptr + memoff.v))
@ -286,30 +270,20 @@ struct MemberOffset {
#endif #endif
#undef Assert #undef Assert
#define Assert(b) \ #define Assert(b) do { if(!(b)) { break_debugger(); } } while(0)
do { \
if (!(b)) { \
break_debugger(); \
} \
} while (0)
#if !defined(LOG_NOT_IMPLEMENTED) #if !defined(LOG_NOT_IMPLEMENTED)
#define LOG_NOT_IMPLEMENTED \ # define LOG_NOT_IMPLEMENTED printf("\nFATAL ERROR: Not implemented yet.\n"); Assert(false); exit(1);
printf("\nFATAL ERROR: Not implemented yet.\n"); \
Assert(false); \
exit(1);
#endif #endif
///////////////////////// /////////////////////////
//~ Bit patterns //~ Bit patterns
#define AlignUpToPow2(bytes_to_align, alignment_bytes) \ #define AlignUpToPow2(bytes_to_align, alignment_bytes) (((bytes_to_align) + (alignment_bytes - 1)) & ~(alignment_bytes - 1))
(((bytes_to_align) + (alignment_bytes - 1)) & ~(alignment_bytes - 1))
inline_function F32 absolute_value_F32(F32 f) { inline_function F32
union { absolute_value_F32(F32 f)
U32 u; {
F32 f; union { U32 u; F32 f; } x;
} x;
x.f = f; x.f = f;
x.u = x.u & ~SignF32; x.u = x.u & ~SignF32;
return x.f; return x.f;
@ -317,17 +291,4 @@ inline_function F32 absolute_value_F32(F32 f) {
// TODO(anton): Understand rjf's bit patterns // 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

26
src/hf/bsplines_and_grid.c
View File

@ -2,15 +2,15 @@ global BSplineCtx g_bspline_ctx = {0};
global Grid g_grid = {0}; global Grid g_grid = {0};
global U32 g_debug_bspline_matrix = 0; global U32 g_debug_bspline_matrix = 0;
function F64 function F64
bspline_recursion(F64 x, U32 k, U32 i) 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(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): // TODO(anton):
// This is like a hack to get the last bspline to be 1 at the last point. // 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; 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); 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("D:\\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("D:\\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("D:\\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("D:\\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("D:\\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("D:\\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("D:\\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("D:\\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("D:\\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\\dBspline9.dat"), dBspl9, num_grid_points, "%13.6e\n");
} }
{ {

42
src/hf/file_io.c
View File

@ -1,9 +1,10 @@
function void write_array_binary_F64(String8 path_to_file, F64 *values, function void
U32 array_size) { 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); OS_Handle file_handle = OS_file_open(OS_AccessFlag_Write | OS_AccessFlag_CreateNew,
path_to_file);
{ {
ArenaTemp scratch = scratch_get(0, 0); ArenaTemp scratch = scratch_get(0, 0);
String8List list = {0}; String8List list = {0};
@ -14,8 +15,7 @@ function void write_array_binary_F64(String8 path_to_file, F64 *values,
OS_file_write(scratch.arena, file_handle, 0, list, 0); OS_file_write(scratch.arena, file_handle, 0, list, 0);
String8List log_list = {0}; String8List log_list = {0};
str8_list_push(scratch.arena, &log_list, str8_list_push(scratch.arena, &log_list, str8_lit("Wrote binary array data to"));
str8_lit("Wrote binary array data to"));
str8_list_push(scratch.arena, &log_list, path_to_file); str8_list_push(scratch.arena, &log_list, path_to_file);
StringJoin join = {0}; StringJoin join = {0};
join.sep = str8_lit(" "); 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); OS_file_close(file_handle);
} }
function void write_string_list_to_file(Arena *arena, String8 path, function void
String8List *list) { write_string_list_to_file(Arena *arena, String8 path, String8List *list)
{
OS_Handle file_handle = OS_Handle file_handle = OS_file_open(OS_AccessFlag_Write | OS_AccessFlag_CreateNew,
OS_file_open(OS_AccessFlag_Write | OS_AccessFlag_CreateNew, path); path);
OS_file_write(arena, file_handle, 0, *list, 0); OS_file_write(arena, file_handle, 0, *list, 0);
U32 debug = 1; U32 debug = 1;
if (debug) { if(debug)
{
String8List log_list = {0}; String8List log_list = {0};
str8_list_push(arena, &log_list, str8_lit("Wrote array to")); str8_list_push(arena, &log_list, str8_lit("Wrote array to"));
str8_list_push(arena, &log_list, path); 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); 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); 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);
scratch_release(scratch); scratch_release(scratch);
} }
function String8 get_eigenvector_filename(Arena *arena, U32 energy_idx,
U32 angular_momentum_idx) { function String8
String8 out = str8_pushf(arena, "E:\\dev\\hf_again\\out\\eigvec_n%i_l%i.dat", get_eigenvector_filename(Arena *arena, U32 energy_idx, U32 angular_momentum_idx)
energy_idx, 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; return out;
} }

14
src/hf/file_io.h
View File

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

11
src/hf/hf_base.h
View File

@ -1,18 +1,18 @@
#ifndef HF_BASE_H #ifndef HF_BASE_H
#define HF_BASE_H #define HF_BASE_H
#include <mkl_cblas.h>
#include <mkl_lapack.h>
// Complex number with double precision // Complex number with double precision
typedef struct Z64 Z64; typedef struct Z64 Z64;
struct Z64 { struct Z64
{
F64 re; F64 re;
F64 im; F64 im;
}; };
typedef struct Mat_F64 Mat_F64; typedef struct Mat_F64 Mat_F64;
struct Mat_F64 { struct Mat_F64
{
U32 size1; U32 size1;
U32 size2; U32 size2;
F64 **matrix; F64 **matrix;
@ -20,7 +20,8 @@ struct Mat_F64 {
}; };
typedef struct GaussLegendre GaussLegendre; typedef struct GaussLegendre GaussLegendre;
struct GaussLegendre { struct GaussLegendre
{
U32 order; U32 order;
F64 *weights; F64 *weights;
F64 *abscissae; F64 *abscissae;

99
src/hf/tests.c
View File

@ -2,14 +2,21 @@
typedef F64 (*func_F64)(F64); 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; 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 w = g_gauss_legendre.weights[i];
F64 z = g_gauss_legendre.abscissae[i]; F64 z = g_gauss_legendre.abscissae[i];
F64 shift = (z*prefac)+((a+b)*0.5); F64 shift = (z*prefac)+((a+b)*0.5);
@ -18,56 +25,60 @@ function F64 gq_integration(F64 a, F64 b, func_F64 func) {
} }
return gq_sum; 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. // 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); F64 pi = 4.0*atan(1.0);
String8 pi_out = str8_pushf( String8 pi_out = str8_pushf(scratch.arena, "Pi from standard library 4.0*atan(1.0) = %.16f \n", pi);
scratch.arena, "Pi from standard library 4.0*atan(1.0) = %.16f \n", pi);
LOG(pi_out.str); LOG(pi_out.str);
F64 a = 0.0; F64 a = 0.0;
F64 b = 2.0*pi; F64 b = 2.0*pi;
F64 gq_sum_single_interval = gq_integration(a, b, cos2); F64 gq_sum_single_interval = gq_integration(a, b, cos2);
String8 out = str8_pushf( String8 out = str8_pushf(scratch.arena, "Integration result for a single interval: %.16f, Reference: %.16f \n", gq_sum_single_interval, pi);
scratch.arena,
"Integration result for a single interval: %.16f, Reference: %.16f \n",
gq_sum_single_interval, pi);
LOG(out.str); LOG(out.str);
// Test with several smaller intervals instead. // Test with several smaller intervals instead.
F64 aa[10] = {0.0}; F64 aa[10] = {0.0};
F64 bb[10] = {0.0}; F64 bb[10] = {0.0};
F64 delta = b/10.0; F64 delta = b/10.0;
for (U32 i = 0; i < 10; i++) { for(U32 i = 0; i < 10; i++)
{
aa[i] = i*delta; aa[i] = i*delta;
bb[i] = (i+1)*delta; bb[i] = (i+1)*delta;
String8 intervals = String8 intervals = str8_pushf(scratch.arena, "%i, a=%f, b=%f \n", i, aa[i], bb[i]);
str8_pushf(scratch.arena, "%i, a=%f, b=%f \n", i, aa[i], bb[i]);
LOG(intervals.str); LOG(intervals.str);
} }
F64 gq_sum_several_intervals = 0.0; 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); gq_sum_several_intervals += gq_integration(aa[i], bb[i], cos2);
} }
out = str8_pushf( out = str8_pushf(scratch.arena, "Integration result for ten intervals: %.16f, Reference: %.16f \n", gq_sum_several_intervals, pi);
scratch.arena,
"Integration result for ten intervals: %.16f, Reference: %.16f \n",
gq_sum_several_intervals, pi);
LOG(out.str); LOG(out.str);
scratch_release(scratch); scratch_release(scratch);
} }
} }
function void mkl_things(void) {
function void mkl_things(void)
{
OS_InitReceipt os_receipt = OS_init(); OS_InitReceipt os_receipt = OS_init();
OS_InitGfxReceipt os_gfx_receipt = OS_gfx_init(os_receipt); OS_InitGfxReceipt os_gfx_receipt = OS_gfx_init(os_receipt);
@ -93,8 +104,7 @@ function void mkl_things(void) {
main_A[14] = (Z64){-0.75, 5.23}; main_A[14] = (Z64){-0.75, 5.23};
main_A[15] = (Z64){ 4.59, 5.41}; main_A[15] = (Z64){ 4.59, 5.41};
LOG("\n\n---- Calling Intel MKL zgeev test (Using Z64 instead of " LOG("\n\n---- Calling Intel MKL zgeev test (Using Z64 instead of MKL_Complex16 etc) ---- \n\n");
"MKL_Complex16 etc) ---- \n\n");
{ {
S32 n = N, lda = N, ldvl = N, ldvr = N, info, lwork; S32 n = N, lda = N, ldvl = N, ldvr = N, info, lwork;
@ -106,25 +116,29 @@ function void mkl_things(void) {
Z64 *vl = PushArray(arena, Z64, N*N); Z64 *vl = PushArray(arena, Z64, N*N);
Z64 *vr = PushArray(arena, Z64, N*N); Z64 *vr = PushArray(arena, Z64, N*N);
Z64 *a = PushArray(arena, Z64, N*N); Z64 *a = PushArray(arena, Z64, N*N);
for (U32 j = 0; j < N; j++) { for(U32 j = 0; j < N; j++)
for (U32 i = 0; i < N; i++) { {
for(U32 i = 0; i < N; i++)
{
U32 index = i*N+j; U32 index = i*N+j;
a[index] = main_A[index]; a[index] = main_A[index];
} }
} }
/* Executable statements */ /* Executable statements */
LOG( " ZGEEV Example Program Results\n" ); LOG( " ZGEEV Example Program Results\n" );
/* Query and allocate the optimal workspace */ /* Query and allocate the optimal workspace */
lwork = -1; lwork = -1;
zgeev("Vectors", "Vectors", &n, a, &lda, w, vl, &ldvl, vr, &ldvr, &wkopt, zgeev( "Vectors", "Vectors", &n, a, &lda, w, vl, &ldvl, vr, &ldvr,
&lwork, rwork, &info); &wkopt, &lwork, rwork, &info );
lwork = (S32)wkopt.re; lwork = (S32)wkopt.re;
work = (Z64*)malloc( lwork*sizeof(Z64) ); work = (Z64*)malloc( lwork*sizeof(Z64) );
/* Solve eigenproblem */ /* Solve eigenproblem */
zgeev("Vectors", "Vectors", &n, a, &lda, w, vl, &ldvl, vr, &ldvr, work, zgeev( "Vectors", "Vectors", &n, a, &lda, w, vl, &ldvl, vr, &ldvr,
&lwork, rwork, &info); work, &lwork, rwork, &info );
/* Check for convergence */ /* Check for convergence */
if( info > 0 ) { if( info > 0 ) {
LOG( "The algorithm failed to compute eigenvalues.\n" ); LOG( "The algorithm failed to compute eigenvalues.\n" );
@ -140,25 +154,34 @@ function void mkl_things(void) {
free( (void*)work ); free( (void*)work );
} /* End of ZGEEV Example */ } /* End of ZGEEV Example */
LOG("\n\n--- End of EntryPoint, exiting program. \n\n"); 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); 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 i = 0; i < test_mat.size1; i++)
for (U32 j = 0; j < test_mat.size2; j++) { {
for(U32 j = 0; j < test_mat.size2; j++)
{
F64 val = i * test_mat.size2 + j; F64 val = i * test_mat.size2 + j;
mat_F64_set(&test_mat, i, j, val); mat_F64_set(&test_mat, i, j, val);
} }
} }
for (U32 i = 0; i < test_mat.size1; i++) { for(U32 i = 0; i < test_mat.size1; i++)
for (U32 j = 0; j < test_mat.size2; j++) { {
for(U32 j = 0; j < test_mat.size2; j++)
{
F64 val = mat_F64_get(&test_mat, i, j); F64 val = mat_F64_get(&test_mat, i, j);
String8 out_str = str8_pushf(scratch.arena, " %2.2f", val); String8 out_str = str8_pushf(scratch.arena, " %2.2f", val);
LOG(out_str.str); 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_zgeev();
test_mkl_dsyevd(); test_mkl_dsyevd();
} }
/* function void */ /* function void */
@ -196,14 +223,12 @@ function void testing_MKL() {
/* printf( " DGEEV Example Program Results\n" ); */ /* printf( " DGEEV Example Program Results\n" ); */
/* /1* Query and allocate the optimal workspace *1/ */ /* /1* Query and allocate the optimal workspace *1/ */
/* lwork = -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 ); */ /* &wkopt, &lwork, &info ); */
/* lwork = (MKL_INT)wkopt; */ /* lwork = (MKL_INT)wkopt; */
/* work = (double*)malloc( lwork*sizeof(double) ); */ /* work = (double*)malloc( lwork*sizeof(double) ); */
/* /1* Solve eigenproblem *1/ */ /* /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 ); */ /* work, &lwork, &info ); */
/* /1* Check for convergence *1/ */ /* /1* Check for convergence *1/ */
/* if( info > 0 ) { */ /* if( info > 0 ) { */

230
src/main.c
View File

@ -1,27 +1,40 @@
#define ENABLE_LOGGING 1
#if ENABLE_LOGGING
#define LOG(msg) { OutputDebugString(msg); }
#else
#define LOG(msg)
#endif
// --- // ---
// Header includes // Header includes
#include <stdlib.h>
#include <math.h>
#include <mkl_cblas.h>
#include "base/base_inc.h" #include "base/base_inc.h"
#include "os/os_inc.h" #include "os/os_inc.h"
#include "hf/bsplines_and_grid.h"
#include "hf/file_io.h" #include "hf/file_io.h"
#include "hf/hf_base.h" #include "hf/hf_base.h"
#include "hf/bsplines_and_grid.h"
// --- // ---
// .C includes // .C includes
#include "base/base_inc.c" #include "base/base_inc.c"
#include "os/os_entry_point.c"
#include "os/os_inc.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/file_io.c"
#include "hf/hf_base.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; typedef struct Eigensolution_F64 Eigensolution_F64;
struct Eigensolution_F64 { struct Eigensolution_F64
{
F64 *eigenvalues_re; F64 *eigenvalues_re;
F64 *eigenvalues_im; F64 *eigenvalues_im;
Mat_F64 right_eigenvectors; Mat_F64 right_eigenvectors;
@ -29,38 +42,47 @@ struct Eigensolution_F64 {
}; };
typedef struct Orbital Orbital; typedef struct Orbital Orbital;
struct Orbital { struct Orbital
{
U32 n; U32 n;
U32 l; U32 l;
U32 j; U32 j;
Eigensolution_F64 eigensolution; Eigensolution_F64 eigensolution;
}; };
typedef struct Atom Atom; typedef struct Atom Atom;
struct Atom { struct Atom
{
U32 N; U32 N;
Orbital *orbitals; Orbital *orbitals;
}; };
typedef struct SortPair_F64 SortPair_F64; typedef struct SortPair_F64 SortPair_F64;
struct SortPair_F64 { struct SortPair_F64
{
F64 value; F64 value;
U64 original_index; U64 original_index;
}; };
global Arena* g_base_arena = 0; global Arena* g_base_arena = 0;
global Arena* g_filename_arena = 0; global Arena* g_filename_arena = 0;
#define NUM_ANGULAR_MOMENTA 3 #define NUM_ANGULAR_MOMENTA 3
global F64 angular_momenta[NUM_ANGULAR_MOMENTA] = {0.0, 1.0, 2.0}; global F64 angular_momenta[NUM_ANGULAR_MOMENTA] = {0.0, 1.0, 2.0};
global F64 *temp_wavefunction_F64; 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; F64 pivot = array[high].value;
U64 i = low - 1; U64 i = low - 1;
SortPair_F64 temp = {0}; SortPair_F64 temp = {0};
for (U64 j = low; j < high; j++) { for (U64 j = low; j < high; j++)
if (array[j].value <= pivot) { {
if(array[j].value <= pivot)
{
i += 1; i += 1;
temp = array[i]; temp = array[i];
array[i] = array[j]; array[i] = array[j];
@ -76,46 +98,61 @@ U64 qsort_partition_F64(SortPair_F64 *array, U64 size, U64 low, U64 high) {
return final_pivot_pos; return final_pivot_pos;
} }
function void qsort_F64(SortPair_F64 *array, U64 size, U64 low, U64 high) { function void qsort_F64(SortPair_F64 *array, U64 size, U64 low, U64 high)
if (low < high) { {
if(low < high)
{
U64 pivot_index = qsort_partition_F64(array, size, low, high); U64 pivot_index = qsort_partition_F64(array, size, low, high);
qsort_F64(array, size, low, pivot_index - 1); qsort_F64(array, size, low, pivot_index - 1);
qsort_F64(array, size, pivot_index + 1, high); 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)); 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].value = array[i];
pairs[i].original_index = 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; array[i] = pairs[i].value;
indices[i] = pairs[i].original_index; indices[i] = pairs[i].original_index;
} }
free(pairs); 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)); 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]; U64 original_index = indices[i];
temp[i] = array[original_index]; temp[i] = array[original_index];
} }
for (U64 i = 0; i < size; i++) { for(U64 i = 0; i < size; i++)
{
array[i] = temp[i]; array[i] = temp[i];
} }
free(temp); free(temp);
} }
/* Auxiliary routine: printing a matrix */ /* Auxiliary routine: printing a matrix */
function void print_eigenvalues(char *desc, int n, F64 *wr, F64 *wi) { function void
print_eigenvalues( char* desc, int n, F64* wr, F64* wi)
{
ArenaTemp scratch = scratch_get(0,0); ArenaTemp scratch = scratch_get(0,0);
int i, j; int i, j;
String8 newline = str8_lit("\n"); String8 newline = str8_lit("\n");
@ -123,8 +160,7 @@ function void print_eigenvalues(char *desc, int n, F64 *wr, F64 *wi) {
LOG(header.str); LOG(header.str);
//printf("\n %s \n", desc); //printf("\n %s \n", desc);
for( j = 0; j < n; j++ ) { for( j = 0; j < n; j++ ) {
String8 outstr = String8 outstr = str8_pushf(scratch.arena, " (%4.5f, %4.5f)\n", wr[j], wi[j]);
str8_pushf(scratch.arena, " (%4.5f, %4.5f)\n", wr[j], wi[j]);
LOG(outstr.str); 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 );
} }
@ -134,17 +170,19 @@ function void print_eigenvalues(char *desc, int n, F64 *wr, F64 *wi) {
scratch_release(scratch); scratch_release(scratch);
} }
function void set_up_first_matrices(Mat_F64 *H, Mat_F64 *H_l, Mat_F64 *B_inv) { 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 = // We work in units hbar = 1, bohr radius a0 = 1, electron mass m_e = 1, and charge e = 1,
// -0.5*d^2/dr^2 + l(l+1)/(2r^2) - Z/r // 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); ArenaTemp scratch = scratch_get(0, 0);
/* { */ /* { */
/* LOG("Setting up matrices.\n"); */ /* LOG("Setting up matrices.\n"); */
/* String8 setuplog = str8_pushf(scratch.arena, */ /* String8 setuplog = str8_pushf(scratch.arena, */
/* "Num knotpoints: %i, Bspline order k = %i, Matrix size1 = N-k-2 /* "Num knotpoints: %i, Bspline order k = %i, Matrix size1 = N-k-2 = %i, size2 = %i\n", */
* = %i, size2 = %i\n", */
/* N, k, mat_size1, mat_size2); */ /* N, k, mat_size1, mat_size2); */
/* LOG(setuplog.str); */ /* LOG(setuplog.str); */
/* } */ /* } */
@ -154,25 +192,25 @@ function void set_up_first_matrices(Mat_F64 *H, Mat_F64 *H_l, Mat_F64 *B_inv) {
U32 k = g_bspline_ctx.order; U32 k = g_bspline_ctx.order;
// Skipping first bspline // Skipping first bspline
for (U32 i = 0; i < H->size1; i++) { for(U32 i = 0; i < H->size1; i++)
for (U32 j = 0; j < H->size2; j++) { {
U32 bspl_index_i = for(U32 j = 0; j < H->size2; j++)
i + 1; // The second Bspline has index 1 in our array etc. {
U32 bspl_index_i = i+1; // The second Bspline has index 1 in our array etc.
U32 bspl_index_j = j+1; U32 bspl_index_j = j+1;
// This logic assumes 1-indexed bsplines // This logic assumes 1-indexed bsplines
F64 abs_index_diff = F64 abs_index_diff = abs((F64)(bspl_index_i+1) - (F64)(bspl_index_j+1));
fabs((F64)(bspl_index_i + 1) - (F64)(bspl_index_j + 1)); if(!(abs_index_diff > ((F64)k-1.0)))
if (!(abs_index_diff > ((F64)k - 1.0))) { {
// We do Gaussian quadrature between each knot point, // We do Gaussian quadrature between each knot point,
// so we need to figure out where to start. // so we need to figure out where to start.
// We start integration in the first shared knotpoint, which is the // We start integration in the first shared knotpoint, which is the one of the highest index.
// one of the highest index. U32 start_knotpoint_index = bspl_index_i < bspl_index_j ?
U32 start_knotpoint_index = bspl_index_j : bspl_index_i;
bspl_index_i < bspl_index_j ? bspl_index_j : bspl_index_i;
// And we integrate over the next k knotpoints. // And we integrate over the next k knotpoints.
U32 end_knotpoint_index = U32 end_knotpoint_index = bspl_index_i < bspl_index_j ?
bspl_index_i < bspl_index_j ? bspl_index_i + k : bspl_index_j + k; bspl_index_i+k : bspl_index_j+k;
F64 term1 = 0.0; F64 term1 = 0.0;
F64 term2 = 0.0; F64 term2 = 0.0;
@ -180,14 +218,18 @@ function void set_up_first_matrices(Mat_F64 *H, Mat_F64 *H_l, Mat_F64 *B_inv) {
F64 Bmat_term = 0.0; F64 Bmat_term = 0.0;
for(U32 knotpoint_idx = start_knotpoint_index; for(U32 knotpoint_idx = start_knotpoint_index;
knotpoint_idx < end_knotpoint_index; knotpoint_idx++) { knotpoint_idx < end_knotpoint_index;
knotpoint_idx++)
{
F64 a = t[knotpoint_idx]; F64 a = t[knotpoint_idx];
F64 b = t[knotpoint_idx+1]; F64 b = t[knotpoint_idx+1];
F64 prefac = 0.5*(b-a); F64 prefac = 0.5*(b-a);
// Only integrate non-zero intervals // Only integrate non-zero intervals
if (prefac > 1e-16) { if(prefac > 1e-16)
for (U32 gq_i = 0; gq_i < g_gauss_legendre.order; gq_i++) { {
for(U32 gq_i = 0; gq_i < g_gauss_legendre.order; gq_i++)
{
F64 w = g_gauss_legendre.weights[gq_i]; F64 w = g_gauss_legendre.weights[gq_i];
F64 z = g_gauss_legendre.abscissae[gq_i]; F64 z = g_gauss_legendre.abscissae[gq_i];
F64 r = (z*prefac)+((a+b)*0.5); F64 r = (z*prefac)+((a+b)*0.5);
@ -206,12 +248,8 @@ function void set_up_first_matrices(Mat_F64 *H, Mat_F64 *H_l, Mat_F64 *B_inv) {
F64 H_term_sum = 0.5*term1 + (-Z)*term3; F64 H_term_sum = 0.5*term1 + (-Z)*term3;
F64 H_l_term = 0.5*term2; F64 H_l_term = 0.5*term2;
/* String8 debug = str8_pushf(scratch.arena, /* 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", */
* "(i=%i,j=%i,t_i=%4.4f,t_i=%4.4f,term1=%.4e,term2=%.4e,term3=%.4e,term_sum=%.4e) /* bspl_index_i, bspl_index_j, t[bspl_index_i+k-1],t[bspl_index_j+k-1],term1,term2,term3,term_sum); */
* \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); */ /* LOG(debug.str); */
mat_F64_set(H, i, j, H_term_sum); mat_F64_set(H, i, j, H_term_sum);
mat_F64_set(H_l, i, j, H_l_term); mat_F64_set(H_l, i, j, H_l_term);
@ -223,37 +261,42 @@ function void set_up_first_matrices(Mat_F64 *H, Mat_F64 *H_l, Mat_F64 *B_inv) {
//LOG("\n"); //LOG("\n");
} }
LOG(str8_pushf(scratch.arena, "H.size1=N-k-2=%i, last bspline index=%i \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);
H->size1, g_bspline_ctx.num_bsplines - 1)
.str);
scratch_release(scratch); scratch_release(scratch);
print_mat_F64(H); //print_mat_F64(&H);
LOG("\n"); //LOG("\n");
// print_mat_F64(B); //print_mat_F64(&B);
} }
} }
function void compute_wf_norm_F64(F64 *coeffs, U64 coeff_size, U64 n, U64 l) {
function void
compute_wf_norm_F64(F64 *coeffs, U64 coeff_size, U64 n, U64 l)
{
ArenaTemp scratch = scratch_get(0,0); ArenaTemp scratch = scratch_get(0,0);
// Gauss legendre integration // Gauss legendre integration
// //
F64 norm = 0.0; 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 a = g_grid.points[i];
F64 b = g_grid.points[i+1]; F64 b = g_grid.points[i+1];
F64 prefac = 0.5*(b-a); F64 prefac = 0.5*(b-a);
// Only integrate non-zero intervals // Only integrate non-zero intervals
if (prefac > 1e-16) { if(prefac > 1e-16)
for (U32 gq_i = 0; gq_i < g_gauss_legendre.order; gq_i++) { {
for(U32 gq_i = 0; gq_i < g_gauss_legendre.order; gq_i++)
{
F64 w = g_gauss_legendre.weights[gq_i]; F64 w = g_gauss_legendre.weights[gq_i];
F64 z = g_gauss_legendre.abscissae[gq_i]; F64 z = g_gauss_legendre.abscissae[gq_i];
F64 r = (z*prefac)+((a+b)*0.5); F64 r = (z*prefac)+((a+b)*0.5);
F64 term_prefac = (prefac*w); F64 term_prefac = (prefac*w);
F64 wf_at_r = 0.0; F64 wf_at_r = 0.0;
for (U64 j = 0; j < coeff_size; j++) { for(U64 j = 0; j < coeff_size; j++)
{
wf_at_r += coeffs[j]*compute_bspline_F64(r, j + 1); wf_at_r += coeffs[j]*compute_bspline_F64(r, j + 1);
} }
@ -261,13 +304,15 @@ function void compute_wf_norm_F64(F64 *coeffs, U64 coeff_size, U64 n, U64 l) {
} }
} }
} }
String8 out = String8 out = str8_pushf(scratch.arena, "n:%i, l:%i norm: %.2f \n", n, l, norm);
str8_pushf(scratch.arena, "n:%i, l:%i norm: %.2f \n", n, l, norm);
LOG(out.str); LOG(out.str);
scratch_release(scratch); scratch_release(scratch);
} }
function void EntryPoint(void) {
function void
EntryPoint(void)
{
OS_InitReceipt os_receipt = OS_init(); OS_InitReceipt os_receipt = OS_init();
OS_InitGfxReceipt os_gfx_receipt = OS_gfx_init(os_receipt); OS_InitGfxReceipt os_gfx_receipt = OS_gfx_init(os_receipt);
@ -281,18 +326,16 @@ function void EntryPoint(void) {
set_up_grid(g_base_arena); set_up_grid(g_base_arena);
temp_wavefunction_F64 = (F64 *)PushArray(g_base_arena, F64, g_grid.num_steps); 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, write_array_binary_F64(str8_lit(grid_file_path_bin), g_grid.points, g_grid.num_steps);
g_grid.num_steps); write_array_F64(str8_lit(grid_file_path), g_grid.points, g_grid.num_steps, "%13.6e\n");
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. //- The BSpline context is the knotpoints and the BSpline order etc.
set_up_bspline_context(g_base_arena); set_up_bspline_context(g_base_arena);
write_array_F64(str8_lit(knotpoints_file_path), g_bspline_ctx.knotpoints, write_array_F64(str8_lit(knotpoints_file_path),
g_bspline_ctx.num_knotpoints, "%13.6e\n"); g_bspline_ctx.knotpoints, g_bspline_ctx.num_knotpoints,
"%13.6e\n");
//- Then we generate the BSplines and save them off for reference and //- Then we generate the BSplines and save them off for reference and debugging.
// debugging.
set_up_bsplines_at_points_and_write_matrix_F64(g_base_arena); set_up_bsplines_at_points_and_write_matrix_F64(g_base_arena);
U32 N = g_bspline_ctx.num_knotpoints; U32 N = g_bspline_ctx.num_knotpoints;
@ -310,12 +353,12 @@ function void EntryPoint(void) {
set_up_first_matrices(&H_base, &H_l_base, &B_inv); set_up_first_matrices(&H_base, &H_l_base, &B_inv);
// Our problem is Hc = EBc, but we want to solve B^-1Hc = Ec, // 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 // so we invert the B matrix and compute the product A = B^-1H before calling zgeev
// zgeev
mat_invert_F64(&B_inv); mat_invert_F64(&B_inv);
// For each angular momentum // 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); mat_F64_copy_to_dst(&H, &H_base);
F64 l = angular_momenta[ang_mom_idx]; F64 l = angular_momenta[ang_mom_idx];
if(l > 1e-16) { if(l > 1e-16) {
@ -331,12 +374,13 @@ function void EntryPoint(void) {
// Multiply to get A = B^-1 H // Multiply to get A = B^-1 H
{ {
S32 n = A.size1; S32 n = A.size1;
cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n, n, n, 1.0, cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
B_inv.data, n, H.data, n, 0.0, A.data, n); n, n, n, 1.0, B_inv.data, n, H.data, n, 0.0, A.data, n);
LOG("Matrix A: \n"); //LOG("Matrix A: \n");
print_mat_F64(&A); //print_mat_F64(&A);
} }
// This arena is used to push results from f. ex eigenvalue computations. // This arena is used to push results from f. ex eigenvalue computations.
Arena *mkl_arena = m_make_arena(); Arena *mkl_arena = m_make_arena();
Eigensolution_F64 eigensolution = {0}; Eigensolution_F64 eigensolution = {0};
@ -365,7 +409,8 @@ function void EntryPoint(void) {
work = (F64 *)malloc(lwork * sizeof(F64) ); work = (F64 *)malloc(lwork * sizeof(F64) );
dgeev("Vectors", "Vectors", &size1, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr, dgeev("Vectors", "Vectors", &size1, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr,
work, &lwork, &info); work, &lwork, &info);
if (info > 0) { if(info > 0)
{
LOG("Failed to compute eigenvalues in dgeev\n"); LOG("Failed to compute eigenvalues in dgeev\n");
exit( 1 ); exit( 1 );
} }
@ -377,33 +422,38 @@ function void EntryPoint(void) {
sort_by_indices_F64(wi, 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; U32 i = 0;
F64 energy = -1000.0; F64 energy = -1000.0;
U32 counter = 0; U32 counter = 0;
while (energy < 0.0) { while(energy < 0.0)
{
energy = wr[i]; energy = wr[i];
U64 energy_index = sorted_indices[i]; U64 energy_index = sorted_indices[i];
U64 n = 1 + i; U64 n = 1 + i;
if (ang_mom_idx > 0) { if(ang_mom_idx > 0)
{
n = 2 + i; n = 2 + i;
} }
// compute_wf_norm_F64(eigensolution.right_eigenvectors.matrix[energy_index],
// size1, n, ang_mom_idx);
write_array_F64( //compute_wf_norm_F64(eigensolution.right_eigenvectors.matrix[energy_index], size1, n, ang_mom_idx);
get_eigenvector_filename(g_filename_arena, n, ang_mom_idx),
eigensolution.right_eigenvectors.matrix[energy_index], size1, write_array_F64(get_eigenvector_filename(g_filename_arena, n, ang_mom_idx),
"%13.6e\n"); eigensolution.right_eigenvectors.matrix[energy_index], size1, "%13.6e\n");
i += 1; i += 1;
counter += 1; counter += 1;
if (counter > 10) { if(counter > 10)
{
break; break;
} }
} }
free( (void*)work ); free( (void*)work );
} }
} }
} }