TREX/Sherman-Morrison
1
0
Fork 0
mirror of https://github.com/TREX-CoE/Sherman-Morrison.git synced 2024-12-26 06:15:08 +01:00
Code Issues Projects Releases Wiki Activity

Woodbury debugging...

Browse Source
This commit is contained in:
Francois Coppens 2021-06-11 08:46:39 +02:00
parent 61844da5d3
commit 573947fe2d
3 changed files with 17 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
include/Helpers.hpp
View File

@ -117,7 +117,7 @@ template <typename T> void matMul(T *A, T *B, T *C, unsigned int M) {
}
template <typename T1, typename T2, typename T3>
static inline void matMul2(T1 *A, T2 *B, T3 *C, unsigned int M, unsigned int N, unsigned int P) {
void matMul2(T1 *A, T2 *B, T3 *C, unsigned int M, unsigned int N, unsigned int P) {
for(unsigned int i = 0; i < M; i++) {
for(unsigned int j = 0; j < P; j++) {
C[i * P + j] = 0;

14
src/Woodbury.cpp
View File

@ -19,12 +19,26 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
// Construct V from Updates_index
unsigned int V[2 * Dim]; // 2 x Dim matrix stored in row-major order
std::memset(V, 0, 2 * Dim * sizeof(unsigned int));
V[Updates_index[0] - 1] = 1;
V[Dim + Updates_index[1] - 1] = 1;
showMatrix2(Slater_inv, Dim, Dim, "Slater_inv");
showMatrix2(Updates, 2, Dim, "Updates");
showMatrix2(Updates_index, 1, 2, "Updates_index");
showMatrix2(V, 2, Dim, "V");
// Compute C
double C[2 * Dim];
matMul2(Slater_inv, Updates, C, Dim, Dim, 2);
int A[6] = {1,2,3,4,5,6};
int Y[12] = {7,8,9,10,11,12,13,14,15,16,17,18};
int Z[8] = {0,0,0,0,0,0,0,0};
matMul2(A,Y,Z,2,3,4);
showMatrix2(Z, 2,4,"Z");
showMatrix2(C, 2, Dim, "C = S_inv * Updates");
// Compute B
double B[4];
matMul2(V, C, B, 2, Dim, 2);

2
tests/test_h5.cpp
View File

@ -52,6 +52,8 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
showMatrix(slater_inverse, dim, "OLD Inverse");
#endif
showMatrix2(slater_matrix, dim, dim, "Slater");
// Transform replacement updates in 'updates[]' into additive updates in 'u[]'
for (j = 0; j < nupdates; j++) {
for (i = 0; i < dim; i++) {