mirror of
https://github.com/TREX-CoE/Sherman-Morrison.git
synced 2024-07-23 03:07:46 +02:00
New implementations for Sherman Morisson
- SM2: Slagel splitting implementation
(http://hdl.handle.net/10919/52966)
- SM3: Keep close to zero denominators for later
- Update tests to show the squared residual
This commit is contained in:
Pablo Oliveira
parent
a68554366c
commit
f096db6411
8
Makefile
8
Makefile
@ -1,11 +1,11 @@
|
||||
## Compilers
|
||||
ARCH =
|
||||
CXX = icpc
|
||||
FC = ifort
|
||||
ARCH =
|
||||
CXX = clang++
|
||||
FC = flang-7
|
||||
H5CXX = h5c++
|
||||
|
||||
## Compiler flags
|
||||
CXXFLAGS = -O0 -debug full -traceback
|
||||
CXXFLAGS = -O0
|
||||
FFLAGS = $(CXXFLAGS)
|
||||
H5CXXFLAGS = $(CXXFLAGS) -fPIC
|
||||
FLIBS = -lstdc++
|
||||
|
||||
@ -4,15 +4,18 @@
|
||||
#include "Helpers.hpp"
|
||||
#include <iostream>
|
||||
|
||||
void SM(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
// Naïve Sherman Morrison
|
||||
void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double *Updates, unsigned int *Updates_index) {
|
||||
|
||||
std::cerr << "Called SM1 with " << N_updates << " updates" << std::endl;
|
||||
double C[Dim];
|
||||
double D[Dim];
|
||||
|
||||
unsigned int l = 0;
|
||||
// For each update
|
||||
for (unsigned int l = 0; l < N_updates; l++) {
|
||||
|
||||
while (l < N_updates) {
|
||||
// C = A^{-1} x U_l
|
||||
for (unsigned int i = 0; i < Dim; i++) {
|
||||
C[i] = 0;
|
||||
@ -23,8 +26,8 @@ void SM(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
// Denominator
|
||||
double den = 1 + C[Updates_index[l] - 1];
|
||||
if (den < 1e-6) {
|
||||
std::cerr << "Breakdown condition triggered" << std::endl;
|
||||
if (fabs(den) < 1e-6) {
|
||||
std::cerr << "Breakdown condition triggered at " << l << std::endl;
|
||||
}
|
||||
double iden = 1 / den;
|
||||
|
||||
@ -40,5 +43,133 @@ void SM(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
Slater_inv[i * Dim + j] -= update;
|
||||
}
|
||||
}
|
||||
|
||||
l += 1;
|
||||
}
|
||||
}
|
||||
|
||||
// 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) {
|
||||
|
||||
std::cerr << "Called SM2 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) < 1e-6) {
|
||||
std::cerr << "Breakdown condition triggered at " << l << std::endl;
|
||||
|
||||
// U_l = U_l / 2 (do the split)
|
||||
for (unsigned int j = 0; j < Dim; j++) {
|
||||
later_updates[later*Dim+j] = Updates[l*Dim+j]/2.0;
|
||||
C[j] /= 2.0;
|
||||
}
|
||||
later_index[later] = Updates_index[l];
|
||||
later++;
|
||||
|
||||
den = 1+C[Updates_index[l]-1];
|
||||
}
|
||||
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 (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) {
|
||||
|
||||
std::cerr << "Called SM2 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) < 1e-6) {
|
||||
std::cerr << "Breakdown condition triggered at " << 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 (later > 0) {
|
||||
SM3(Slater_inv, Dim, later, later_updates, later_index);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
void SM(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
double *Updates, unsigned int *Updates_index) {
|
||||
SM2(Slater_inv, Dim, N_updates, Updates, Updates_index);
|
||||
}
|
||||
|
||||
@ -1,116 +0,0 @@
|
||||
#include <iostream>
|
||||
#include <string>
|
||||
#include "hdf5/serial/hdf5.h"
|
||||
#include "hdf5/serial/H5Cpp.h"
|
||||
|
||||
#include "SM_MaponiA3.hpp"
|
||||
#include "SM_Standard.hpp"
|
||||
#include "Helpers.hpp"
|
||||
|
||||
using namespace H5;
|
||||
//#define DEBUG
|
||||
|
||||
const H5std_string FILE_NAME( "datasets.hdf5" );
|
||||
|
||||
void read_int(H5File file, std::string key, unsigned int * data) {
|
||||
DataSet ds = file.openDataSet(key);
|
||||
ds.read(data, PredType::STD_U32LE);
|
||||
ds.close();
|
||||
}
|
||||
|
||||
void read_double(H5File file, std::string key, double * data) {
|
||||
DataSet ds = file.openDataSet(key);
|
||||
ds.read(data, PredType::IEEE_F64LE);
|
||||
ds.close();
|
||||
}
|
||||
|
||||
int test_cycle(H5File file, int cycle) {
|
||||
|
||||
/* Read the data */
|
||||
|
||||
std::string group = "cycle_" + std::to_string(cycle);
|
||||
|
||||
unsigned int dim, nupdates, col, i, j;
|
||||
read_int(file, group + "/slater_matrix_dim", &dim);
|
||||
read_int(file, group + "/nupdates", &nupdates);
|
||||
|
||||
double * slater_matrix = new double[dim*dim];
|
||||
read_double(file, group + "/slater_matrix", slater_matrix);
|
||||
|
||||
double * slater_inverse = new double[dim*dim];
|
||||
read_double(file, group + "/slater_inverse", slater_inverse);
|
||||
slater_inverse = transpose(slater_inverse, dim);
|
||||
|
||||
unsigned int * col_update_index = new unsigned int[nupdates];
|
||||
read_int(file, group + "/col_update_index", col_update_index);
|
||||
|
||||
double * updates = new double[nupdates*dim];
|
||||
read_double(file, group + "/updates", updates);
|
||||
|
||||
double * u = new double[nupdates*dim];
|
||||
|
||||
|
||||
/* Test */
|
||||
#ifdef DEBUG
|
||||
showMatrix(slater_matrix, dim, "OLD Slater");
|
||||
#endif
|
||||
|
||||
#ifdef DEBUG
|
||||
showMatrix(slater_inverse, dim, "OLD Inverse");
|
||||
#endif
|
||||
|
||||
for (j = 0; j < nupdates; j++) {
|
||||
for (i = 0; i < dim; i++) {
|
||||
col = col_update_index[j];
|
||||
u[i + j*dim] = updates[i + j*dim] - slater_matrix[i*dim + (col - 1)];
|
||||
slater_matrix[i*dim + (col - 1)] = updates[i + j*dim];
|
||||
}
|
||||
}
|
||||
|
||||
//MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index);
|
||||
SM(slater_inverse, dim, nupdates, u, col_update_index);
|
||||
|
||||
#ifdef DEBUG
|
||||
showMatrix(slater_matrix, dim, "NEW Slater");
|
||||
#endif
|
||||
|
||||
#ifdef DEBUG
|
||||
showMatrix(slater_inverse, dim, "NEW Inverse");
|
||||
#endif
|
||||
|
||||
double * res = new double[dim*dim] {0};
|
||||
matMul(slater_matrix, slater_inverse, res, dim);
|
||||
bool ok = is_identity(res, dim, 0.5e-4);
|
||||
|
||||
#ifdef DEBUG
|
||||
showMatrix(res, dim, "Result");
|
||||
#endif
|
||||
|
||||
delete [] res, updates, u, col_update_index,
|
||||
slater_matrix, slater_inverse;
|
||||
|
||||
return ok;
|
||||
}
|
||||
|
||||
int main(int argc, char **argv) {
|
||||
if (argc != 2) {
|
||||
std::cerr << "Execute from within 'datasets/'" << std::endl;
|
||||
std::cerr << "usage: test_external_h5 <cycle>" << std::endl;
|
||||
return 1;
|
||||
}
|
||||
int cycle = std::stoi(argv[1]);
|
||||
H5File file(FILE_NAME, H5F_ACC_RDONLY);
|
||||
|
||||
bool ok = test_cycle(file, cycle);
|
||||
|
||||
if (ok) {
|
||||
std::cerr << "ok -- cycle " << std::to_string(cycle)
|
||||
<< std::endl;
|
||||
}
|
||||
else {
|
||||
std::cerr << "failed -- cycle " << std::to_string(cycle)
|
||||
<< std::endl;
|
||||
}
|
||||
return ok;
|
||||
}
|
||||
|
||||
@ -12,6 +12,18 @@ using namespace H5;
|
||||
|
||||
const H5std_string FILE_NAME( "datasets.hdf5" );
|
||||
|
||||
|
||||
double residual2(double * A, unsigned int Dim) {
|
||||
double res = 0.0;
|
||||
for (unsigned int i = 0; i < Dim; i++) {
|
||||
for (unsigned int j = 0; j < Dim; j++) {
|
||||
double delta = (A[i * Dim + j] - (i == j));
|
||||
res += delta*delta;
|
||||
}
|
||||
}
|
||||
return res;
|
||||
}
|
||||
|
||||
void read_int(H5File file, std::string key, unsigned int * data) {
|
||||
DataSet ds = file.openDataSet(key);
|
||||
ds.read(data, PredType::STD_U32LE);
|
||||
@ -66,8 +78,8 @@ int test_cycle(H5File file, int cycle) {
|
||||
}
|
||||
}
|
||||
|
||||
MaponiA3(slater_inverse, dim, nupdates, u, col_update_index);
|
||||
//SM(slater_inverse, dim, nupdates, updates, col_update_index);
|
||||
//MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index);
|
||||
SM(slater_inverse, dim, nupdates, u, col_update_index);
|
||||
|
||||
#ifdef DEBUG
|
||||
showMatrix(slater_matrix, dim, "NEW Slater");
|
||||
@ -79,7 +91,10 @@ int test_cycle(H5File file, int cycle) {
|
||||
|
||||
double * res = new double[dim*dim] {0};
|
||||
matMul(slater_matrix, slater_inverse, res, dim);
|
||||
bool ok = is_identity(res, dim, 0.5e-4);
|
||||
bool ok = is_identity(res, dim, 1e-4);
|
||||
|
||||
double res2 = residual2(res, dim);
|
||||
std::cerr << "Residual = " << res2 << std::endl;
|
||||
|
||||
#ifdef DEBUG
|
||||
showMatrix(res, dim, "Result");
|
||||
|
||||
Loading…
Reference in New Issue
Block a user