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

Optimize out Identity matrix

Browse Source
This commit is contained in:
Pablo Oliveira 2021-02-09 14:44:54 +01:00
parent 6f34f485d3
commit b575cda8d6
1 changed files with 10 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
SM_MaponiA3.cpp
View File

@ -8,7 +8,6 @@ void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M, unsigned int
unsigned int k, l, lbar, i, j, tmp = M; unsigned int k, l, lbar, i, j, tmp = M;
unsigned int *p = new unsigned int[M+1]; unsigned int *p = new unsigned int[M+1];
double *Id = new double[M*M];
double alpha, beta; double alpha, beta;
double *breakdown = new double[M+1]; double *breakdown = new double[M+1];
double *Al = new double[M*M]; double *Al = new double[M*M];
@ -26,14 +25,6 @@ void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M, unsigned int
} }
} }
// Initialize identity matrix
for (i = 0; i < M; i++) {
for (j = 0; j < M; j++) {
if (i != j) Id[i*M+j] = 0;
else Id[i*M+j] = 1;
}
}
// Initialize ylk with zeros // Initialize ylk with zeros
for (l = 0; l < M; l++) { for (l = 0; l < M; l++) {
for (k = 0; k < M+1; k++) { for (k = 0; k < M+1; k++) {
@ -80,14 +71,21 @@ void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M, unsigned int
// pointer 'copy' that points to whereever 'Slater_inv' points to now. // pointer 'copy' that points to whereever 'Slater_inv' points to now.
double *copy = Slater_inv; double *copy = Slater_inv;
double *U = new double[M*M];
// Construct A-inverse from A0-inverse and the ylk // Construct A-inverse from A0-inverse and the ylk
for (l = 0; l < M; l++) { for (l = 0; l < M; l++) {
k = l+1; k = l+1;
double * U = outProd(ylk[l][p[k]], (Id + (p[k]-1)*M), M); for (unsigned int i = 0; i < M; i++) {
for (unsigned int j = 0; j < M; j++) {
U[i*M+j] = ylk[l][p[k]][i+1] * ((p[k]-1) == j);
}
}
beta = 1 + ylk[l][p[k]][p[k]]; beta = 1 + ylk[l][p[k]][p[k]];
for (i = 0; i < M; i++) { for (i = 0; i < M; i++) {
for (j = 0; j < M; j++) { for (j = 0; j < M; j++) {
Al[i*M+j] = Id[i*M+j] - U[i*M+j] / beta; Al[i*M+j] = (i == j) - U[i*M+j] / beta;
} }
} }
Slater_inv = matMul(Al, Slater_inv, M); Slater_inv = matMul(Al, Slater_inv, M);
@ -106,7 +104,7 @@ void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M, unsigned int
} }
delete [] ylk[l]; delete [] ylk[l];
} }
delete [] Id, Al; delete [] Al, U;
delete [] p, breakdown; delete [] p, breakdown;
} }