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

Added random matrix initialisation of arbitrary size. Added function to calculate determinant to test matrix invertability of A and A0.

Browse Source
This commit is contained in:
François Coppens 2021-01-29 12:53:20 +01:00
parent 8f54ca1124
commit 7d3ffe1d2a
1 changed files with 56 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

66
SM-MaponiA3.cpp
View File

@ -3,6 +3,9 @@
#include <iostream> #include <iostream>
#include <string> #include <string>
#include <cmath>
#include <cstdlib>
#include <ctime>
using namespace std; using namespace std;
uint getMaxIndex(double* arr, uint size); uint getMaxIndex(double* arr, uint size);
@ -12,9 +15,12 @@ template<typename T>void showMatrix(T** matrix, uint size, string name);
template<typename T>void showMatrixT(T** matrix, uint size, string name); template<typename T>void showMatrixT(T** matrix, uint size, string name);
template<typename T>T** matMul(T** A, T** B, uint size); template<typename T>T** matMul(T** A, T** B, uint size);
template<typename T1, typename T2>T1** outProd(T1* vec1, T2* vec2, uint size); template<typename T1, typename T2>T1** outProd(T1* vec1, T2* vec2, uint size);
template<typename T>T matDet(T** A, int M);
int main() { int main() {
srand((unsigned) time(0));
uint randRange = 1; // to get random integers in range [-randRange, randRange]
uint M = 3; uint M = 3;
uint i, j, k, l, lbar, tmp; uint i, j, k, l, lbar, tmp;
double alpha, beta; double alpha, beta;
@ -49,8 +55,6 @@ int main() {
// Initialize all matrices with zeros // Initialize all matrices with zeros
for (i = 0; i < M; i++) { for (i = 0; i < M; i++) {
for (j = 0; j < M; j++) { for (j = 0; j < M; j++) {
A[i][j] = 0;
//Ainv[i][j] = 0;
A0[i][j] = 0; A0[i][j] = 0;
A0inv[i][j] = 0; A0inv[i][j] = 0;
Ar[i][j] = 0; Ar[i][j] = 0;
@ -67,10 +71,17 @@ int main() {
} }
} }
// Initialize A // // Initialize A with M=3 and Eq. (17) from paper
A[0][0] = 1; A[0][1] = 1; A[0][2] = -1; // A[0][0] = 1; A[0][1] = 1; A[0][2] = -1;
A[1][0] = 1; A[1][1] = 1; A[1][2] = 0; // A[1][0] = 1; A[1][1] = 1; A[1][2] = 0;
A[2][0] = -1; A[2][1] = 0; A[2][2] = -1; // A[2][0] = -1; A[2][1] = 0; A[2][2] = -1;
// Fill A with random numbers from [-1,1]
for (i = 0; i < M; i++) {
for (j = 0; j < M; j++) {
A[i][j] = rand()%(2*randRange+1)-randRange;
}
}
// Define identity matrix, A0, A0inv and p // Define identity matrix, A0, A0inv and p
p[0] = 0; p[0] = 0;
@ -100,7 +111,7 @@ int main() {
} }
// showVector(ylk[0][k], M+1, "y0k"); // showVector(ylk[0][k], M+1, "y0k");
} }
showMatrixT(ylk[0], M+1, "y0k"); // showMatrixT(ylk[0], M+1, "y0k");
// Calculate all the ylk from the y0k // Calculate all the ylk from the y0k
// showVector(p, M+1, "p"); // showVector(p, M+1, "p");
@ -131,9 +142,8 @@ int main() {
// showVector(ylk[l][p[k]], M+1, "ylk"); // showVector(ylk[l][p[k]], M+1, "ylk");
} }
} }
showMatrixT(ylk[1], M+1, "y1k"); // showMatrixT(ylk[1], M+1, "y1k");
showMatrixT(ylk[2], M+1, "y2k"); // showMatrixT(ylk[2], M+1, "y2k");
// EVERYTHING WORKS UPTO HERE
// Construct A-inverse from A0-inverse and the ylk // Construct A-inverse from A0-inverse and the ylk
double** U; double** U;
@ -258,3 +268,39 @@ T1** outProd(T1* vec1, T2* vec2, uint size) {
} }
return C; return C;
} }
template<typename T>
T matDet(T** A, 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;
}