Sherman-Morrison/include/SM_Helpers.hpp

128 lines
3.4 KiB
C++

// SM_Helpers.hpp
// Some usefull helper functions to support the Maponi algorithm.
#include <iostream>
#include <string>
#include <cstring>
#include <cmath>
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;
}