TREX/Sherman-Morrison
1
0
Fork 0
mirror of https://github.com/TREX-CoE/Sherman-Morrison.git synced 2025-01-13 06:28:35 +01:00
Code Issues Projects Releases Wiki Activity

- Updated cleaning target in Makefile

Browse Source 
- Updated gitignore
- Added python script to compute basic statistical observables from testruns
- Added test bash script to automate test_h5: the cycles in the dataset are not continues and test_h5 crashes if a cycle number is not present in the dataset.
- Added common threshold for all 3 standard SM algos.
This commit is contained in:
Francois Coppens 2021-04-09 17:21:31 +02:00
parent 5175797a6d
commit c2f643ad8f
9 changed files with 175 additions and 87 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
.gitignore vendored
View File

@ -1,6 +1,10 @@
*.o *.o
*.mod *.mod
*.dbg *.dbg
*.cmdx
*.cmod
*.ilm
*.stb
.vscode .vscode
Slater* Slater*
Updates* Updates*

2
Makefile
View File

@ -50,7 +50,7 @@ EXEC := $(BIN_DIR)/cMaponiA3_test_3x3_3 \
all: $(EXEC) all: $(EXEC)
clean: clean:
@rm -vrf $(OBJ_DIR) *.dbg @rm -vrf $(OBJ_DIR) *.dbg *.cmdx *.cmod *.ilm *.stb
distclean: clean distclean: clean
@rm -vrf $(BIN_DIR) \ @rm -vrf $(BIN_DIR) \

9
smvars.sh
View File

