mirror of
https://github.com/TREX-CoE/Sherman-Morrison.git
synced 2024-12-26 06:15:08 +01:00
commit
ab661ad785
@ -82,58 +82,107 @@ T1 *outProd(T1 *vec1, T2 *vec2, unsigned int M) {
|
||||
return C;
|
||||
}
|
||||
|
||||
template <typename T> T matDet(T **A, unsigned int M) {
|
||||
int det = 0, p, h, k, i, j;
|
||||
T **temp = new T *[M];
|
||||
for (int i = 0; i < M; i++)
|
||||
temp[i] = new T[M];
|
||||
if (M == 1) {
|
||||
return A[0][0];
|
||||
} else if (M == 2) {
|
||||
det = (A[0][0] * A[1][1] - A[0][1] * A[1][0]);
|
||||
return det;
|
||||
} else {
|
||||
for (p = 0; p < M; p++) {
|
||||
h = 0;
|
||||
k = 0;
|
||||
for (i = 1; i < M; i++) {
|
||||
for (j = 0; j < M; j++) {
|
||||
if (j == p) {
|
||||
continue;
|
||||
}
|
||||
temp[h][k] = A[i][j];
|
||||
k++;
|
||||
if (k == M - 1) {
|
||||
h++;
|
||||
k = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
det = det + A[0][p] * pow(-1, p) * matDet(temp, M - 1);
|
||||
}
|
||||
return det;
|
||||
}
|
||||
delete[] temp;
|
||||
}
|
||||
// // This flat version doesn't work. Get's stuck in an infinite recursion loop.
|
||||
// template <typename T> T determinant(T *A, unsigned int M) {
|
||||
// std::cout << "determinant() called..." << std::endl;
|
||||
// T det = 0;
|
||||
// int p, h, k, i, j;
|
||||
// T *temp = new T[M * M];
|
||||
// if (M == 1) {
|
||||
// return A[0];
|
||||
// } else if (M == 2) {
|
||||
// det = (A[0] * A[3] - A[1] * A[2]);
|
||||
// return det;
|
||||
// } else {
|
||||
// for (p = 0; p < M; p++) {
|
||||
// h = 0;
|
||||
// k = 0;
|
||||
// for (i = 1; i < M; i++) {
|
||||
// for (j = 0; j < M; j++) {
|
||||
// if (j == p) {
|
||||
// continue;
|
||||
// }
|
||||
// temp[h * M + k] = A[i * M + j];
|
||||
// k++;
|
||||
// if (k == M - 1) {
|
||||
// h++;
|
||||
// k = 0;
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// det = det + A[p] * pow(-1, p) * determinant(temp, M - 1);
|
||||
// }
|
||||
// return det;
|
||||
// }
|
||||
// delete temp;
|
||||
// }
|
||||
|
||||
// // This version also gets stuck in a recursion loop
|
||||
// template <typename T> T determinant(T **A, unsigned int M) {
|
||||
// int p, h, k, i, j;
|
||||
// T det = 0;
|
||||
// T **temp = new T *[M];
|
||||
// for (int i = 0; i < M; i++) {
|
||||
// temp[i] = new T[M];
|
||||
// }
|
||||
// if (M == 1) {
|
||||
// return A[0][0];
|
||||
// } else if (M == 2) {
|
||||
// det = (A[0][0] * A[1][1] - A[0][1] * A[1][0]);
|
||||
// return det;
|
||||
// } else {
|
||||
// for (p = 0; p < M; p++) {
|
||||
// h = 0;
|
||||
// k = 0;
|
||||
// for (i = 1; i < M; i++) {
|
||||
// for (j = 0; j < M; j++) {
|
||||
// if (j == p) {
|
||||
// continue;
|
||||
// }
|
||||
// temp[h][k] = A[i][j];
|
||||
// k++;
|
||||
// if (k == M - 1) {
|
||||
// h++;
|
||||
// k = 0;
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// det = det + A[0][p] * pow(-1, p) * determinant(temp, M - 1);
|
||||
// }
|
||||
// return det;
|
||||
// }
|
||||
// delete[] temp;
|
||||
// }
|
||||
|
||||
template <typename T> bool is_identity(T *A, unsigned int M, double tolerance) {
|
||||
for (unsigned int i = 0; i < M; i++) {
|
||||
for (unsigned int j = 0; j < M; j++) {
|
||||
if (i == j && fabs(A[i * M + j] - 1) > tolerance)
|
||||
if (i == j && std::fabs(A[i * M + j] - 1) > tolerance) {
|
||||
return false;
|
||||
if (i != j && fabs(A[i * M + j]) > tolerance)
|
||||
}
|
||||
if (i != j && std::fabs(A[i * M + j]) > tolerance) {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
template <typename T> T norm_max(T * A, unsigned int Dim) {
|
||||
template <typename T>
|
||||
bool is_identity2(T *A, unsigned int M, double tolerance) {
|
||||
double det = determinant(A, M);
|
||||
if (det - 1 > tolerance) {
|
||||
return false;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
template <typename T> T norm_max(T *A, unsigned int Dim) {
|
||||
T res = 0;
|
||||
for (unsigned int i = 0; i < Dim; i++) {
|
||||
for (unsigned int j = 0; j < Dim; j++) {
|
||||
T delta = A[i * Dim + j];
|
||||
delta = fabs(delta);
|
||||
delta = std::fabs(delta);
|
||||
if (delta > res) {
|
||||
res = delta;
|
||||
}
|
||||
@ -142,23 +191,23 @@ template <typename T> T norm_max(T * A, unsigned int Dim) {
|
||||
return res;
|
||||
}
|
||||
|
||||
template <typename T> T norm_frobenius2(T * A, unsigned int Dim) {
|
||||
template <typename T> T norm_frobenius2(T *A, unsigned int Dim) {
|
||||
T res = 0;
|
||||
for (unsigned int i = 0; i < Dim; i++) {
|
||||
for (unsigned int j = 0; j < Dim; j++) {
|
||||
T delta = A[i * Dim + j];
|
||||
res += delta*delta;
|
||||
res += delta * delta;
|
||||
}
|
||||
}
|
||||
return res;
|
||||
}
|
||||
|
||||
template <typename T> T residual_max(T * A, unsigned int Dim) {
|
||||
template <typename T> T residual_max(T *A, unsigned int Dim) {
|
||||
T res = 0;
|
||||
for (unsigned int i = 0; i < Dim; i++) {
|
||||
for (unsigned int j = 0; j < Dim; j++) {
|
||||
T delta = A[i * Dim + j] - (i == j);
|
||||
delta = fabs(delta);
|
||||
delta = std::fabs(delta);
|
||||
if (delta > res) {
|
||||
res = delta;
|
||||
}
|
||||
@ -167,12 +216,12 @@ template <typename T> T residual_max(T * A, unsigned int Dim) {
|
||||
return res;
|
||||
}
|
||||
|
||||
template <typename T> T residual_frobenius2(T * A, unsigned int Dim) {
|
||||
template <typename T> T residual_frobenius2(T *A, unsigned int Dim) {
|
||||
T res = 0;
|
||||
for (unsigned int i = 0; i < Dim; i++) {
|
||||
for (unsigned int j = 0; j < Dim; j++) {
|
||||
T delta = A[i * Dim + j] - (i == j);
|
||||
res += delta*delta;
|
||||
res += delta * delta;
|
||||
}
|
||||
}
|
||||
return res;
|
||||
|
@ -87,7 +87,7 @@ void vfc_hashmap_destroy(vfc_hashmap_t map) {
|
||||
|
||||
// allocate and initialize the map
|
||||
vfc_hashmap_t vfc_hashmap_create() {
|
||||
vfc_hashmap_t map = (vfc_hashmap_t) calloc(1, sizeof(struct vfc_hashmap_st));
|
||||
vfc_hashmap_t map = (vfc_hashmap_t)calloc(1, sizeof(struct vfc_hashmap_st));
|
||||
|
||||
if (map == NULL) {
|
||||
return NULL;
|
||||
@ -96,7 +96,7 @@ vfc_hashmap_t vfc_hashmap_create() {
|
||||
map->capacity = (size_t)(1 << map->nbits);
|
||||
map->mask = map->capacity - 1;
|
||||
// an item is now a value and a key
|
||||
map->items = (size_t *) calloc(map->capacity, 2 * sizeof(size_t));
|
||||
map->items = (size_t *)calloc(map->capacity, 2 * sizeof(size_t));
|
||||
if (map->items == NULL) {
|
||||
vfc_hashmap_destroy(map);
|
||||
return NULL;
|
||||
@ -160,7 +160,7 @@ static void maybe_rehash_map(vfc_hashmap_t map) {
|
||||
map->nbits++;
|
||||
map->capacity = (size_t)(1 << map->nbits);
|
||||
map->mask = map->capacity - 1;
|
||||
map->items = (size_t *) calloc(map->capacity, 2 * sizeof(size_t));
|
||||
map->items = (size_t *)calloc(map->capacity, 2 * sizeof(size_t));
|
||||
map->nitems = 0;
|
||||
map->n_deleted_items = 0;
|
||||
for (ii = 0; ii < old_capacity; ii++) {
|
||||
|
@ -1,37 +1,33 @@
|
||||
/*
|
||||
* This file defines "vfc_probes", a hashtable-based structure which can be used
|
||||
* to place "probes" in a code and store the different values of test variables.
|
||||
* These test results can then be exported in a CSV file, and used to generate a
|
||||
* Verificarlo test report.
|
||||
*/
|
||||
|
||||
* This file defines "vfc_probes", a hashtable-based structure which can be used
|
||||
* to place "probes" in a code and store the different values of test variables.
|
||||
* These test results can then be exported in a CSV file, and used to generate a
|
||||
* Verificarlo test report.
|
||||
*/
|
||||
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <stdio.h>
|
||||
|
||||
#include "vfc_hashmap.h"
|
||||
|
||||
#define VAR_NAME(var) #var // Simply returns the name of var into a string
|
||||
|
||||
|
||||
/*
|
||||
* A probe containing a double value as well as its key, which is needed when
|
||||
* dumping the probes
|
||||
*/
|
||||
* A probe containing a double value as well as its key, which is needed when
|
||||
* dumping the probes
|
||||
*/
|
||||
|
||||
struct vfc_probe_node {
|
||||
char * key;
|
||||
char *key;
|
||||
double value;
|
||||
};
|
||||
|
||||
typedef struct vfc_probe_node vfc_probe_node;
|
||||
|
||||
|
||||
|
||||
/*
|
||||
* The probes structure. It simply acts as a wrapper for a Verificarlo hashmap.
|
||||
*/
|
||||
* The probes structure. It simply acts as a wrapper for a Verificarlo hashmap.
|
||||
*/
|
||||
|
||||
struct vfc_probes {
|
||||
vfc_hashmap_t map;
|
||||
@ -39,10 +35,9 @@ struct vfc_probes {
|
||||
|
||||
typedef struct vfc_probes vfc_probes;
|
||||
|
||||
|
||||
/*
|
||||
* Initialize an empty vfc_probes instance
|
||||
*/
|
||||
* Initialize an empty vfc_probes instance
|
||||
*/
|
||||
|
||||
vfc_probes vfc_init_probes() {
|
||||
vfc_probes probes;
|
||||
@ -51,20 +46,18 @@ vfc_probes vfc_init_probes() {
|
||||
return probes;
|
||||
}
|
||||
|
||||
|
||||
|
||||
/*
|
||||
* Free all probes
|
||||
*/
|
||||
* Free all probes
|
||||
*/
|
||||
|
||||
void vfc_free_probes(vfc_probes * probes) {
|
||||
void vfc_free_probes(vfc_probes *probes) {
|
||||
|
||||
// Before freeing the map, iterate manually over all items to free the keys
|
||||
vfc_probe_node * probe = NULL;
|
||||
for(int i = 0; i < probes->map->capacity; i++) {
|
||||
probe = (vfc_probe_node*) get_value_at(probes->map->items, i);
|
||||
if(probe != NULL) {
|
||||
if(probe->key != NULL) {
|
||||
vfc_probe_node *probe = NULL;
|
||||
for (int i = 0; i < probes->map->capacity; i++) {
|
||||
probe = (vfc_probe_node *)get_value_at(probes->map->items, i);
|
||||
if (probe != NULL) {
|
||||
if (probe->key != NULL) {
|
||||
free(probe->key);
|
||||
}
|
||||
}
|
||||
@ -73,14 +66,12 @@ void vfc_free_probes(vfc_probes * probes) {
|
||||
vfc_hashmap_free(probes->map);
|
||||
}
|
||||
|
||||
|
||||
|
||||
/*
|
||||
* Helper function to generate the key from test and variable name
|
||||
*/
|
||||
* Helper function to generate the key from test and variable name
|
||||
*/
|
||||
|
||||
char * gen_probe_key(char * testName, char * varName) {
|
||||
char * key = (char *) malloc(strlen(testName) + strlen(varName) + 2);
|
||||
char *gen_probe_key(char *testName, char *varName) {
|
||||
char *key = (char *)malloc(strlen(testName) + strlen(varName) + 2);
|
||||
strcpy(key, testName);
|
||||
strcat(key, ",");
|
||||
strcat(key, varName);
|
||||
@ -88,42 +79,33 @@ char * gen_probe_key(char * testName, char * varName) {
|
||||
return key;
|
||||
}
|
||||
|
||||
|
||||
|
||||
/*
|
||||
* Helper function to detect forbidden character ',' in the keys
|
||||
*/
|
||||
* Helper function to detect forbidden character ',' in the keys
|
||||
*/
|
||||
|
||||
void validate_probe_key(char * str) {
|
||||
void validate_probe_key(char *str) {
|
||||
unsigned int len = strlen(str);
|
||||
|
||||
for(unsigned int i=0; i<len; i++) {
|
||||
if(str[i] == ',') {
|
||||
fprintf(
|
||||
stderr,
|
||||
for (unsigned int i = 0; i < len; i++) {
|
||||
if (str[i] == ',') {
|
||||
fprintf(stderr,
|
||||
"Error [verificarlo]: One of your probes has a ',' in its test \
|
||||
or variable name (\"%s\"), which is forbidden\n",
|
||||
str
|
||||
);
|
||||
str);
|
||||
exit(1);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
/*
|
||||
* Add a new probe. If an issue with the key is detected (forbidden characters or
|
||||
* a duplicate key), an error will be thrown.
|
||||
*/
|
||||
* Add a new probe. If an issue with the key is detected (forbidden characters
|
||||
* or a duplicate key), an error will be thrown.
|
||||
*/
|
||||
|
||||
int vfc_put_probe(
|
||||
vfc_probes * probes,
|
||||
char * testName, char * varName,
|
||||
double val
|
||||
) {
|
||||
int vfc_put_probe(vfc_probes *probes, char *testName, char *varName,
|
||||
double val) {
|
||||
|
||||
if(probes == NULL) {
|
||||
if (probes == NULL) {
|
||||
return 1;
|
||||
}
|
||||
|
||||
@ -133,100 +115,85 @@ int vfc_put_probe(
|
||||
validate_probe_key(varName);
|
||||
|
||||
// Get the key, which is : testName + "," + varName
|
||||
char * key = gen_probe_key(testName, varName);
|
||||
char *key = gen_probe_key(testName, varName);
|
||||
|
||||
// Look for a duplicate key
|
||||
vfc_probe_node * oldProbe = (vfc_probe_node*) vfc_hashmap_get(
|
||||
probes->map, vfc_hashmap_str_function(key)
|
||||
);
|
||||
vfc_probe_node *oldProbe = (vfc_probe_node *)vfc_hashmap_get(
|
||||
probes->map, vfc_hashmap_str_function(key));
|
||||
|
||||
if(oldProbe != NULL) {
|
||||
if(strcmp(key, oldProbe->key) == 0) {
|
||||
fprintf(
|
||||
stderr,
|
||||
if (oldProbe != NULL) {
|
||||
if (strcmp(key, oldProbe->key) == 0) {
|
||||
fprintf(stderr,
|
||||
"Error [verificarlo]: you have a duplicate error with one of \
|
||||
your probes (\"%s\"). Please make sure to use different names.\n",
|
||||
key
|
||||
);
|
||||
key);
|
||||
exit(1);
|
||||
}
|
||||
}
|
||||
|
||||
// Insert the element in the hashmap
|
||||
vfc_probe_node * newProbe = (vfc_probe_node*) malloc(sizeof(vfc_probe_node));
|
||||
vfc_probe_node *newProbe = (vfc_probe_node *)malloc(sizeof(vfc_probe_node));
|
||||
newProbe->key = key;
|
||||
newProbe->value = val;
|
||||
|
||||
vfc_hashmap_insert(
|
||||
probes->map, vfc_hashmap_str_function(key), newProbe
|
||||
);
|
||||
vfc_hashmap_insert(probes->map, vfc_hashmap_str_function(key), newProbe);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
|
||||
/*
|
||||
* Remove (free) an element from the hash table
|
||||
*/
|
||||
* Remove (free) an element from the hash table
|
||||
*/
|
||||
|
||||
int vfc_remove_probe(vfc_probes * probes, char * testName, char * varName) {
|
||||
int vfc_remove_probe(vfc_probes *probes, char *testName, char *varName) {
|
||||
|
||||
if(probes == NULL) {
|
||||
if (probes == NULL) {
|
||||
return 1;
|
||||
}
|
||||
|
||||
// Get the key, which is : testName + "," + varName
|
||||
char * key = gen_probe_key(testName, varName);
|
||||
char *key = gen_probe_key(testName, varName);
|
||||
|
||||
vfc_hashmap_remove(probes->map, vfc_hashmap_str_function(key));
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
|
||||
/*
|
||||
* Return the number of probes stored in the hashmap
|
||||
*/
|
||||
* Return the number of probes stored in the hashmap
|
||||
*/
|
||||
|
||||
unsigned int vfc_num_probes(vfc_probes * probes) {
|
||||
unsigned int vfc_num_probes(vfc_probes *probes) {
|
||||
return vfc_hashmap_num_items(probes->map);
|
||||
}
|
||||
|
||||
|
||||
|
||||
/*
|
||||
* Dump probes in a .csv file (the double values are converted to hex), then
|
||||
* free it.
|
||||
*/
|
||||
* Dump probes in a .csv file (the double values are converted to hex), then
|
||||
* free it.
|
||||
*/
|
||||
|
||||
int vfc_dump_probes(vfc_probes * probes) {
|
||||
int vfc_dump_probes(vfc_probes *probes) {
|
||||
|
||||
if(probes == NULL) {
|
||||
if (probes == NULL) {
|
||||
return 1;
|
||||
}
|
||||
|
||||
// Get export path from the VFC_PROBES_OUTPUT env variable
|
||||
char* exportPath = getenv("VFC_PROBES_OUTPUT");
|
||||
if(!exportPath) {
|
||||
printf(
|
||||
"Warning [verificarlo]: VFC_PROBES_OUTPUT is not set, probes will \
|
||||
not be dumped\n"
|
||||
);
|
||||
char *exportPath = getenv("VFC_PROBES_OUTPUT");
|
||||
if (!exportPath) {
|
||||
printf("Warning [verificarlo]: VFC_PROBES_OUTPUT is not set, probes will \
|
||||
not be dumped\n");
|
||||
vfc_free_probes(probes);
|
||||
return 0;
|
||||
}
|
||||
|
||||
FILE * fp = fopen(exportPath, "w");
|
||||
FILE *fp = fopen(exportPath, "w");
|
||||
|
||||
if(fp == NULL) {
|
||||
fprintf(
|
||||
stderr,
|
||||
if (fp == NULL) {
|
||||
fprintf(stderr,
|
||||
"Error [verificarlo]: impossible to open the CSV file to save your \
|
||||
probes (\"%s\")\n",
|
||||
exportPath
|
||||
);
|
||||
exportPath);
|
||||
exit(1);
|
||||
}
|
||||
|
||||
@ -234,15 +201,11 @@ int vfc_dump_probes(vfc_probes * probes) {
|
||||
fprintf(fp, "test,variable,value\n");
|
||||
|
||||
// Iterate over all table elements
|
||||
vfc_probe_node * probe = NULL;
|
||||
for(int i = 0; i < probes->map->capacity; i++) {
|
||||
probe = (vfc_probe_node*) get_value_at(probes->map->items, i);
|
||||
if(probe != NULL) {
|
||||
fprintf(
|
||||
fp, "%s,%a\n",
|
||||
probe->key,
|
||||
probe->value
|
||||
);
|
||||
vfc_probe_node *probe = NULL;
|
||||
for (int i = 0; i < probes->map->capacity; i++) {
|
||||
probe = (vfc_probe_node *)get_value_at(probes->map->items, i);
|
||||
if (probe != NULL) {
|
||||
fprintf(fp, "%s,%a\n", probe->key, probe->value);
|
||||
}
|
||||
}
|
||||
|
||||
|
@ -25,10 +25,10 @@ void selectLargestDenominator(unsigned int l, unsigned int N_updates,
|
||||
for (unsigned int j = lbar; j < N_updates + 1; j++) {
|
||||
index = p[j];
|
||||
component = Updates_index[index - 1];
|
||||
breakdown = abs(1 + ylk[l][index][component]);
|
||||
breakdown = std::fabs(1 + ylk[l][index][component]);
|
||||
#ifdef DEBUG
|
||||
std::cout << "Inside selectLargestDenominator()" << std::endl;
|
||||
std::cout << "breakdown = abs(1 + ylk[" << l << "][" << index << "]["
|
||||
std::cout << "breakdown = fabs(1 + ylk[" << l << "][" << index << "]["
|
||||
<< component << "]) = " << breakdown << std::endl;
|
||||
std::cout << std::endl;
|
||||
#endif
|
||||
|
@ -76,7 +76,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
<< "] = " << beta << std::endl;
|
||||
std::cerr << std::endl;
|
||||
#endif
|
||||
if (fabs(beta) < threshold()) {
|
||||
if (std::fabs(beta) < threshold()) {
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
}
|
||||
@ -184,7 +184,7 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
}
|
||||
}
|
||||
#ifdef DEBUG
|
||||
showVector(ylk[0][k], Dim+1, "");
|
||||
showVector(ylk[0][k], Dim + 1, "");
|
||||
#endif
|
||||
}
|
||||
|
||||
@ -208,7 +208,7 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
<< "] = " << beta << std::endl;
|
||||
std::cerr << std::endl;
|
||||
#endif
|
||||
if (fabs(beta) < threshold()) {
|
||||
if (std::fabs(beta) < threshold()) {
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
for (unsigned int i = 1; i < Dim + 1; i++) {
|
||||
@ -237,7 +237,10 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
(i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] * ibeta;
|
||||
}
|
||||
}
|
||||
matMul(Al, last, next, Dim); tmp = next; next = last; last = tmp;
|
||||
matMul(Al, last, next, Dim);
|
||||
tmp = next;
|
||||
next = last;
|
||||
last = tmp;
|
||||
#ifdef DEBUG
|
||||
showMatrix(last, Dim, "last");
|
||||
#endif
|
||||
@ -275,17 +278,17 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
}
|
||||
|
||||
extern "C" {
|
||||
void MaponiA3_f(double **linSlater_inv, unsigned int *Dim,
|
||||
void MaponiA3_f(double **linSlater_inv, unsigned int *Dim,
|
||||
unsigned int *N_updates, double **linUpdates,
|
||||
unsigned int **Updates_index) {
|
||||
MaponiA3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
extern "C" {
|
||||
void MaponiA3S_f(double **linSlater_inv, unsigned int *Dim,
|
||||
void MaponiA3S_f(double **linSlater_inv, unsigned int *Dim,
|
||||
unsigned int *N_updates, double **linUpdates,
|
||||
unsigned int **Updates_index) {
|
||||
MaponiA3S(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
@ -23,7 +23,7 @@ void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
// Denominator
|
||||
double den = 1 + C[Updates_index[l] - 1];
|
||||
if (fabs(den) < threshold()) {
|
||||
if (std::fabs(den) < threshold()) {
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
}
|
||||
@ -71,7 +71,7 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
// Denominator
|
||||
double den = 1 + C[Updates_index[l] - 1];
|
||||
if (fabs(den) < threshold()) {
|
||||
if (std::fabs(den) < threshold()) {
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
|
||||
@ -131,7 +131,7 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
// Denominator
|
||||
double den = 1 + C[Updates_index[l] - 1];
|
||||
if (fabs(den) < threshold()) {
|
||||
if (std::fabs(den) < threshold()) {
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
|
||||
@ -173,7 +173,8 @@ 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)
|
||||
// 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 " << N_updates << " updates" << std::endl;
|
||||
@ -197,7 +198,7 @@ void SM4(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
|
||||
// Denominator
|
||||
double den = 1 + C[Updates_index[l] - 1];
|
||||
if (fabs(den) < threshold()) {
|
||||
if (std::fabs(den) < threshold()) {
|
||||
std::cerr << "Breakdown condition triggered at " << Updates_index[l]
|
||||
<< std::endl;
|
||||
|
||||
@ -237,27 +238,23 @@ void SM4(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
|
||||
}
|
||||
|
||||
extern "C" {
|
||||
void SM1_f(double **linSlater_inv, unsigned int *Dim,
|
||||
unsigned int *N_updates, double **linUpdates,
|
||||
unsigned int **Updates_index) {
|
||||
void SM1_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
|
||||
double **linUpdates, unsigned int **Updates_index) {
|
||||
SM1(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
||||
}
|
||||
|
||||
void SM2_f(double **linSlater_inv, unsigned int *Dim,
|
||||
unsigned int *N_updates, double **linUpdates,
|
||||
unsigned int **Updates_index) {
|
||||
SM2(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
||||
}
|
||||
|
||||
void SM3_f(double **linSlater_inv, unsigned int *Dim,
|
||||
unsigned int *N_updates, double **linUpdates,
|
||||
unsigned int **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);
|
||||
}
|
||||
}
|
||||
|
||||
void SM2_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
|
||||
double **linUpdates, unsigned int **Updates_index) {
|
||||
SM2(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
|
||||
}
|
||||
|
||||
void SM3_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
|
||||
double **linUpdates, unsigned int **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);
|
||||
}
|
||||
}
|
||||
|
@ -1,6 +1,6 @@
|
||||
// main.cpp
|
||||
#include "SM_Maponi.hpp"
|
||||
#include "SM_Helpers.hpp"
|
||||
#include "SM_Maponi.hpp"
|
||||
|
||||
int main() {
|
||||
|
||||
@ -9,24 +9,31 @@ int main() {
|
||||
|
||||
// Declare, allocate all vectors and matrices and fill them with zeros
|
||||
unsigned int *Ar_index = new unsigned int[M];
|
||||
double *A = new double[M*M]; // The matrix to be inverted
|
||||
double *A0 = new double[M*M]; // A diagonal matrix with the digonal elements of A
|
||||
double *Ar = new double[M*M]; // The update matrix
|
||||
double *A0_inv = new double[M*M]; // The inverse
|
||||
double *A = new double[M * M]; // The matrix to be inverted
|
||||
double *A0 =
|
||||
new double[M * M]; // A diagonal matrix with the digonal elements of A
|
||||
double *Ar = new double[M * M]; // The update matrix
|
||||
double *A0_inv = new double[M * M]; // The inverse
|
||||
|
||||
// Fill with zeros
|
||||
for (i = 0; i < M; i++) {
|
||||
for (j = 0; j < M; j++) {
|
||||
A0[i*M+j] = 0;
|
||||
Ar[i*M+j] = 0;
|
||||
A0_inv[i*M+j] = 0;
|
||||
A0[i * M + j] = 0;
|
||||
Ar[i * M + j] = 0;
|
||||
A0_inv[i * M + j] = 0;
|
||||
}
|
||||
}
|
||||
|
||||
// Initialize A with M=3 and fill acc. to Eq. (17) from paper
|
||||
A[0] = 1; A[3] = 1; A[6] = -1;
|
||||
A[1] = 1; A[4] = 1; A[7] = 0;
|
||||
A[2] = -1; A[5] = 0; A[8] = -1;
|
||||
A[0] = 1;
|
||||
A[3] = 1;
|
||||
A[6] = -1;
|
||||
A[1] = 1;
|
||||
A[4] = 1;
|
||||
A[7] = 0;
|
||||
A[2] = -1;
|
||||
A[5] = 0;
|
||||
A[8] = -1;
|
||||
|
||||
showMatrix(A, M, "A");
|
||||
|
||||
@ -34,20 +41,21 @@ int main() {
|
||||
// the inverse of A0_inv of diagonal matrix A0_inv
|
||||
// and the update matrix Ar
|
||||
for (i = 0; i < M; i++) {
|
||||
A0[i*M + i] = A[i*M + i];
|
||||
A0_inv[i*M + i] = 1.0/A[i*M + i];
|
||||
Ar_index[i] = i+1; // ! First column needs to start with 1 !
|
||||
A0[i * M + i] = A[i * M + i];
|
||||
A0_inv[i * M + i] = 1.0 / A[i * M + i];
|
||||
Ar_index[i] = i + 1; // ! First column needs to start with 1 !
|
||||
for (j = 0; j < M; j++) {
|
||||
Ar[i*M + j] = A[i*M + j] - A0[i*M + j];
|
||||
Ar[i * M + j] = A[i * M + j] - A0[i * M + j];
|
||||
}
|
||||
}
|
||||
|
||||
// Define pointers dim and n_updates to use in Sherman-Morrison(...) function call
|
||||
// Define pointers dim and n_updates to use in Sherman-Morrison(...) function
|
||||
// call
|
||||
MaponiA3(A0_inv, M, M, Ar, Ar_index);
|
||||
showMatrix(A0_inv, M, "A0_inv");
|
||||
|
||||
// Deallocate all vectors and matrices
|
||||
delete [] A, A0, A0_inv, Ar, Ar_index;
|
||||
delete[] A, A0, A0_inv, Ar, Ar_index;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
@ -1,22 +1,22 @@
|
||||
#include "hdf5/serial/hdf5.h"
|
||||
#include "hdf5/serial/H5Cpp.h"
|
||||
#include "hdf5/serial/hdf5.h"
|
||||
|
||||
#include "SM_Helpers.hpp"
|
||||
#include "SM_Maponi.hpp"
|
||||
#include "SM_Standard.hpp"
|
||||
#include "SM_Helpers.hpp"
|
||||
|
||||
using namespace H5;
|
||||
// #define DEBUG
|
||||
|
||||
const H5std_string FILE_NAME( "dataset.hdf5" );
|
||||
const H5std_string FILE_NAME("dataset.hdf5");
|
||||
|
||||
void read_int(H5File file, std::string key, unsigned int * data) {
|
||||
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) {
|
||||
void read_double(H5File file, std::string key, double *data) {
|
||||
DataSet ds = file.openDataSet(key);
|
||||
ds.read(data, PredType::IEEE_F64LE);
|
||||
ds.close();
|
||||
@ -32,19 +32,19 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
|
||||
read_int(file, group + "/slater_matrix_dim", &dim);
|
||||
read_int(file, group + "/nupdates", &nupdates);
|
||||
|
||||
double * slater_matrix = new double[dim*dim];
|
||||
double *slater_matrix = new double[dim * dim];
|
||||
read_double(file, group + "/slater_matrix", slater_matrix);
|
||||
|
||||
double * slater_inverse = new double[dim*dim];
|
||||
double *slater_inverse = new double[dim * dim];
|
||||
read_double(file, group + "/slater_inverse", slater_inverse);
|
||||
|
||||
unsigned int * col_update_index = new unsigned int[nupdates];
|
||||
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];
|
||||
double *updates = new double[nupdates * dim];
|
||||
read_double(file, group + "/updates", updates);
|
||||
|
||||
double * u = new double[nupdates*dim];
|
||||
double *u = new double[nupdates * dim];
|
||||
|
||||
/* Test */
|
||||
#ifdef DEBUG
|
||||
@ -55,8 +55,9 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
|
||||
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];
|
||||
u[i + j * dim] =
|
||||
updates[i + j * dim] - slater_matrix[i * dim + (col - 1)];
|
||||
slater_matrix[i * dim + (col - 1)] = updates[i + j * dim];
|
||||
}
|
||||
}
|
||||
|
||||
@ -93,20 +94,33 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
|
||||
showMatrix(slater_inverse, dim, "NEW Inverse");
|
||||
#endif
|
||||
|
||||
double * res = new double[dim*dim] {0};
|
||||
double *res = new double[dim * dim]{0};
|
||||
matMul(slater_matrix, slater_inverse, res, dim);
|
||||
bool ok = is_identity(res, dim, tolerance);
|
||||
|
||||
double res_max = residual_max(res, dim);
|
||||
double res2 = residual_frobenius2(res, dim);
|
||||
std::cout << "Residual = " << version << " " << cycle << " " << res_max << " " << res2 << std::endl;
|
||||
// double det;
|
||||
// double **tmp = new double *[dim];
|
||||
// for (int i = 0; i < dim; i++) {
|
||||
// tmp[i] = new double[dim];
|
||||
// for (int j = 0; j < dim; j++) {
|
||||
// tmp[i][j] = res[i * dim + j];
|
||||
// }
|
||||
// }
|
||||
// det = determinant(tmp, dim);
|
||||
// delete[] tmp;
|
||||
// std::cout << "Residual = " << version << " " << cycle << " " << res_max <<
|
||||
// " "
|
||||
// << res2 << " " << det << std::endl;
|
||||
std::cout << "Residual = " << version << " " << cycle << " " << res_max << " "
|
||||
<< res2 << std::endl;
|
||||
|
||||
#ifdef DEBUG
|
||||
showMatrix(res, dim, "Result");
|
||||
#endif
|
||||
|
||||
delete [] res, updates, u, col_update_index,
|
||||
slater_matrix, slater_inverse;
|
||||
delete[] res, updates, u, col_update_index, slater_matrix, slater_inverse;
|
||||
|
||||
return ok;
|
||||
}
|
||||
@ -114,7 +128,9 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
|
||||
int main(int argc, char **argv) {
|
||||
if (argc != 5) {
|
||||
std::cerr << "Execute from within 'datasets/'" << std::endl;
|
||||
std::cerr << "usage: test_h5 <version> <start cycle> <stop cycle> <tolerance>" << std::endl;
|
||||
std::cerr
|
||||
<< "usage: test_h5 <version> <start cycle> <stop cycle> <tolerance>"
|
||||
<< std::endl;
|
||||
return 1;
|
||||
}
|
||||
std::string version(argv[1]);
|
||||
@ -124,15 +140,12 @@ int main(int argc, char **argv) {
|
||||
H5File file(FILE_NAME, H5F_ACC_RDONLY);
|
||||
|
||||
bool ok;
|
||||
for (int cycle = start_cycle; cycle < stop_cycle+1; cycle++) {
|
||||
for (int cycle = start_cycle; cycle < stop_cycle + 1; cycle++) {
|
||||
ok = test_cycle(file, cycle, version, tolerance);
|
||||
if (ok) {
|
||||
std::cerr << "ok -- cycle " << std::to_string(cycle)
|
||||
<< std::endl;
|
||||
}
|
||||
else {
|
||||
std::cerr << "failed -- cycle " << std::to_string(cycle)
|
||||
<< std::endl;
|
||||
std::cerr << "ok -- cycle " << std::to_string(cycle) << std::endl;
|
||||
} else {
|
||||
std::cerr << "failed -- cycle " << std::to_string(cycle) << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
|
@ -3,31 +3,30 @@
|
||||
// cycles in a CSV file, instead of accepting a start and an end cycle (which
|
||||
// makes it easier to select the exact cycles we are interested in with vfc_ci).
|
||||
|
||||
#include <hdf5/serial/hdf5.h>
|
||||
#include <hdf5/serial/H5Cpp.h>
|
||||
#include <hdf5/serial/hdf5.h>
|
||||
|
||||
#include <vector>
|
||||
#include <fstream>
|
||||
#include <sstream>
|
||||
#include <vector>
|
||||
|
||||
|
||||
#include "SM_Helpers.hpp"
|
||||
#include "SM_Maponi.hpp"
|
||||
#include "SM_Standard.hpp"
|
||||
#include "SM_Helpers.hpp"
|
||||
#include "vfc_probe.h"
|
||||
|
||||
using namespace H5;
|
||||
// #define DEBUG
|
||||
|
||||
const H5std_string FILE_NAME( "datasets/ci_dataset.hdf5" );
|
||||
const H5std_string FILE_NAME("datasets/ci_dataset.hdf5");
|
||||
|
||||
void read_int(H5File file, std::string key, unsigned int * data) {
|
||||
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) {
|
||||
void read_double(H5File file, std::string key, double *data) {
|
||||
DataSet ds = file.openDataSet(key);
|
||||
ds.read(data, PredType::IEEE_F64LE);
|
||||
ds.close();
|
||||
@ -42,14 +41,15 @@ std::vector<int> get_cycles_list(std::string path) {
|
||||
std::string cycle_str;
|
||||
std::vector<int> cycles_list = {};
|
||||
|
||||
while(string_stream >> cycle_str) {
|
||||
while (string_stream >> cycle_str) {
|
||||
cycles_list.push_back(std::stoi(cycle_str));
|
||||
}
|
||||
|
||||
return cycles_list;
|
||||
}
|
||||
|
||||
int test_cycle(H5File file, int cycle, std::string version, vfc_probes * probes) {
|
||||
int test_cycle(H5File file, int cycle, std::string version,
|
||||
vfc_probes *probes) {
|
||||
|
||||
/* Read the data */
|
||||
|
||||
@ -59,11 +59,12 @@ int test_cycle(H5File file, int cycle, std::string version, vfc_probes * probes)
|
||||
// being zero-padded. This is used when calling vfc_put_probe later on.
|
||||
std::string zero_padded_group = std::to_string(cycle);
|
||||
zero_padded_group = "cycle_" +
|
||||
std::string(5 - zero_padded_group.length(), '0') + zero_padded_group;
|
||||
std::string(5 - zero_padded_group.length(), '0') +
|
||||
zero_padded_group;
|
||||
|
||||
try{
|
||||
try {
|
||||
file.openGroup(group);
|
||||
} catch(H5::Exception& e){
|
||||
} catch (H5::Exception &e) {
|
||||
std::cerr << "group " << group << "not found" << std::endl;
|
||||
return 0;
|
||||
}
|
||||
@ -72,21 +73,20 @@ int test_cycle(H5File file, int cycle, std::string version, vfc_probes * probes)
|
||||
read_int(file, group + "/slater_matrix_dim", &dim);
|
||||
read_int(file, group + "/nupdates", &nupdates);
|
||||
|
||||
|
||||
double * slater_matrix = new double[dim*dim];
|
||||
double *slater_matrix = new double[dim * dim];
|
||||
read_double(file, group + "/slater_matrix", slater_matrix);
|
||||
|
||||
double * slater_inverse = new double[dim*dim];
|
||||
double *slater_inverse = new double[dim * dim];
|
||||
read_double(file, group + "/slater_inverse", slater_inverse);
|
||||
//slater_inverse = transpose(slater_inverse, dim);
|
||||
// slater_inverse = transpose(slater_inverse, dim);
|
||||
|
||||
unsigned int * col_update_index = new unsigned int[nupdates];
|
||||
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];
|
||||
double *updates = new double[nupdates * dim];
|
||||
read_double(file, group + "/updates", updates);
|
||||
|
||||
double * u = new double[nupdates*dim];
|
||||
double *u = new double[nupdates * dim];
|
||||
|
||||
/* Test */
|
||||
#ifdef DEBUG
|
||||
@ -100,8 +100,9 @@ int test_cycle(H5File file, int cycle, std::string version, vfc_probes * probes)
|
||||
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];
|
||||
u[i + j * dim] =
|
||||
updates[i + j * dim] - slater_matrix[i * dim + (col - 1)];
|
||||
slater_matrix[i * dim + (col - 1)] = updates[i + j * dim];
|
||||
}
|
||||
}
|
||||
|
||||
@ -128,7 +129,7 @@ int test_cycle(H5File file, int cycle, std::string version, vfc_probes * probes)
|
||||
showMatrix(slater_inverse, dim, "NEW Inverse");
|
||||
#endif
|
||||
|
||||
double * res = new double[dim*dim] {0};
|
||||
double *res = new double[dim * dim]{0};
|
||||
matMul(slater_matrix, slater_inverse, res, dim);
|
||||
bool ok = is_identity(res, dim, 1e-3);
|
||||
|
||||
@ -139,11 +140,11 @@ int test_cycle(H5File file, int cycle, std::string version, vfc_probes * probes)
|
||||
showMatrix(res, dim, "Result");
|
||||
#endif
|
||||
|
||||
vfc_put_probe(probes, &(zero_padded_group)[0], &("res_max_" + version)[0], res_max);
|
||||
vfc_put_probe(probes, &(zero_padded_group)[0], &("res_max_" + version)[0],
|
||||
res_max);
|
||||
vfc_put_probe(probes, &(zero_padded_group)[0], &("res2_" + version)[0], res2);
|
||||
|
||||
delete [] res, updates, u, col_update_index,
|
||||
slater_matrix, slater_inverse;
|
||||
delete[] res, updates, u, col_update_index, slater_matrix, slater_inverse;
|
||||
|
||||
return ok;
|
||||
}
|
||||
@ -165,12 +166,9 @@ int main(int argc, char **argv) {
|
||||
for (int i = 0; i < cycles_list.size(); i++) {
|
||||
ok = test_cycle(file, cycles_list[i], version, &probes);
|
||||
if (ok) {
|
||||
std::cout << "ok -- cycle " << std::to_string(i)
|
||||
<< std::endl;
|
||||
}
|
||||
else {
|
||||
std::cerr << "failed -- cycle " << std::to_string(i)
|
||||
<< std::endl;
|
||||
std::cout << "ok -- cycle " << std::to_string(i) << std::endl;
|
||||
} else {
|
||||
std::cerr << "failed -- cycle " << std::to_string(i) << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
|
15
tools/restyle.sh
Executable file
15
tools/restyle.sh
Executable file
@ -0,0 +1,15 @@
|
||||
#!/usr/bin/env bash
|
||||
|
||||
STYLE='--style=LLVM'
|
||||
FORMATER='clang-format -i'
|
||||
|
||||
if [[ -z $SMVARS ]]
|
||||
then
|
||||
echo '$SMVARS is not set. Please source '/path/to/Sherman-Morrison/smvars.sh''
|
||||
exit 1
|
||||
fi
|
||||
|
||||
for ext in c cpp h hpp
|
||||
do
|
||||
find $SMROOT -type f -iname "*.${ext}" -exec echo "$FORMATER $STYLE" {} \;
|
||||
done
|
Loading…
Reference in New Issue
Block a user