mirror of
https://github.com/TREX-CoE/Sherman-Morrison.git
synced 2024-08-25 14:41:45 +02:00
Merge pull request #32 from pablooliveira/sm4
First implementation of SM4
This commit is contained in:
commit
8595a598cc
@ -10,3 +10,8 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
|||||||
// 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);
|
||||||
|
|
||||||
|
// Sherman Morrison (SM3+SM2), leaving zero denominators for later (SM3), and
|
||||||
|
// when none are left falling back on Splitting (SM2)
|
||||||
|
void SM4(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||||
|
double *Updates, unsigned int *Updates_index);
|
||||||
|
@ -172,6 +172,70 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
// Sherman Morrison, mix between SM3 + SM2
|
||||||
|
// Leave zero denominators for later (SM3), and when none are left then split (SM2)
|
||||||
|
void SM4(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||||
|
double *Updates, unsigned int *Updates_index) {
|
||||||
|
std::cerr << "Called SM4 with updates " << N_updates << std::endl;
|
||||||
|
double C[Dim];
|
||||||
|
double D[Dim];
|
||||||
|
|
||||||
|
double later_updates[Dim * N_updates];
|
||||||
|
unsigned int later_index[N_updates];
|
||||||
|
unsigned int later = 0;
|
||||||
|
|
||||||
|
unsigned int l = 0;
|
||||||
|
// For each update
|
||||||
|
while (l < N_updates) {
|
||||||
|
// C = A^{-1} x U_l
|
||||||
|
for (unsigned int i = 0; i < Dim; i++) {
|
||||||
|
C[i] = 0;
|
||||||
|
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) < threshold()) {
|
||||||
|
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||||
|
<< std::endl;
|
||||||
|
|
||||||
|
for (unsigned int j = 0; j < Dim; j++) {
|
||||||
|
later_updates[later * Dim + j] = Updates[l * Dim + j];
|
||||||
|
}
|
||||||
|
later_index[later] = Updates_index[l];
|
||||||
|
later++;
|
||||||
|
l += 1;
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
double iden = 1 / den;
|
||||||
|
|
||||||
|
// D = v^T x A^{-1}
|
||||||
|
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++) {
|
||||||
|
double update = C[i] * D[j] * iden;
|
||||||
|
Slater_inv[i * Dim + j] -= update;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
l += 1;
|
||||||
|
}
|
||||||
|
|
||||||
|
// If all the updates have failed, fall back on splitting (SM2)
|
||||||
|
if (later == N_updates) {
|
||||||
|
SM2(Slater_inv, Dim, later, later_updates, later_index);
|
||||||
|
}
|
||||||
|
// If some have failed, make a recursive call
|
||||||
|
else if (later > 0) {
|
||||||
|
SM4(Slater_inv, Dim, later, later_updates, later_index);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
extern "C" {
|
extern "C" {
|
||||||
void SM1_f(double **linSlater_inv, unsigned int *Dim,
|
void SM1_f(double **linSlater_inv, unsigned int *Dim,
|
||||||
unsigned int *N_updates, double **linUpdates,
|
unsigned int *N_updates, double **linUpdates,
|
||||||
@ -190,4 +254,9 @@ extern "C" {
|
|||||||
unsigned int **Updates_index) {
|
unsigned int **Updates_index) {
|
||||||
SM3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
SM3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
||||||
}
|
}
|
||||||
|
void SM4_f(double **linSlater_inv, unsigned int *Dim,
|
||||||
|
unsigned int *N_updates, double **linUpdates,
|
||||||
|
unsigned int **Updates_index) {
|
||||||
|
SM4(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
@ -36,5 +36,12 @@ module Sherman_Morrison
|
|||||||
real(c_double), dimension(:,:), allocatable, intent(in) :: Updates
|
real(c_double), dimension(:,:), allocatable, intent(in) :: Updates
|
||||||
real(c_double), dimension(:,:), allocatable, intent(in out) :: Slater_inv
|
real(c_double), dimension(:,:), allocatable, intent(in out) :: Slater_inv
|
||||||
end subroutine SM3
|
end subroutine SM3
|
||||||
|
subroutine SM4(Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="SM4_f")
|
||||||
|
use, intrinsic :: iso_c_binding, only : c_int, c_double
|
||||||
|
integer(c_int), intent(in) :: dim, n_updates
|
||||||
|
integer(c_int), dimension(:), allocatable, intent(in) :: Updates_index
|
||||||
|
real(c_double), dimension(:,:), allocatable, intent(in) :: Updates
|
||||||
|
real(c_double), dimension(:,:), allocatable, intent(in out) :: Slater_inv
|
||||||
|
end subroutine SM4
|
||||||
end interface
|
end interface
|
||||||
end module Sherman_Morrison
|
end module Sherman_Morrison
|
||||||
|
@ -104,6 +104,8 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
|
|||||||
SM2(slater_inverse, dim, nupdates, u, col_update_index);
|
SM2(slater_inverse, dim, nupdates, u, col_update_index);
|
||||||
} else if (version == "sm3") {
|
} else if (version == "sm3") {
|
||||||
SM3(slater_inverse, dim, nupdates, u, col_update_index);
|
SM3(slater_inverse, dim, nupdates, u, col_update_index);
|
||||||
|
} else if (version == "sm4") {
|
||||||
|
SM4(slater_inverse, dim, nupdates, u, col_update_index);
|
||||||
} else {
|
} else {
|
||||||
std::cerr << "Unknown version " << version << std::endl;
|
std::cerr << "Unknown version " << version << std::endl;
|
||||||
exit(1);
|
exit(1);
|
||||||
|
@ -137,6 +137,8 @@ int test_cycle(H5File file, int cycle, std::string version, vfc_probes * probes)
|
|||||||
SM2(slater_inverse, dim, nupdates, u, col_update_index);
|
SM2(slater_inverse, dim, nupdates, u, col_update_index);
|
||||||
} else if (version == "sm3") {
|
} else if (version == "sm3") {
|
||||||
SM3(slater_inverse, dim, nupdates, u, col_update_index);
|
SM3(slater_inverse, dim, nupdates, u, col_update_index);
|
||||||
|
} else if (version == "sm4") {
|
||||||
|
SM4(slater_inverse, dim, nupdates, u, col_update_index);
|
||||||
} else {
|
} else {
|
||||||
std::cerr << "Unknown version " << version << std::endl;
|
std::cerr << "Unknown version " << version << std::endl;
|
||||||
exit(1);
|
exit(1);
|
||||||
|
Loading…
Reference in New Issue
Block a user