@ -2,12 +2,13 @@
ENV=$1 ENV=$1
## Set Sherman-Morrison root dir
PWD=$(pwd) PWD=$(pwd)
SRCDIR=$(dirname $BASH_SOURCE) SRCDIR=$(dirname $BASH_SOURCE)
case $SRCDIR in case $SRCDIR in
/*) SMROOT=$SRCDIR ;; ## sourced from absolute path /*) SMROOT=$SRCDIR ;;
*) ## sourced from absolute path *)
if [[ $SRCDIR = . ]] ## check if already in root if [[ $SRCDIR = . ]]
then then
SMROOT=$PWD SMROOT=$PWD
else else
@ -17,6 +18,7 @@ case $SRCDIR in
esac esac
export SMROOT export SMROOT
## Set environment for hdf5-tools and Makefile
case $ENV in case $ENV in
intel) intel)
export HDF5_CXX=icpc export HDF5_CXX=icpc
@ -40,6 +42,7 @@ case $ENV in
;; ;;
esac esac
## Export path, but only once
if [[ -z $SMVARS ]] if [[ -z $SMVARS ]]
then then
export PATH=$SMROOT/bin:$PATH export PATH=$SMROOT/bin:$PATH

2
src/SM_MaponiA3.cpp
View File

@ -102,7 +102,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "] = " << beta << endl; cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "] = " << beta << endl;
cout << endl; cout << endl;
#endif #endif
if (fabs(beta) < 1e-6) { if (fabs(beta) < 1e-3) {
cerr << "Break-down occured." << endl; cerr << "Break-down occured." << endl;
} }

105
src/SM_Standard.cpp
View File

@ -4,10 +4,16 @@
#include "Helpers.hpp" #include "Helpers.hpp"
#include <iostream> #include <iostream>
// Set common break-down threshold
double threshold() {
double const threshold = 1e-3;
return threshold;
}
// Naïve Sherman Morrison // Naïve Sherman Morrison
void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
double *Updates, unsigned int *Updates_index) { double *Updates, unsigned int *Updates_index)
{
std::cerr << "Called SM1 with " << N_updates << " updates" << std::endl; std::cerr << "Called SM1 with " << N_updates << " updates" << std::endl;
double C[Dim]; double C[Dim];
@ -15,30 +21,37 @@ void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
unsigned int l = 0; unsigned int l = 0;
// For each update // For each update
while (l < N_updates) { while (l < N_updates)
{
// C = A^{-1} x U_l // C = A^{-1} x U_l
for (unsigned int i = 0; i < Dim; i++) { for (unsigned int i = 0; i < Dim; i++)
{
C[i] = 0; C[i] = 0;
for (unsigned int j = 0; j < Dim; j++) { for (unsigned int j = 0; j < Dim; j++)
{
C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j]; C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j];
} }
} }
// Denominator // Denominator
double den = 1 + C[Updates_index[l] - 1]; double den = 1 + C[Updates_index[l] - 1];
if (fabs(den) < 1e-6) { if (fabs(den) < threshold())
std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl; {
std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl;
} }
double iden = 1 / den; double iden = 1 / den;
// D = v^T x A^{-1} // D = v^T x A^{-1}
for (unsigned int j = 0; j < Dim; j++) { for (unsigned int j = 0; j < Dim; j++)
{
D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j]; D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];
} }
// A^{-1} = A^{-1} - C x D / den // A^{-1} = A^{-1} - C x D / den
for (unsigned int i = 0; i < Dim; i++) { for (unsigned int i = 0; i < Dim; i++)
for (unsigned int j = 0; j < Dim; j++) { {
for (unsigned int j = 0; j < Dim; j++)
{
double update = C[i] * D[j] * iden; double update = C[i] * D[j] * iden;
Slater_inv[i * Dim + j] -= update; Slater_inv[i * Dim + j] -= update;
} }
@ -51,52 +64,61 @@ void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
// Sherman Morrison, with J. Slagel splitting // Sherman Morrison, with J. Slagel splitting
// http://hdl.handle.net/10919/52966 // http://hdl.handle.net/10919/52966
void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates, void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
double *Updates, unsigned int *Updates_index) { double *Updates, unsigned int *Updates_index)
{
std::cerr << "Called SM2 with updates " << N_updates << std::endl; std::cerr << "Called SM2 with updates " << N_updates << std::endl;
double C[Dim]; double C[Dim];
double D[Dim]; double D[Dim];
double later_updates[Dim*N_updates]; double later_updates[Dim * N_updates];
unsigned int later_index[N_updates]; unsigned int later_index[N_updates];
unsigned int later = 0; unsigned int later = 0;
unsigned int l = 0; unsigned int l = 0;
// For each update // For each update
while (l < N_updates) { while (l < N_updates)
{
// C = A^{-1} x U_l // C = A^{-1} x U_l
for (unsigned int i = 0; i < Dim; i++) { for (unsigned int i = 0; i < Dim; i++)
{
C[i] = 0; C[i] = 0;
for (unsigned int j = 0; j < Dim; j++) { for (unsigned int j = 0; j < Dim; j++)
{
C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j]; C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j];
} }
} }
// Denominator // Denominator
double den = 1 + C[Updates_index[l] - 1]; double den = 1 + C[Updates_index[l] - 1];
if (fabs(den) < 1e-6) { if (fabs(den) < threshold())
std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl; {
std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl;
// U_l = U_l / 2 (do the split) // U_l = U_l / 2 (do the split)
for (unsigned int j = 0; j < Dim; j++) { for (unsigned int j = 0; j < Dim; j++)
later_updates[later*Dim+j] = Updates[l*Dim+j]/2.0; {
later_updates[later * Dim + j] = Updates[l * Dim + j] / 2.0;
C[j] /= 2.0; C[j] /= 2.0;
} }
later_index[later] = Updates_index[l]; later_index[later] = Updates_index[l];
later++; later++;
den = 1+C[Updates_index[l]-1]; den = 1 + C[Updates_index[l] - 1];
} }
double iden = 1 / den; double iden = 1 / den;
// D = v^T x A^{-1} // D = v^T x A^{-1}
for (unsigned int j = 0; j < Dim; j++) { for (unsigned int j = 0; j < Dim; j++)
{
D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j]; D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];
} }
// A^{-1} = A^{-1} - C x D / den // A^{-1} = A^{-1} - C x D / den
for (unsigned int i = 0; i < Dim; i++) { for (unsigned int i = 0; i < Dim; i++)
for (unsigned int j = 0; j < Dim; j++) { {
for (unsigned int j = 0; j < Dim; j++)
{
double update = C[i] * D[j] * iden; double update = C[i] * D[j] * iden;
Slater_inv[i * Dim + j] -= update; Slater_inv[i * Dim + j] -= update;
} }
@ -104,41 +126,48 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
l += 1; l += 1;
} }
if (later > 0) { if (later > 0)
{
SM2(Slater_inv, Dim, later, later_updates, later_index); SM2(Slater_inv, Dim, later, later_updates, later_index);
} }
} }
// Sherman Morrison, leaving zero denominators for later // Sherman Morrison, leaving zero denominators for later
void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
double *Updates, unsigned int *Updates_index) { double *Updates, unsigned int *Updates_index)
{
std::cerr << "Called SM3 with updates " << N_updates << std::endl; std::cerr << "Called SM3 with updates " << N_updates << std::endl;
double C[Dim]; double C[Dim];
double D[Dim]; double D[Dim];
double later_updates[Dim*N_updates]; double later_updates[Dim * N_updates];
unsigned int later_index[N_updates]; unsigned int later_index[N_updates];
unsigned int later = 0; unsigned int later = 0;
unsigned int l = 0; unsigned int l = 0;
// For each update // For each update
while (l < N_updates) { while (l < N_updates)
{
// C = A^{-1} x U_l // C = A^{-1} x U_l
for (unsigned int i = 0; i < Dim; i++) { for (unsigned int i = 0; i < Dim; i++)
{
C[i] = 0; C[i] = 0;
for (unsigned int j = 0; j < Dim; j++) { for (unsigned int j = 0; j < Dim; j++)
{
C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j]; C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j];
} }
} }
// Denominator // Denominator
double den = 1 + C[Updates_index[l] - 1]; double den = 1 + C[Updates_index[l] - 1];
if (fabs(den) < 1e-6) { if (fabs(den) < threshold())
std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl; {
std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl;
for (unsigned int j = 0; j < Dim; j++) { for (unsigned int j = 0; j < Dim; j++)
later_updates[later*Dim+j] = Updates[l*Dim+j]; {
later_updates[later * Dim + j] = Updates[l * Dim + j];
} }
later_index[later] = Updates_index[l]; later_index[later] = Updates_index[l];
later++; later++;
@ -148,13 +177,16 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
double iden = 1 / den; double iden = 1 / den;
// D = v^T x A^{-1} // D = v^T x A^{-1}
for (unsigned int j = 0; j < Dim; j++) { for (unsigned int j = 0; j < Dim; j++)
{
D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j]; D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];
} }
// A^{-1} = A^{-1} - C x D / den // A^{-1} = A^{-1} - C x D / den
for (unsigned int i = 0; i < Dim; i++) { for (unsigned int i = 0; i < Dim; i++)
for (unsigned int j = 0; j < Dim; j++) { {
for (unsigned int j = 0; j < Dim; j++)
{
double update = C[i] * D[j] * iden; double update = C[i] * D[j] * iden;
Slater_inv[i * Dim + j] -= update; Slater_inv[i * Dim + j] -= update;
} }
@ -162,7 +194,8 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
l += 1; l += 1;
} }
if (later > 0) { if (later > 0)
{
SM3(Slater_inv, Dim, later, later_updates, later_index); SM3(Slater_inv, Dim, later, later_updates, later_index);
} }
} }

71
src/det.irp.f
View File

@ -452,7 +452,13 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
!DIR$ ASSUME (LDS >= $n) !DIR$ ASSUME (LDS >= $n)
integer :: i,j integer :: i,j
double precision :: zj, zj1, zj2, zj3 double precision :: zj, zj1, zj2, zj3
101 format ('Updating column ', I2)
102 format ('#BREAKDOWN_OCCURED: TH=1.d-3')
write(502, '("In subroutine DET_UPDATE21()")')
write(502,101) l
do i=1,$n do i=1,$n
u(i) = m(i) - S(i,l) u(i) = m(i) - S(i,l)
enddo enddo
@ -469,9 +475,9 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
d_inv = 1.d0/d d_inv = 1.d0/d
d = d+zj d = d+zj
lambda = d*d_inv lambda = d*d_inv
if (dabs(lambda) < 1.d-6) then if (dabs(lambda) < 1.d-3) then
d = 0.d0 d = 0.d0
write(502,"('#BREAKDOWN_OCCURED')") write(502,102)
return return
end if end if
@ -970,7 +976,7 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d)
!DIR$ ASSUME (MOD(LDS,$IRP_ALIGN/8) == 0) !DIR$ ASSUME (MOD(LDS,$IRP_ALIGN/8) == 0)
!DIR$ ASSUME (n>150) !DIR$ ASSUME (n>150)
integer :: i,j,n4 integer :: i,j,n4
double precision :: zl double precision :: zl
!DIR$ NOPREFETCH !DIR$ NOPREFETCH
@ -1168,7 +1174,7 @@ END_PROVIDER
!! Some usefull formats for output !! Some usefull formats for output
10001 format ('#START_PACKET') 10001 format ('#START_PACKET')
10008 format ('#CYCLE_ID: ', I4) 10008 format ('#CYCLE_ID: ', I5)
10000 format ('#SLATER_MATRIX_DIM: ', I3) 10000 format ('#SLATER_MATRIX_DIM: ', I3)
10002 format ('#NUPDATES: ', I2) 10002 format ('#NUPDATES: ', I2)
10003 format ('#SLATER_MATRIX: (i (outer), j (inner)), slater_matrix_alpha(i,j), ddet * slater_matrix_alpha_inv_det(i,j) / ddet') 10003 format ('#SLATER_MATRIX: (i (outer), j (inner)), slater_matrix_alpha(i,j), ddet * slater_matrix_alpha_inv_det(i,j) / ddet')
@ -1177,8 +1183,9 @@ END_PROVIDER
10006 format ('#COL_UPDATE_COMP_(',I0.2,'): ', E23.15) 10006 format ('#COL_UPDATE_COMP_(',I0.2,'): ', E23.15)
10007 format ('#END_PACKET',/) 10007 format ('#END_PACKET',/)
open (unit = 501, file = "dataset.dat") !! slightly cleaner output open (unit = 501, file = "dataset.dat") !! Only Slater matrices plus updates
open (unit = 502, file = "dataset.fulltrace.dat") open (unit = 502, file = "dataset.fulltrace.dat")
write(502,'("In subroutine BLD_DET_ALPHA_VALUE_CURR(...)")')
if (ifirst == 0) then !! If this is the first time we enter this subroutine if (ifirst == 0) then !! If this is the first time we enter this subroutine
ifirst = 1 ifirst = 1
@ -1191,7 +1198,6 @@ END_PROVIDER
PROVIDE mo_value PROVIDE mo_value
if (det_i /= det_alpha_order(1) ) then if (det_i /= det_alpha_order(1) ) then
! write(*,*) "det_i: ", det_i ! write(*,*) "det_i: ", det_i
! if (det_i == -1 ) then
!! determin number of updates, the updates and !! determin number of updates, the updates and
n_to_do = 0 n_to_do = 0
@ -1225,38 +1231,6 @@ END_PROVIDER
end do end do
end do end do
! print n_to_do (number of updates to do)
! to_do (array of the columns that need to be swapped)
! mo_list_alpha_curr (list of orbitals to build the current determinant)
! mo_list_alpha_prev (list of the previous determinant)
!
! slater_matrix_alpha (slater matrix that needs to be inverted. This is the initial matrix)
! slater_matrix_alpha_inv_det = inverse of the slater matrix divided by the determinant (ddet: l. 1216)
!
!
!
! 1 2 3 4
! --------
! 2 4 6 8
! 2 4 10 8
!
! mo_list(:) (2,4,10,8)
! mo_list_prev(:) (2,4,6,8)
! n_todo 1
! to_do (3)
!
!
! 1 2 3 4
! --------
! 2 4 6 8
! 2 5 10 8
!
! mo_list(:) (2,5,10,8)
! mo_list_prev(:) (2,4,6,8)
! n_todo 2
! to_do (2,3)
ddet = 0.d0 ddet = 0.d0
if (n_to_do < shiftl(elec_alpha_num,1)) then if (n_to_do < shiftl(elec_alpha_num,1)) then
@ -1291,6 +1265,7 @@ END_PROVIDER
slater_matrix_alpha, & slater_matrix_alpha, &
slater_matrix_alpha_inv_det, & slater_matrix_alpha_inv_det, &
ddet) ddet)
write(502, '("BACK in subroutine BLD_DET_ALPHA_VALUE_CURR(...)")')
if (ddet /= 0.d0) then if (ddet /= 0.d0) then
det_alpha_value_curr = ddet det_alpha_value_curr = ddet
else else
@ -1305,15 +1280,9 @@ END_PROVIDER
endif endif
enddo enddo
endif endif
write(501,10007)
write(502,10007)
cycle_id = cycle_id + 1
else else
ddet = 0.d0 ddet = 0.d0
endif endif
@ -1337,12 +1306,24 @@ END_PROVIDER
slater_matrix_alpha_inv_det(k,i) = mo_value(i,mo_list_alpha_curr(k)) slater_matrix_alpha_inv_det(k,i) = mo_value(i,mo_list_alpha_curr(k))
enddo enddo
enddo enddo
if (det_i /= det_alpha_order(1)) then
write(502,"('#NO_INVERSE_FOUND: ')") cycle_id
else
write(502,"('#FIRST_DETERMINANT_IN_SUPERCYCLE')")
endif
call invert(slater_matrix_alpha_inv_det,elec_alpha_num_8,elec_alpha_num,ddet) call invert(slater_matrix_alpha_inv_det,elec_alpha_num_8,elec_alpha_num,ddet)
endif endif
write(501,10007)
write(502,10007)
cycle_id = cycle_id + 1
ASSERT (ddet /= 0.d0) ASSERT (ddet /= 0.d0)
det_alpha_value_curr = ddet det_alpha_value_curr = ddet
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, det_beta_value_curr ] BEGIN_PROVIDER [ double precision, det_beta_value_curr ]

41
tests/run_test.sh Executable file
View File

@ -0,0 +1,41 @@
#!/usr/bin/env bash
## IMPORTANT: start cycle will always be 1 because of the intention of this
## script and the way it is written now. In fact, var. $START is not used at
## the moment.
## Call: $ ./run_test.sh <{sm1|sm2|sm3|maponia3}> <start cycle> <stop cycle>
## Example: ./run_test.sh sm2 1 500
TEST=test_h5
ALGO=$1
START=$2
STOP=$3
BLOCKSIZE=329
BLOCKS=$((STOP / BLOCKSIZE))
LEFT=$((STOP % BLOCKSIZE))
if [[ $LEFT -ne 0 ]]
then
BSTART=0
BSTOP=0
LAST=0
for ((i=1; i<=$BLOCKS; i++))
do
BSTART=$(((i-1)*BLOCKSIZE+1))
BSTOP=$((i*BLOCKSIZE-1))
$TEST $ALGO $BSTART $BSTOP
LAST=$i
done
LSTART=$((LAST*BLOCKSIZE+1))
LSTOP=$((LSTART+LEFT-1))
$TEST $ALGO $LSTART $LSTOP
else
for ((i=1; i<=$BLOCKS; i++))
do
BSTART=$(((i-1)*BLOCKSIZE+1))
BSTOP=$((i*BLOCKSIZE-1))
$TEST $ALGO $BSTART $BSTOP
done
fi

2
tests/test_h5.cpp
View File

@ -13,7 +13,7 @@ using namespace H5;
const H5std_string FILE_NAME( "dataset.hdf5" ); const H5std_string FILE_NAME( "dataset.hdf5" );
double residual_max(double * A, unsigned int Dim) { double residual_max(double * A, unsigned int Dim) {
double max= 0.0; double max = 0.0;
for (unsigned int i = 0; i < Dim; i++) { for (unsigned int i = 0; i < Dim; i++) {
for (unsigned int j = 0; j < Dim; j++) { for (unsigned int j = 0; j < Dim; j++) {
double delta = (A[i * Dim + j] - (i == j)); double delta = (A[i * Dim + j] - (i == j));

26
tools/stats.py Executable file
View File

@ -0,0 +1,26 @@
#!/usr/bin/env python
import sys
import numpy as np
filename = sys.argv[1]
Cycle, Rmax, R2 = np.loadtxt(filename, usecols = (3, 4, 5), unpack = True)
print ("========================================================")
print ("Statistics for file: ", filename)
print ("========================================================")
print ("Max element residual")
print ("--------------------")
print ("MINIMUM, Cycle#: ", np.nanmin(Rmax), ", ", Cycle[np.nanargmin(Rmax)])
print ("MAXIMUM, Cycle#: ", np.nanmax(Rmax), ", " , Cycle[np.nanargmax(Rmax)])
print ("MEAN: ", np.nanmean(Rmax))
print ("STD: ", np.nanstd(Rmax))
print ("")
print ("Frobenius norm squared residual")
print ("--------------------------------")
print ("MINIMUM, Cycle#: ", np.nanmin(R2), ", ", Cycle[np.nanargmin(R2)])
print ("MAXIMUM, Cycle#: ", np.nanmax(R2), ", " , Cycle[np.nanargmax(R2)])
print ("MEAN: ", np.nanmean(R2))
print ("STD: ", np.nanstd(R2))
print ("")