mirror of
https://github.com/TREX-CoE/Sherman-Morrison.git
synced 2024-07-22 18:57:40 +02:00
Merge pull request #25 from fmgjcoppens/testing
- Updated cleaning target in Makefile
This commit is contained in:
François Coppens
GitHub
commit
23dfb25d9a
4
.gitignore
vendored
4
.gitignore
vendored
@ -1,6 +1,10 @@
|
||||
*.o
|
||||
*.mod
|
||||
*.dbg
|
||||
*.cmdx
|
||||
*.cmod
|
||||
*.ilm
|
||||
*.stb
|
||||
.vscode
|
||||
Slater*
|
||||
Updates*
|
||||
|
||||
2
Makefile
2
Makefile
@ -50,7 +50,7 @@ EXEC := $(BIN_DIR)/cMaponiA3_test_3x3_3 \
|
||||
all: $(EXEC)
|
||||
|
||||
clean:
|
||||
@rm -vrf $(OBJ_DIR) *.dbg
|
||||
@rm -vrf $(OBJ_DIR) *.dbg *.cmdx *.cmod *.ilm *.stb
|
||||
|
||||
distclean: clean
|
||||
@rm -vrf $(BIN_DIR) \
|
||||
|
||||
@ -2,12 +2,13 @@
|
||||
|
||||
ENV=$1
|
||||
|
||||
## Set Sherman-Morrison root dir
|
||||
PWD=$(pwd)
|
||||
SRCDIR=$(dirname $BASH_SOURCE)
|
||||
case $SRCDIR in
|
||||
/*) SMROOT=$SRCDIR ;; ## sourced from absolute path
|
||||
*) ## sourced from absolute path
|
||||
if [[ $SRCDIR = . ]] ## check if already in root
|
||||
/*) SMROOT=$SRCDIR ;;
|
||||
*)
|
||||
if [[ $SRCDIR = . ]]
|
||||
then
|
||||
SMROOT=$PWD
|
||||
else
|
||||
@ -17,6 +18,7 @@ case $SRCDIR in
|
||||
esac
|
||||
export SMROOT
|
||||
|
||||
## Set environment for hdf5-tools and Makefile
|
||||
case $ENV in
|
||||
intel)
|
||||
export HDF5_CXX=icpc
|
||||
@ -40,6 +42,7 @@ case $ENV in
|
||||
;;
|
||||
esac
|
||||
|
||||
## Export path, but only once
|
||||
if [[ -z $SMVARS ]]
|
||||
then
|
||||
export PATH=$SMROOT/bin:$PATH
|
||||
|
||||
@ -102,7 +102,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
|
||||
cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "] = " << beta << endl;
|
||||
cout << endl;
|
||||
#endif
|
||||
if (fabs(beta) < 1e-6) {
|
||||
if (fabs(beta) < 1e-3) {
|
||||
cerr << "Break-down occured." << endl;
|
||||
}
|
||||
|
||||
|
||||
@ -4,10 +4,16 @@
|
||||
#include "Helpers.hpp"
|
||||
#include <iostream>
|
||||
|
||||
// Set common break-down threshold
|
||||
double threshold() {
|
||||
double const threshold = 1e-3;
|
||||
return threshold;
|
||||
}
|
||||
|
||||
// Naïve Sherman Morrison
|
||||
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;
|
||||
double C[Dim];
|
||||
@ -15,30 +21,37 @@ void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
unsigned int l = 0;
|
||||
// For each update
|
||||
while (l < N_updates) {
|
||||
while (l < N_updates)
|
||||
{
|
||||
// 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;
|
||||
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];
|
||||
}
|
||||
}
|
||||
|
||||
// Denominator
|
||||
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;
|
||||
}
|
||||
double iden = 1 / den;
|
||||
|
||||
// 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];
|
||||
}
|
||||
|
||||
// A^{-1} = A^{-1} - C x D / den
|
||||
for (unsigned int i = 0; i < Dim; i++) {
|
||||
for (unsigned int j = 0; j < Dim; j++) {
|
||||
for (unsigned int i = 0; i < Dim; i++)
|
||||
{
|
||||
for (unsigned int j = 0; j < Dim; j++)
|
||||
{
|
||||
double update = C[i] * D[j] * iden;
|
||||
Slater_inv[i * Dim + j] -= update;
|
||||
}
|
||||
@ -51,7 +64,8 @@ void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
// Sherman Morrison, with J. Slagel splitting
|
||||
// http://hdl.handle.net/10919/52966
|
||||
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;
|
||||
double C[Dim];
|
||||
@ -63,22 +77,27 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
unsigned int l = 0;
|
||||
// For each update
|
||||
while (l < N_updates) {
|
||||
while (l < N_updates)
|
||||
{
|
||||
// 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;
|
||||
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];
|
||||
}
|
||||
}
|
||||
|
||||
// Denominator
|
||||
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;
|
||||
|
||||
// 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;
|
||||
C[j] /= 2.0;
|
||||
}
|
||||
@ -90,13 +109,16 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double iden = 1 / den;
|
||||
|
||||
// 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];
|
||||
}
|
||||
|
||||
// A^{-1} = A^{-1} - C x D / den
|
||||
for (unsigned int i = 0; i < Dim; i++) {
|
||||
for (unsigned int j = 0; j < Dim; j++) {
|
||||
for (unsigned int i = 0; i < Dim; i++)
|
||||
{
|
||||
for (unsigned int j = 0; j < Dim; j++)
|
||||
{
|
||||
double update = C[i] * D[j] * iden;
|
||||
Slater_inv[i * Dim + j] -= update;
|
||||
}
|
||||
@ -104,14 +126,16 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
l += 1;
|
||||
}
|
||||
|
||||
if (later > 0) {
|
||||
if (later > 0)
|
||||
{
|
||||
SM2(Slater_inv, Dim, later, later_updates, later_index);
|
||||
}
|
||||
}
|
||||
|
||||
// Sherman Morrison, leaving zero denominators for later
|
||||
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;
|
||||
double C[Dim];
|
||||
@ -123,21 +147,26 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
unsigned int l = 0;
|
||||
// For each update
|
||||
while (l < N_updates) {
|
||||
while (l < N_updates)
|
||||
{
|
||||
// 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;
|
||||
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];
|
||||
}
|
||||
}
|
||||
|
||||
// Denominator
|
||||
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;
|
||||
|
||||
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_index[later] = Updates_index[l];
|
||||
@ -148,13 +177,16 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double iden = 1 / den;
|
||||
|
||||
// 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];
|
||||
}
|
||||
|
||||
// A^{-1} = A^{-1} - C x D / den
|
||||
for (unsigned int i = 0; i < Dim; i++) {
|
||||
for (unsigned int j = 0; j < Dim; j++) {
|
||||
for (unsigned int i = 0; i < Dim; i++)
|
||||
{
|
||||
for (unsigned int j = 0; j < Dim; j++)
|
||||
{
|
||||
double update = C[i] * D[j] * iden;
|
||||
Slater_inv[i * Dim + j] -= update;
|
||||
}
|
||||
@ -162,7 +194,8 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
l += 1;
|
||||
}
|
||||
|
||||
if (later > 0) {
|
||||
if (later > 0)
|
||||
{
|
||||
SM3(Slater_inv, Dim, later, later_updates, later_index);
|
||||
}
|
||||
}
|
||||
|
||||
@ -453,6 +453,12 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
|
||||
integer :: i,j
|
||||
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
|
||||
u(i) = m(i) - S(i,l)
|
||||
enddo
|
||||
@ -469,9 +475,9 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
|
||||
d_inv = 1.d0/d
|
||||
d = d+zj
|
||||
lambda = d*d_inv
|
||||
if (dabs(lambda) < 1.d-6) then
|
||||
if (dabs(lambda) < 1.d-3) then
|
||||
d = 0.d0
|
||||
write(502,"('#BREAKDOWN_OCCURED')")
|
||||
write(502,102)
|
||||
return
|
||||
end if
|
||||
|
||||
@ -1168,7 +1174,7 @@ END_PROVIDER
|
||||
|
||||
!! Some usefull formats for output
|
||||
10001 format ('#START_PACKET')
|
||||
10008 format ('#CYCLE_ID: ', I4)
|
||||
10008 format ('#CYCLE_ID: ', I5)
|
||||
10000 format ('#SLATER_MATRIX_DIM: ', I3)
|
||||
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')
|
||||
@ -1177,8 +1183,9 @@ END_PROVIDER
|
||||
10006 format ('#COL_UPDATE_COMP_(',I0.2,'): ', E23.15)
|
||||
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")
|
||||
write(502,'("In subroutine BLD_DET_ALPHA_VALUE_CURR(...)")')
|
||||
|
||||
if (ifirst == 0) then !! If this is the first time we enter this subroutine
|
||||
ifirst = 1
|
||||
@ -1191,7 +1198,6 @@ END_PROVIDER
|
||||
PROVIDE mo_value
|
||||
if (det_i /= det_alpha_order(1) ) then
|
||||
! write(*,*) "det_i: ", det_i
|
||||
! if (det_i == -1 ) then
|
||||
|
||||
!! determin number of updates, the updates and
|
||||
n_to_do = 0
|
||||
@ -1225,38 +1231,6 @@ END_PROVIDER
|
||||
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
|
||||
|
||||
if (n_to_do < shiftl(elec_alpha_num,1)) then
|
||||
@ -1291,6 +1265,7 @@ END_PROVIDER
|
||||
slater_matrix_alpha, &
|
||||
slater_matrix_alpha_inv_det, &
|
||||
ddet)
|
||||
write(502, '("BACK in subroutine BLD_DET_ALPHA_VALUE_CURR(...)")')
|
||||
if (ddet /= 0.d0) then
|
||||
det_alpha_value_curr = ddet
|
||||
else
|
||||
@ -1306,14 +1281,8 @@ END_PROVIDER
|
||||
enddo
|
||||
endif
|
||||
|
||||
write(501,10007)
|
||||
write(502,10007)
|
||||
cycle_id = cycle_id + 1
|
||||
|
||||
else
|
||||
|
||||
ddet = 0.d0
|
||||
|
||||
endif
|
||||
|
||||
|
||||
@ -1337,12 +1306,24 @@ END_PROVIDER
|
||||
slater_matrix_alpha_inv_det(k,i) = mo_value(i,mo_list_alpha_curr(k))
|
||||
enddo
|
||||
enddo
|
||||
call invert(slater_matrix_alpha_inv_det,elec_alpha_num_8,elec_alpha_num,ddet)
|
||||
|
||||
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)
|
||||
endif
|
||||
|
||||
write(501,10007)
|
||||
write(502,10007)
|
||||
cycle_id = cycle_id + 1
|
||||
|
||||
ASSERT (ddet /= 0.d0)
|
||||
|
||||
det_alpha_value_curr = ddet
|
||||
|
||||
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, det_beta_value_curr ]
|
||||
|
||||
41
tests/run_test.sh
Executable file
41
tests/run_test.sh
Executable 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
|
||||
26
tools/stats.py
Executable file
26
tools/stats.py
Executable 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 ("")
|
||||
Loading…
Reference in New Issue
Block a user