TREX/Sherman-Morrison
1
0
Fork 0
mirror of https://github.com/TREX-CoE/Sherman-Morrison.git synced 2024-11-04 05:03:59 +01:00
Code Issues Projects Releases Wiki Activity

Split the program in a header file and an implementation file and included them in main. Does not compile. Get the following error: ld: /tmp/icpcnmVxvw.o: in function void showMatrix<int>(int**, unsigned int, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >)'

Browse Source
This commit is contained in:
François Coppens 2021-02-02 17:06:32 +01:00
parent 9d47bf9bb3
commit 7d55d7db77
4 changed files with 111 additions and 117 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
Makefile
View File

@ -5,12 +5,12 @@ CXXFLAGS=-O0 -debug full -traceback
# ARCH=-xCORE-AVX2 # ARCH=-xCORE-AVX2
DEPS = SM-MaponiA3.cpp DEPS = SM-MaponiA3.cpp
OBJ = SM-MaponiA3.o OBJ = SM-MaponiA3.o main.o
%.o: %.c $(DEPS) %.o: %.cpp $(DEPS)
$(CXX) $(ARCH) -c -o $@ $< $(CFLAGS) $(CXX) $(ARCH) -c -o $@ $< $(CFLAGS)
SM-MaponiA3: $(OBJ) Sherman-Morrison: $(OBJ)
$(CXX) $(ARCH) -o $@ $^ $(CFLAGS) $(CXX) $(ARCH) -o $@ $^ $(CFLAGS)
clean: clean:

123
SM-MaponiA3.cpp
View File

@ -1,107 +1,5 @@
// Algorithm 3 from P. Maponi, // SM-MaponiA3.cpp
// p. 283, doi:10.1016/j.laa.2006.07.007 #include "SM-MaponiA3.hpp"
#include <iostream>
#include <string>
#include <cmath>
#include <cstdlib>
#include <ctime>
using namespace std;
uint getMaxIndex(double *arr, uint size);
template<typename T>void showScalar(T scalar, string name);
template<typename T>void showVector(T *vector, uint size, string name);
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>T **matMul(T **A, T **B, uint size);
template<typename T1, typename T2>T1 **outProd(T1 *vec1, T2 *vec2, uint size);
template<typename T>T matDet(T **A, int M);
void Sherman_Morrison(int **Slater0, double **Slater_inv, uint *Dim, uint *N_updates, int **Updates, uint *Updates_index);
int main() {
srand((unsigned) time(0));
uint randRange = 1; // to get random integers in range [-randRange, randRange]
uint M = 3; // Dimension of the Slater-matrix
uint i, j; // Indices for iterators
// Declare and allocate all vectors and matrices
uint *indices_of_updates = new uint[M];
int **A = new int*[M]; // The matrix to be inverted
int **A0 = new int*[M]; // A diagonal matrix with the digonal elements of A
int **Ar = new int*[M]; // The update matrix
double **A0inv = new double*[M]; // Inverse of A0
for (i = 0; i < M; i++) {
A[i] = new int[M];
A0[i] = new int[M];
Ar[i] = new int[M];
A0inv[i] = new double[M];
}
// Initialize all matrices with zeros
for (i = 0; i < M; i++) {
for (j = 0; j < M; j++) {
A0[i][j] = 0;
Ar[i][j] = 0;
A0inv[i][j] = 0;
}
}
// Initialize A with M=3 and Eq. (17) from paper
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[2][0] = -1; A[2][1] = 0; A[2][2] = -1;
// // Fill A with random numbers from [-randRange,randRange]
// // and check if A and A0 are invertable
// do {
// for (i = 0; i < M; i++) {
// for (j = 0; j < M; j++) {
// A[i][j] = rand()%(2*randRange+1)-randRange;
// }
// }
// for (i = 0; i < M; i++) {
// A0[i][i] = A[i][i];
// }
// } while (matDet(A, M) == 0 || matDet(A0, M) == 0);
showMatrix(A, M, "A");
// Init Ar: the update matrix
for (i = 0; i < M; i++) {
for (j = 0; j < M; j++) {
Ar[i][j] = A[i][j] - A0[i][j];
}
}
// Define A0inv
for (i = 0; i < M; i++) {
A0inv[i][i] = 1.0/A[i][i];
}
showMatrix(A0inv, M, "A0inv");
// Populate indices_of_updates
for (i = 0; i < M; i++) {
indices_of_updates[i] = i;
}
uint *dim = new uint(M);
Sherman_Morrison(A0, A0inv, dim, dim, Ar, indices_of_updates);
showMatrixT(A0inv, M, "A0inv");
// Deallocate all vectors and matrices
for (i = 0; i < M; i++) {
delete [] A[i];
delete [] A0[i];
delete [] A0inv[i];
delete [] Ar[i];
}
delete [] A, A0, A0inv, Ar;
return 0;
}
uint getMaxIndex(double *arr, uint size) { uint getMaxIndex(double *arr, uint size) {
uint i; uint i;
@ -223,8 +121,6 @@ T matDet(T** A, int M) {
} }
void Sherman_Morrison(int **Slater0, double **Slater_inv, uint *Dim, uint *N_updates, int **Updates, uint *Updates_index) { void Sherman_Morrison(int **Slater0, double **Slater_inv, uint *Dim, uint *N_updates, int **Updates, uint *Updates_index) {
cout << "We just entered Sherman-Morrison" << endl;
showMatrix(Slater_inv, *Dim, "Slater_inv");
uint k, l, lbar, i, j, tmp, M = *Dim; uint k, l, lbar, i, j, tmp, M = *Dim;
uint *p = new uint[M+1]; uint *p = new uint[M+1];
double *breakdown = new double[M+1]; double *breakdown = new double[M+1];
@ -292,13 +188,11 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, uint *Dim, uint *N_upd
} }
} }
cout << "Just before the final construction of the inverse..." << endl;
showMatrix(Slater_inv, M, "Slater_inv");
// Construct A-inverse from A0-inverse and the ylk // Construct A-inverse from A0-inverse and the ylk
double** U; double **U;
double** Al = new double*[M]; double **Al = new double*[M];
for (i = 0; i < M; i++) Al[i] = new double[M]; for (i = 0; i < M; i++) Al[i] = new double[M];
for (l = 0; l < M; l++) { for (l = 0; l < M; l++) {
k = l+1; k = l+1;
U = outProd(ylk[l][p[k]], Id[p[k]-1], M); U = outProd(ylk[l][p[k]], Id[p[k]-1], M);
@ -310,10 +204,11 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, uint *Dim, uint *N_upd
} }
Slater_inv = matMul(Al, Slater_inv, M); Slater_inv = matMul(Al, Slater_inv, M);
} }
cout << "Just after the final construction of the inverse..." << endl;
showMatrix(Slater_inv, M, "Slater_inv"); delete [] p, breakdown;
for (i = 0; i < M; i++) { for (i = 0; i < M; i++) {
delete [] Id[i];
delete [] U[i]; delete [] U[i];
delete [] Al[i]; delete [] Al[i];
} }
@ -324,4 +219,4 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, uint *Dim, uint *N_upd
} }
delete [] ylk[l]; delete [] ylk[l];
} }
} }

