mirror of
https://github.com/TREX-CoE/Sherman-Morrison.git
synced 2024-07-23 03:07:46 +02:00
Problem fixed: template function definitions (implementations) need to be fully known in the header file (.hpp) and not in the implementation file (.cpp).
This commit is contained in:
François Coppens
parent
7d55d7db77
commit
68d27b6df6
1
.gitignore
vendored
1
.gitignore
vendored
@ -1,3 +1,4 @@
|
||||
*.o
|
||||
SM-MaponiA3
|
||||
Sherman-Morrison
|
||||
.vscode
|
||||
127
SM-MaponiA3.cpp
127
SM-MaponiA3.cpp
@ -1,128 +1,9 @@
|
||||
// SM-MaponiA3.cpp
|
||||
#include "SM-MaponiA3.hpp"
|
||||
|
||||
uint getMaxIndex(double *arr, uint size) {
|
||||
uint i;
|
||||
uint max = arr[0];
|
||||
uint maxi = 0;
|
||||
for (i = 1; i < size; i++) {
|
||||
if (arr[i] > max) {
|
||||
max = arr[i];
|
||||
maxi = i;
|
||||
}
|
||||
}
|
||||
return maxi;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
void showScalar(T scalar, string name) {
|
||||
cout << name << " = " << scalar << endl << endl;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
void showVector(T* vector, uint size, string name) {
|
||||
cout << name << " = " << endl;
|
||||
for (uint i = 0; i < size; i++) {
|
||||
cout << "[ " << vector[i] << " ]" << endl;
|
||||
}
|
||||
cout << endl;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
void showMatrix(T** matrix, uint size, string name) {
|
||||
cout << name << " = " << endl;
|
||||
for (uint i = 0; i < size; i++) {
|
||||
cout << "[ ";
|
||||
for (uint j = 0; j < size; j++) {
|
||||
cout << matrix[i][j] << " ";
|
||||
}
|
||||
cout << " ]" << endl;
|
||||
}
|
||||
cout << endl;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
void showMatrixT(T** matrix, uint size, string name) {
|
||||
cout << name << " = " << endl;
|
||||
for (uint i = 0; i < size; i++) {
|
||||
cout << "[ ";
|
||||
for (uint j = 0; j < size; j++) {
|
||||
cout << matrix[j][i] << " ";
|
||||
}
|
||||
cout << " ]" << endl;
|
||||
}
|
||||
cout << endl;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
T** matMul(T** A, T** B, uint size) {
|
||||
T** C = new T*[size];
|
||||
for (uint i = 0; i < size; i++) {
|
||||
C[i] = new T[size];
|
||||
}
|
||||
for (uint i = 0; i < size; i++) {
|
||||
for (uint j = 0; j < size; j++) {
|
||||
for (uint k = 0; k < size; k++) {
|
||||
C[i][j] += A[i][k] * B[k][j];
|
||||
}
|
||||
}
|
||||
}
|
||||
return C;
|
||||
}
|
||||
|
||||
template<typename T1, typename T2>
|
||||
T1** outProd(T1* vec1, T2* vec2, uint size) {
|
||||
T1** C = new T1*[size];
|
||||
for (uint i = 0; i < size; i++) {
|
||||
C[i] = new T1[size];
|
||||
}
|
||||
for (uint i = 0; i < size; i++) {
|
||||
for (uint j = 0; j < size; j++) {
|
||||
C[i][j] = vec1[i+1] * vec2[j];
|
||||
}
|
||||
}
|
||||
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;
|
||||
}
|
||||
|
||||
void Sherman_Morrison(int **Slater0, double **Slater_inv, uint *Dim, uint *N_updates, int **Updates, uint *Updates_index) {
|
||||
uint k, l, lbar, i, j, tmp, M = *Dim;
|
||||
uint *p = new uint[M+1];
|
||||
void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, unsigned int *N_updates, int **Updates, unsigned int *Updates_index) {
|
||||
unsigned int k, l, lbar, i, j, tmp, M = *Dim;
|
||||
unsigned int *p = new unsigned int[M+1];
|
||||
double *breakdown = new double[M+1];
|
||||
double alpha, beta;
|
||||
|
||||
@ -219,4 +100,4 @@ void Sherman_Morrison(int **Slater0, double **Slater_inv, uint *Dim, uint *N_upd
|
||||
}
|
||||
delete [] ylk[l];
|
||||
}
|
||||
}
|
||||
}
|
||||
130
SM-MaponiA3.hpp
130
SM-MaponiA3.hpp
@ -4,12 +4,124 @@
|
||||
#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);
|
||||
void Sherman_Morrison(int **Slater0, double **Slater_inv, unsigned int *Dim, unsigned int *N_updates, int **Updates, unsigned int *Updates_index);
|
||||
|
||||
template<typename T>
|
||||
unsigned int getMaxIndex(T *vector, unsigned int size) {
|
||||
unsigned int i;
|
||||
unsigned int max = vector[0];
|
||||
unsigned int maxi = 0;
|
||||
for (i = 1; i < size; i++) {
|
||||
if (vector[i] > max) {
|
||||
max = vector[i];
|
||||
maxi = i;
|
||||
}
|
||||
}
|
||||
return maxi;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
void showScalar(T scalar, string name) {
|
||||
cout << name << " = " << scalar << endl << endl;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
void showVector(T* vector, unsigned int size, string name) {
|
||||
cout << name << " = " << endl;
|
||||
for (unsigned int i = 0; i < size; i++) {
|
||||
cout << "[ " << vector[i] << " ]" << endl;
|
||||
}
|
||||
cout << endl;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
void showMatrix(T** matrix, unsigned int size, string name) {
|
||||
cout << name << " = " << endl;
|
||||
for (unsigned int i = 0; i < size; i++) {
|
||||
cout << "[ ";
|
||||
for (unsigned int j = 0; j < size; j++) {
|
||||
cout << matrix[i][j] << " ";
|
||||
}
|
||||
cout << " ]" << endl;
|
||||
}
|
||||
cout << endl;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
void showMatrixT(T** matrix, unsigned int size, string name) {
|
||||
cout << name << " = " << endl;
|
||||
for (unsigned int i = 0; i < size; i++) {
|
||||
cout << "[ ";
|
||||
for (unsigned int j = 0; j < size; j++) {
|
||||
cout << matrix[j][i] << " ";
|
||||
}
|
||||
cout << " ]" << endl;
|
||||
}
|
||||
cout << endl;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
T** matMul(T** A, T** B, unsigned int size) {
|
||||
T** C = new T*[size];
|
||||
for (unsigned int i = 0; i < size; i++) {
|
||||
C[i] = new T[size];
|
||||
}
|
||||
for (unsigned int i = 0; i < size; i++) {
|
||||
for (unsigned int j = 0; j < size; j++) {
|
||||
for (unsigned int k = 0; k < size; k++) {
|
||||
C[i][j] += A[i][k] * B[k][j];
|
||||
}
|
||||
}
|
||||
}
|
||||
return C;
|
||||
}
|
||||
|
||||
template<typename T1, typename T2>
|
||||
T1** outProd(T1* vec1, T2* vec2, unsigned int size) {
|
||||
T1** C = new T1*[size];
|
||||
for (unsigned int i = 0; i < size; i++) {
|
||||
C[i] = new T1[size];
|
||||
}
|
||||
for (unsigned int i = 0; i < size; i++) {
|
||||
for (unsigned int j = 0; j < size; j++) {
|
||||
C[i][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;
|
||||
}
|
||||
15
main.cpp
15
main.cpp
@ -1,20 +1,19 @@
|
||||
// Algorithm 3 from P. Maponi,
|
||||
// p. 283, doi:10.1016/j.laa.2006.07.007
|
||||
|
||||
#include "SM-MaponiA3.hpp"
|
||||
#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
|
||||
unsigned int randRange = 1; // to get random integers in range [-randRange, randRange]
|
||||
unsigned int M = 3; // Dimension of the Slater-matrix
|
||||
unsigned int i, j; // Indices for iterators
|
||||
|
||||
// Declare and allocate all vectors and matrices
|
||||
uint *Ar_index = new uint[M];
|
||||
unsigned int *Ar_index = new unsigned int[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
|
||||
@ -64,8 +63,8 @@ int main() {
|
||||
}
|
||||
}
|
||||
|
||||
uint *dim = new uint(M);
|
||||
uint *n_updates = new uint(M);
|
||||
unsigned int *dim = new unsigned int(M);
|
||||
unsigned int *n_updates = new unsigned int(M);
|
||||
Sherman_Morrison(A0, A0_inv, dim, n_updates, Ar, Ar_index);
|
||||
|
||||
showMatrix(A0_inv, M, "A0_inv");
|
||||
|
||||
Loading…
Reference in New Issue
Block a user