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

Regularised main loop.

Browse Source
This commit is contained in:
François Coppens 2021-03-04 18:33:28 +01:00
parent 9c82092cff
commit 066e50914c
1 changed files with 12 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
src/SM_MaponiA3.cpp
View File

@ -62,32 +62,32 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
// showMatrix(Slater_inv, Dim, "Slater_inv"); // showMatrix(Slater_inv, Dim, "Slater_inv");
// Calculate all the ylk from the y0k // Calculate all the ylk from the y0k
for (l = 1; l < N_updates; l++) { for (l = 0; l < N_updates; l++) {
// Select update with largest break-down val // Select update with largest break-down val
lbar = l; max = 0; lbar = l+1; max = 0;
for (j = l; j < N_updates + 1; j++) { for (j = l+1; j < N_updates + 1; j++) {
component = Updates_index[p[j] - 1]; component = Updates_index[p[j] - 1];
breakdown = abs(1 + ylk[l - 1][p[j]][component]); breakdown = abs(1 + ylk[l+1 - 1][p[j]][component]);
if (breakdown > max) { if (breakdown > max) {
max = breakdown; max = breakdown;
lbar = j; lbar = j;
} }
} tmp = p[l]; p[l] = p[lbar]; p[lbar] = tmp; } tmp = p[l+1]; p[l+1] = p[lbar]; p[lbar] = tmp;
component = Updates_index[p[l] - 1]; component = Updates_index[p[l+1] - 1];
beta = 1 + ylk[l - 1][p[l]][component]; beta = 1 + ylk[l+1 - 1][p[l+1]][component];
if (fabs(beta) < 1e-6) { if (beta == 0) {
cout << "Break-down occured. Exiting..." << endl; cout << "Break-down occured. Exiting..." << endl;
exit(1); exit(1);
} }
for (k = l + 1; k < N_updates + 1; k++) { for (k = l+1 + 1; k < N_updates + 1; k++) {
alpha = ylk[l - 1][p[k]][component] / beta; alpha = ylk[l+1 - 1][p[k]][component] / beta;
cout << "( l, k, p[k], component ) = (" << l << ", " << k << ", " << p[k] << ", " << component << ")" << endl; cout << "( l, k, p[k], component ) = (" << l << ", " << k << ", " << p[k] << ", " << component << ")" << endl;
for (i = 1; i < Dim + 1; i++) { for (i = 1; i < Dim + 1; i++) {
cout << "ylk[" << l << "][p[" << k << "]][" << i << "] = ylk[" << l - 1 << "][p[" << k << "]][" << i << "] - alpha * ylk[" << l - 1 << "][p[" << l << "]][" << i << "]" << endl; cout << "ylk[" << l << "][p[" << k << "]][" << i << "] = ylk[" << l - 1 << "][p[" << k << "]][" << i << "] - alpha * ylk[" << l - 1 << "][p[" << l << "]][" << i << "]" << endl;
ylk[l][p[k]][i] = ylk[l - 1][p[k]][i] ylk[l+1][p[k]][i] = ylk[l+1 - 1][p[k]][i]
- alpha * ylk[l - 1][p[l]][i]; - alpha * ylk[l+1 - 1][p[l+1]][i];
} }
} }
} }