2021-04-14 15:32:54 +02:00
|
|
|
// SM_Helpers.hpp
|
2021-02-03 12:13:09 +01:00
|
|
|
// Some usefull helper functions to support the Maponi algorithm.
|
2021-04-15 14:38:42 +02:00
|
|
|
#include <cmath>
|
|
|
|
#include <cstring>
|
2021-02-03 12:13:09 +01:00
|
|
|
#include <iostream>
|
|
|
|
#include <string>
|
|
|
|
|
2021-04-15 18:13:26 +02:00
|
|
|
// #define DEBUG
|
|
|
|
|
|
|
|
#ifndef THRESHOLD
|
|
|
|
#define THRESHOLD 1e-3
|
|
|
|
#endif
|
|
|
|
double threshold();
|
|
|
|
|
2021-04-14 16:19:49 +02:00
|
|
|
void Switch(unsigned int *p, unsigned int l, unsigned int lbar);
|
|
|
|
|
2021-04-15 14:38:42 +02:00
|
|
|
void selectLargestDenominator(unsigned int l, unsigned int N_updates,
|
|
|
|
unsigned int *Updates_index, unsigned int *p,
|
|
|
|
double ***ylk);
|
2021-04-14 16:19:49 +02:00
|
|
|
|
2021-04-15 14:38:42 +02:00
|
|
|
template <typename T> void showScalar(T scalar, std::string name) {
|
|
|
|
std::cout << name << " = " << scalar << std::endl << std::endl;
|
2021-02-03 12:13:09 +01:00
|
|
|
}
|
|
|
|
|
2021-04-15 14:38:42 +02:00
|
|
|
template <typename T>
|
2021-04-14 15:32:54 +02:00
|
|
|
void showVector(T *vector, unsigned int size, std::string name) {
|
2021-04-15 14:38:42 +02:00
|
|
|
std::cout << name << " = " << std::endl;
|
|
|
|
for (unsigned int i = 0; i < size; i++) {
|
|
|
|
std::cout << "[ " << vector[i] << " ]" << std::endl;
|
|
|
|
}
|
|
|
|
std::cout << std::endl;
|
2021-02-03 12:13:09 +01:00
|
|
|
}
|
|
|
|
|
2021-04-15 14:38:42 +02:00
|
|
|
template <typename T>
|
2021-04-14 15:32:54 +02:00
|
|
|
void showMatrix(T *matrix, unsigned int M, std::string name) {
|
2021-04-15 14:38:42 +02:00
|
|
|
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] << ",";
|
|
|
|
}
|
2021-02-03 12:13:09 +01:00
|
|
|
}
|
2021-04-15 14:38:42 +02:00
|
|
|
std::cout << " ]," << std::endl;
|
|
|
|
}
|
|
|
|
std::cout << "]" << std::endl;
|
|
|
|
std::cout << std::endl;
|
2021-02-03 12:13:09 +01:00
|
|
|
}
|
|
|
|
|
2021-04-15 14:38:42 +02:00
|
|
|
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];
|
The algorithm now works for the following 4x4 example with 2 updates:
S = [1,0,1,-1; 0,1,1,0; -1,0,-1,0; 1,1,1,1]
S_inv = [1,-1,1,1; 1,0,2,1; -1,1,-2,-1; -1,0,-1,0]
u1 = [0,-2,0,0]
u2 = [0,-1,0,0]
upd_idx = [2,4]
To go from Maponi's examples where the number of updates is always equal
to the the dimension of the matrix, and the decomposition is always
diagonal, to cases with a non-diagonal decomposition and a number of
updates unequal to its size, the following changed needed to be made:
* in the calculation of the {y0k} an extra inner for-loop needs to be
added to make it a full matrix-vector multiplication due to the fact
that A0 is not a diagonal matrix
* in some places the use of the update-order vector p needs
the be replaced with that of upd_idx to make sure the correct
component of the ylk is selected and the proper rank-1 matrices are
constructed
* when a matrix is passed from Fortran to C++ with 2D adressing, it is
passed in colum-major order. The passed matrix needs to be transposed
before passing to C++. Doing this inside the algorithm will break
compatibility with called from C/C++.
2021-02-21 18:28:08 +01:00
|
|
|
}
|
2021-04-15 14:38:42 +02:00
|
|
|
}
|
|
|
|
return B;
|
The algorithm now works for the following 4x4 example with 2 updates:
S = [1,0,1,-1; 0,1,1,0; -1,0,-1,0; 1,1,1,1]
S_inv = [1,-1,1,1; 1,0,2,1; -1,1,-2,-1; -1,0,-1,0]
u1 = [0,-2,0,0]
u2 = [0,-1,0,0]
upd_idx = [2,4]
To go from Maponi's examples where the number of updates is always equal
to the the dimension of the matrix, and the decomposition is always
diagonal, to cases with a non-diagonal decomposition and a number of
updates unequal to its size, the following changed needed to be made:
* in the calculation of the {y0k} an extra inner for-loop needs to be
added to make it a full matrix-vector multiplication due to the fact
that A0 is not a diagonal matrix
* in some places the use of the update-order vector p needs
the be replaced with that of upd_idx to make sure the correct
component of the ylk is selected and the proper rank-1 matrices are
constructed
* when a matrix is passed from Fortran to C++ with 2D adressing, it is
passed in colum-major order. The passed matrix needs to be transposed
before passing to C++. Doing this inside the algorithm will break
compatibility with called from C/C++.
2021-02-21 18:28:08 +01:00
|
|
|
}
|
|
|
|
|
2021-04-15 14:38:42 +02:00
|
|
|
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];
|
|
|
|
}
|
2021-03-01 16:08:14 +01:00
|
|
|
}
|
2021-04-15 14:38:42 +02:00
|
|
|
}
|
2021-03-01 16:08:14 +01:00
|
|
|
}
|
|
|
|
|
2021-04-15 14:38:42 +02:00
|
|
|
template <typename T1, typename T2>
|
2021-02-09 13:40:52 +01:00
|
|
|
T1 *outProd(T1 *vec1, T2 *vec2, unsigned int M) {
|
2021-04-15 14:38:42 +02:00
|
|
|
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];
|
2021-02-03 12:13:09 +01:00
|
|
|
}
|
2021-04-15 14:38:42 +02:00
|
|
|
}
|
|
|
|
return C;
|
2021-02-03 12:13:09 +01:00
|
|
|
}
|
|
|
|
|
2021-04-15 14:38:42 +02:00
|
|
|
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++;
|
2021-02-03 12:13:09 +01:00
|
|
|
k = 0;
|
2021-04-15 14:38:42 +02:00
|
|
|
}
|
2021-02-03 12:13:09 +01:00
|
|
|
}
|
2021-04-15 14:38:42 +02:00
|
|
|
}
|
|
|
|
det = det + A[0][p] * pow(-1, p) * matDet(temp, M - 1);
|
2021-02-03 12:13:09 +01:00
|
|
|
}
|
2021-04-15 14:38:42 +02:00
|
|
|
return det;
|
|
|
|
}
|
|
|
|
delete[] temp;
|
2021-02-09 13:40:52 +01:00
|
|
|
}
|
2021-02-16 10:49:15 +01:00
|
|
|
|
2021-04-15 14:38:42 +02:00
|
|
|
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++) {
|
2021-05-12 16:10:27 +02:00
|
|
|
if (i == j && std::fabs(A[i * M + j] - 1) > tolerance)
|
2021-04-15 14:38:42 +02:00
|
|
|
return false;
|
2021-05-12 16:10:27 +02:00
|
|
|
if (i != j && std::fabs(A[i * M + j]) > tolerance)
|
2021-04-15 14:38:42 +02:00
|
|
|
return false;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
return true;
|
2021-02-16 10:49:15 +01:00
|
|
|
}
|
2021-05-07 17:11:04 +02:00
|
|
|
|
2021-05-10 14:50:45 +02:00
|
|
|
template <typename T> T norm_max(T *A, unsigned int Dim) {
|
2021-05-07 17:11:04 +02:00
|
|
|
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];
|
2021-05-12 16:10:27 +02:00
|
|
|
delta = std::fabs(delta);
|
2021-05-07 17:11:04 +02:00
|
|
|
if (delta > res) {
|
|
|
|
res = delta;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
return res;
|
|
|
|
}
|
|
|
|
|
2021-05-10 14:50:45 +02:00
|
|
|
template <typename T> T norm_frobenius2(T *A, unsigned int Dim) {
|
2021-05-07 17:11:04 +02:00
|
|
|
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];
|
2021-05-10 14:50:45 +02:00
|
|
|
res += delta * delta;
|
2021-05-07 17:11:04 +02:00
|
|
|
}
|
|
|
|
}
|
|
|
|
return res;
|
|
|
|
}
|
|
|
|
|
2021-05-10 14:50:45 +02:00
|
|
|
template <typename T> T residual_max(T *A, unsigned int Dim) {
|
2021-05-07 17:11:04 +02:00
|
|
|
T res = 0;
|
|
|
|
for (unsigned int i = 0; i < Dim; i++) {
|
|
|
|
for (unsigned int j = 0; j < Dim; j++) {
|
2021-05-10 14:50:45 +02:00
|
|
|
T delta = A[i * Dim + j] - (i == j);
|
2021-05-12 16:10:27 +02:00
|
|
|
delta = std::fabs(delta);
|
2021-05-07 17:11:04 +02:00
|
|
|
if (delta > res) {
|
|
|
|
res = delta;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
return res;
|
|
|
|
}
|
|
|
|
|
2021-05-10 14:50:45 +02:00
|
|
|
template <typename T> T residual_frobenius2(T *A, unsigned int Dim) {
|
2021-05-07 17:11:04 +02:00
|
|
|
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);
|
2021-05-10 14:50:45 +02:00
|
|
|
res += delta * delta;
|
2021-05-07 17:11:04 +02:00
|
|
|
}
|
|
|
|
}
|
|
|
|
return res;
|
|
|
|
}
|