15
SM-MaponiA3.hpp Normal file
View File

@ -0,0 +1,15 @@
// SM-MaponiA3.hpp
#include <iostream>
#include <cmath>
#include <string>
using namespace std;
uint getMaxIndex(double *arr, uint size);
template<typename T>void showScalar(T scalar, string name);
template<typename T>void showVector(T *vector, uint size, string name);
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>T **matMul(T **A, T **B, uint size);
template<typename T1, typename T2>T1 **outProd(T1 *vec1, T2 *vec2, uint size);
template<typename T>T matDet(T **A, int M);
void Sherman_Morrison(int **Slater0, double **Slater_inv, uint *Dim, uint *N_updates, int **Updates, uint *Updates_index);

84
main.cpp Normal file
View File

@ -0,0 +1,84 @@
// Algorithm 3 from P. Maponi,
// p. 283, doi:10.1016/j.laa.2006.07.007
#include <cstdlib>
#include <ctime>
#include "SM-MaponiA3.hpp"
using namespace std;
int main() {
srand((unsigned) time(0));
uint randRange = 1; // to get random integers in range [-randRange, randRange]
uint M = 3; // Dimension of the Slater-matrix
uint i, j; // Indices for iterators
// Declare and allocate all vectors and matrices
uint *Ar_index = new uint[M];
int **A = new int*[M]; // The matrix to be inverted
int **A0 = new int*[M]; // A diagonal matrix with the digonal elements of A
int **Ar = new int*[M]; // The update matrix
double **A0_inv = new double*[M]; // Inverse of A0
for (i = 0; i < M; i++) {
A[i] = new int[M];
A0[i] = new int[M];
Ar[i] = new int[M];
A0_inv[i] = new double[M];
}
// Initialize all matrices with zeros
for (i = 0; i < M; i++) {
for (j = 0; j < M; j++) {
A0[i][j] = 0;
Ar[i][j] = 0;
A0_inv[i][j] = 0;
}
}
// Initialize A with M=3 and Eq. (17) from paper
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[2][0] = -1; A[2][1] = 0; A[2][2] = -1;
// // Fill A with random numbers from [-randRange,randRange]
// // and check if A and A0 are invertable
// do {
// for (i = 0; i < M; i++) {
// for (j = 0; j < M; j++) {
// A[i][j] = rand()%(2*randRange+1)-randRange;
// }
// }
// for (i = 0; i < M; i++) {
// A0[i][i] = A[i][i];
// }
// } while (matDet(A, M) == 0 || matDet(A0, M) == 0);
showMatrix(A, M, "A");
// Init the update matrix Ar, A0_inv and Ar_index
for (i = 0; i < M; i++) {
A0_inv[i][i] = 1.0/A[i][i];
Ar_index[i] = i;
for (j = 0; j < M; j++) {
Ar[i][j] = A[i][j] - A0[i][j];
}
}
uint *dim = new uint(M);
uint *n_updates = new uint(M);
Sherman_Morrison(A0, A0_inv, dim, n_updates, Ar, Ar_index);
showMatrix(A0_inv, M, "A0_inv");
// Deallocate all vectors and matrices
for (i = 0; i < M; i++) {
delete [] A[i];
delete [] A0[i];
delete [] A0_inv[i];
delete [] Ar[i];
}
delete [] A, A0, A0_inv, Ar, Ar_index;
delete dim, n_updates;
return 0;
}