Sherman-Morrison/include/SM_Helpers.hpp
Francois Coppens b863d40c05 - Setting break-down threshold trough smvars.sh
- Fixed SM3 to escape infinite update-loop
2021-04-15 18:13:26 +02:00

130 lines
3.2 KiB
C++

// SM_Helpers.hpp
// Some usefull helper functions to support the Maponi algorithm.
#include <cmath>
#include <cstring>
#include <iostream>
#include <string>
// #define DEBUG
#ifndef THRESHOLD
#define THRESHOLD 1e-3
#endif
double threshold();
void Switch(unsigned int *p, unsigned int l, unsigned int lbar);
void selectLargestDenominator(unsigned int l, unsigned int N_updates,
unsigned int *Updates_index, unsigned int *p,
double ***ylk);
template <typename T> void showScalar(T scalar, std::string name) {
std::cout << name << " = " << scalar << std::endl << std::endl;
}
template <typename T>
void showVector(T *vector, unsigned int size, std::string name) {
std::cout << name << " = " << std::endl;
for (unsigned int i = 0; i < size; i++) {
std::cout << "[ " << vector[i] << " ]" << std::endl;
}
std::cout << std::endl;
}
template <typename T>
void showMatrix(T *matrix, unsigned int M, std::string name) {
std::cout.precision(17);
std::cout << name << " = [" << std::endl;
for (unsigned int i = 0; i < M; i++) {
std::cout << "[";
for (unsigned int j = 0; j < M; j++) {
if (matrix[i * M + j] >= 0) {
std::cout << " " << matrix[i * M + j] << ",";
} else {
std::cout << " " << matrix[i * M + j] << ",";
}
}
std::cout << " ]," << std::endl;
}
std::cout << "]" << std::endl;
std::cout << std::endl;
}
template <typename T> T *transpose(T *A, unsigned int M) {
T *B = new T[M * M];
for (unsigned int i = 0; i < M; i++) {
for (unsigned int j = 0; j < M; j++) {
B[i * M + j] = A[i + j * M];
}
}
return B;
}
template <typename T> void matMul(T *A, T *B, T *C, unsigned int M) {
memset(C, 0, M * M * sizeof(T));
for (unsigned int i = 0; i < M; i++) {
for (unsigned int j = 0; j < M; j++) {
for (unsigned int k = 0; k < M; k++) {
C[i * M + j] += A[i * M + k] * B[k * M + j];
}
}
}
}
template <typename T1, typename T2>
T1 *outProd(T1 *vec1, T2 *vec2, unsigned int M) {
T1 *C = new T1[M * M];
for (unsigned int i = 0; i < M; i++) {
for (unsigned int j = 0; j < M; j++) {
C[i * M + j] = vec1[i + 1] * vec2[j];
}
}
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;
}
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)
return false;
if (i != j && fabs(A[i * M + j]) > tolerance)
return false;
}
}
return true;
}