From 48c27eb47ea60162d25ab1537eb19a956b8dbf05 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois?= Date: Thu, 18 Feb 2021 19:43:20 +0100 Subject: [PATCH] * Modified Maponi algo 3 so that it can be used with a number of updates that is smaller than the size of the Slater-matrix * Removed the Slater-matrix as an argument, since it is not used in the algo. * Added a manual 4x4 example to debug MaponiA3 to work with a number of updates that is smaller than the size of the Slater-matrix * Added a new Octave script to quickly check if the computes inverse is correct. --- .gitignore | 3 +- Makefile | 2 +- Octave-tests/SMMaponiA3.m | 25 +++-- QMCChem_dataset_test.f90 | 40 ++++---- QMCChem_dataset_test.m | 12 +++ SM_MaponiA3.cpp | 110 ++++++++++++--------- SM_MaponiA3.hpp | 2 +- SM_MaponiA3_mod.f90 | 4 +- cMaponiA3_test.cpp | 2 +- fMaponiA3_test.f90 | 197 +++++++++++++++++++++++++++++++------- fMaponiA3_test.f90.org | 52 ++++++++++ 11 files changed, 341 insertions(+), 108 deletions(-) create mode 100644 QMCChem_dataset_test.m create mode 100644 fMaponiA3_test.f90.org diff --git a/.gitignore b/.gitignore index d36ee77..6d31edd 100644 --- a/.gitignore +++ b/.gitignore @@ -4,4 +4,5 @@ cMaponiA3_test fMaponiA3_test QMCChem_dataset_test - +Slater_* +Updates* diff --git a/Makefile b/Makefile index 4d2890d..5c8d229 100644 --- a/Makefile +++ b/Makefile @@ -38,7 +38,7 @@ clean: @rm -vf *.o *.mod distclean: clean - @rm -vf cMaponiA3_test fMaponiA3_test QMCChem_dataset_test + @rm -vf cMaponiA3_test fMaponiA3_test QMCChem_dataset_test Slater_* Updates.dat ## Linking the C++ example program cMaponiA3_test: $(cMaponiA3_testOBJ) diff --git a/Octave-tests/SMMaponiA3.m b/Octave-tests/SMMaponiA3.m index 5779bd9..27af486 100644 --- a/Octave-tests/SMMaponiA3.m +++ b/Octave-tests/SMMaponiA3.m @@ -2,12 +2,12 @@ ## p. 283, doi:10.1016/j.laa.2006.07.007 clc ## Clear the screen -## Define the matrix to be inverted. This is example 8 from the paper -## In the future this matrix needs to be read from the function call arguments -A=[1,1,-1; ... - 1,1,0; ... - -1,0,-1]; -A0=diag(diag(A)); ## The diagonal part of A +### Define the matrix to be inverted. This is example 8 from the paper +### In the future this matrix needs to be read from the function call arguments +#A=[1,1,-1; ... +# 1,1,0; ... +# -1,0,-1]; +#A0=diag(diag(A)); ## The diagonal part of A ### The modified example that gives all singular updates at some point #A=[1,1,1; ... @@ -17,10 +17,18 @@ A0=diag(diag(A)); ## The diagonal part of A ### A square uniform distributed random integer matrix with entries in [-1,1] #do -# A=randi([-1,1],3,3); +# A=randi([-1,1],4,4); # A0=diag(diag(A)); ## The diagonal part of A #until (det(A)!=0 && det(A0)!=0) ## We need both matrices to be simultaniously non-singular +## Define the matrix to be inverted. This is example 8 from the paper +## In the future this matrix needs to be read from the function call arguments +A=[1,0, 1,-1; ... + 0,-1, 1, -1; ... + -1,0,-1, 0; ... + 1,1, 1, 1]; +A0=diag(diag(A)); ## The diagonal part of A + ### A square uniform distributed random float matrix with entries in (0,1) #do # A=rand(5); @@ -35,6 +43,7 @@ Ainv=zeros(nCols,nCols); ylk=zeros(nCols,nCols,nCols); p=zeros(nCols,1); breakdown=zeros(nCols,1); +cutOff=10e-6 A,A0 printf("Determinant of A is: %d\n",det(A)) @@ -109,4 +118,4 @@ else printf("Still not found. Giving up!\n"); IdTest endif -endif \ No newline at end of file +endif diff --git a/QMCChem_dataset_test.f90 b/QMCChem_dataset_test.f90 index d53bd01..4bebd0d 100644 --- a/QMCChem_dataset_test.f90 +++ b/QMCChem_dataset_test.f90 @@ -20,7 +20,7 @@ program QMCChem_dataset_test !! Write current S and S_inv to file for check in Octave open(unit = 2000, file = "Slater_old.dat") - open(unit = 3000, file = "Slater_inv_old.dat") + open(unit = 3000, file = "Slater_old_inv.dat") do i=1,dim do j=1,dim write(2000,"(E23.15, 1X)", advance="no") S(i,j) @@ -29,38 +29,44 @@ program QMCChem_dataset_test write(2000,*) write(3000,*) end do - call flush(2000) - call flush(3000) close(2000) close(3000) - !! Send S, S_inv and Updates to MaponiA3 algo - call MaponiA3(S, S_inv, dim, n_updates, Updates, Updates_index) + !! Write Updates to file to check + open(unit = 2000, file = "Updates.dat") + do i=1,dim + do j=1,n_updates + write(2000,"(E23.15, 1X)", advance="no") Updates(i,j) + end do + write(2000,*) + end do + close(2000) - !! Update S itself + !! Update S do j=1,n_updates do i=1,dim col = Updates_index(j) - S(i,col) = S(i,col) + Updates(i,col) + S(i,col) = S(i,col) + Updates(i,j) end do end do + !! Send S_inv and Updates to MaponiA3 algo + call MaponiA3(S_inv, dim, n_updates, Updates, Updates_index) + !! Write new S and S_inv to file for check in Octave - open(unit = 2000, file = "Slater_new.dat") - open(unit = 3000, file = "Slater_inv_new.dat") + open(unit = 4000, file = "Slater.dat") + open(unit = 5000, file = "Slater_inv.dat") do i=1,dim do j=1,dim - write(2000,"(E23.15, 1X)", advance="no") S(i,j) - write(3000,"(E23.15, 1X)", advance="no") S_inv(i,j) + write(4000,"(E23.15, 1X)", advance="no") S(i,j) + write(5000,"(E23.15, 1X)", advance="no") S_inv(i,j) end do - write(2000,*) - write(3000,*) + write(4000,*) + write(5000,*) end do - call flush(2000) - call flush(3000) - close(2000) - close(3000) + close(4000) + close(5000) deallocate(S, S_inv, Updates, Updates_index) end program diff --git a/QMCChem_dataset_test.m b/QMCChem_dataset_test.m new file mode 100644 index 0000000..cc075f3 --- /dev/null +++ b/QMCChem_dataset_test.m @@ -0,0 +1,12 @@ +Sold = dlmread('Slater_old.dat'); +Sold_inv = dlmread('Slater_old_inv.dat'); +S = dlmread('Slater.dat'); +S_inv = dlmread('Slater_inv.dat'); +det(Sold*transpose(Sold_inv)) +trace(Sold*transpose(Sold_inv)) +norm(Sold*transpose(Sold_inv)) + +det(S*transpose(S_inv)) +trace(S*transpose(S_inv)) +norm(S*transpose(S_inv)) + diff --git a/SM_MaponiA3.cpp b/SM_MaponiA3.cpp index f24a831..95dd69d 100644 --- a/SM_MaponiA3.cpp +++ b/SM_MaponiA3.cpp @@ -4,54 +4,75 @@ #include "SM_MaponiA3.hpp" #include "Helpers.hpp" -void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M, - unsigned int N_updates, double *Updates, - unsigned int *Updates_index) { +void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, + double *Updates, unsigned int *Updates_index) { - unsigned int k, l, lbar, i, j, tmp = M; - unsigned int *p = new unsigned int[M + 1]; + unsigned int k, l, lbar, i, j, tmp, component; + unsigned int *p = new unsigned int[N_updates + 1]; double alpha, beta; - double *breakdown = new double[M + 1]; - double *Al = new double[M * M]; + double *breakdown = new double[N_updates + 1]; + double *Al = new double[Dim * Dim]; p[0] = 0; - for (i = 0; i < M; i++) { + for (i = 0; i < N_updates; i++) { p[i + 1] = i + 1; } // Declare auxiliary solution matrix ylk - double ***ylk = new double **[M]; - for (l = 0; l < M; l++) { - ylk[l] = new double *[M + 1]; - for (k = 0; k < M + 1; k++) { - ylk[l][k] = new double[M + 1] {0}; + double ***ylk = new double **[N_updates]; + for (l = 0; l < N_updates; l++) { + ylk[l] = new double *[N_updates + 1]; + for (k = 0; k < N_updates + 1; k++) { + ylk[l][k] = new double[Dim + 1] {0}; } } - - // Calculate all the y0k in M^2 multiplications instead of M^3 + + cout << endl; + // Calculate all the y0k in Dim^2 multiplications instead of Dim^3 for (k = 1; k < N_updates + 1; k++) { - for (i = 1; i < M + 1; i++) { + for (i = 1; i < Dim + 1; i++) { + cout << k << "," << i << endl; ylk[0][k][i] = - Slater_inv[(i - 1) * M + (i - 1)] * Updates[(i - 1) * M + (k - 1)]; + Slater_inv[(i - 1) * Dim + (i - 1)] * Updates[(i - 1) * Dim + (k - 1)]; } } + cout << "y0k = " << endl; + for (k = 1; k < N_updates + 1; k++) { + cout << "[ "; + for (i = 1; i < Dim + 1; i++) { + cout << ylk[0][k][i] << " "; + } + cout << " ]" << endl; + } + cout << endl; + + showVector(p, N_updates + 1, "p"); + showVector(Updates_index, N_updates, "Updates_index"); + // Calculate all the ylk from the y0k - for (l = 1; l < M; l++) { - for (j = l; j < M + 1; j++) { - breakdown[j] = abs(1 + ylk[l - 1][p[j]][p[j]]); + for (l = 1; l < N_updates; l++) { + for (j = l; j < N_updates + 1; j++) { + component = Updates_index[p[j] - 1]; + cout << component << endl; + breakdown[j] = abs(1 + ylk[l - 1][p[j]][component]); + } + lbar = getMaxIndex(breakdown, N_updates + 1); + // Reset breakdown back to 0 for next round to avoid case where its first element is always the largest + for (i = 0; i < N_updates + 1; i++) { + breakdown[i] = 0; } - lbar = getMaxIndex(breakdown, M + 1); tmp = p[l]; p[l] = p[lbar]; p[lbar] = tmp; - for (k = l + 1; k < M + 1; k++) { - beta = 1 + ylk[l - 1][p[l]][p[l]]; - if (beta == 0) { - cout << "Break-down condition occured. Exiting..." << endl; - exit; - } - for (i = 1; i < M + 1; i++) { - alpha = ylk[l - 1][p[k]][p[l]] / beta; + component = Updates_index[p[l] - 1]; + beta = 1 + ylk[l - 1][p[l]][component]; + if (beta == 0) { + cout << "Break-down occured. Exiting..." << endl; + exit; + } + for (k = l + 1; k < N_updates + 1; k++) { + alpha = ylk[l - 1][p[k]][component] / beta; + for (i = 1; i < Dim + 1; i++) { ylk[l][p[k]][i] = ylk[l - 1][p[k]][i] - alpha * ylk[l - 1][p[l]][i]; } } @@ -63,26 +84,27 @@ void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M, double *copy = Slater_inv; // Construct A-inverse from A0-inverse and the ylk - for (l = 0; l < M; l++) { + for (l = 0; l < N_updates; l++) { k = l + 1; - beta = 1 + ylk[l][p[k]][p[k]]; - for (i = 0; i < M; i++) { - for (j = 0; j < M; j++) { - Al[i * M + j] = (i == j) - (j == p[k]-1) * ylk[l][p[k]][i + 1] / beta; + component = Updates_index[p[k] - 1]; + beta = 1 + ylk[l][p[k]][component]; + for (i = 0; i < Dim; i++) { + for (j = 0; j < Dim; j++) { + Al[i * Dim + j] = (i == j) - (j == p[k]-1) * ylk[l][p[k]][i + 1] / beta; } } - Slater_inv = matMul(Al, Slater_inv, M); + Slater_inv = matMul(Al, Slater_inv, Dim); } // Assign the new values of 'Slater_inv' to the old values in 'copy[][]' - for (i = 0; i < M; i++) { - for (j = 0; j < M; j++) { - copy[i * M + j] = Slater_inv[i * M + j]; + for (i = 0; i < Dim; i++) { + for (j = 0; j < Dim; j++) { + copy[i * Dim + j] = Slater_inv[i * Dim + j]; } } - for (l = 0; l < M; l++) { - for (k = 0; k < M + 1; k++) { + for (l = 0; l < N_updates; l++) { + for (k = 0; k < N_updates + 1; k++) { delete[] ylk[l][k]; } delete[] ylk[l]; @@ -92,11 +114,11 @@ void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M, } extern "C" { - void MaponiA3_f(double **linSlater0, double **linSlater_inv, unsigned int *Dim, + void MaponiA3_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates, double **linUpdates, unsigned int **Updates_index) { - MaponiA3(*linSlater0, *linSlater_inv, *Dim, - *N_updates, *linUpdates, - *Updates_index); + MaponiA3(*linSlater_inv, *Dim, + *N_updates, *linUpdates, + *Updates_index); } } diff --git a/SM_MaponiA3.hpp b/SM_MaponiA3.hpp index 305a2a9..8d8867a 100644 --- a/SM_MaponiA3.hpp +++ b/SM_MaponiA3.hpp @@ -1,3 +1,3 @@ -void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M, +void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index); diff --git a/SM_MaponiA3_mod.f90 b/SM_MaponiA3_mod.f90 index d0d8a4c..e672372 100644 --- a/SM_MaponiA3_mod.f90 +++ b/SM_MaponiA3_mod.f90 @@ -1,10 +1,10 @@ module Sherman_Morrison interface - subroutine MaponiA3(Slater0, Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="MaponiA3_f") + subroutine MaponiA3(Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="MaponiA3_f") use, intrinsic :: iso_c_binding, only : c_int, c_double integer(c_int), intent(in) :: dim, n_updates integer(c_int), dimension(:), allocatable, intent(in) :: Updates_index - real(c_double), dimension(:,:), allocatable, intent(in) :: Slater0, Updates + real(c_double), dimension(:,:), allocatable, intent(in) :: Updates real(c_double), dimension(:,:), allocatable, intent(in out) :: Slater_inv end subroutine MaponiA3 end interface diff --git a/cMaponiA3_test.cpp b/cMaponiA3_test.cpp index 0f0645e..8e9ffb4 100644 --- a/cMaponiA3_test.cpp +++ b/cMaponiA3_test.cpp @@ -48,7 +48,7 @@ int main() { // Define pointers dim and n_updates to use in Sherman-Morrison(...) function call unsigned int dim = M; unsigned int n_updates = M; - MaponiA3(A0, A0_inv, dim, n_updates, Ar, Ar_index); + MaponiA3(A0_inv, dim, n_updates, Ar, Ar_index); showMatrix(A0_inv, M, "A0_inv"); // Deallocate all vectors and matrices diff --git a/fMaponiA3_test.f90 b/fMaponiA3_test.f90 index b8c519c..16c425c 100644 --- a/fMaponiA3_test.f90 +++ b/fMaponiA3_test.f90 @@ -3,50 +3,181 @@ program Interface_test use, intrinsic :: iso_c_binding, only : c_int, c_double implicit none - integer i, j !! Iterators - integer(c_int) :: dim, N_updates - integer(c_int), dimension(:), allocatable :: Ar_index - real(c_double), dimension(:,:), allocatable :: A, A0, Ar - real(c_double), dimension(:,:), allocatable :: A0_inv + integer i, j, col !! Iterators + integer(c_int) :: Dim, N_updates + integer(c_int), dimension(:), allocatable :: Updates_index + real(c_double), dimension(:,:), allocatable :: S, A, Updates + real(c_double), dimension(:,:), allocatable :: S_inv, A_inv - dim = 3 - N_updates = 3 - allocate(Ar_index(dim), A(dim,dim), A0(dim,dim), Ar(dim,dim), A0_inv(dim,dim)) + Dim = 4 + N_updates = 2 + allocate( S(Dim,Dim), S_inv(Dim,Dim), A(Dim,Dim), A_inv(Dim,Dim), & + Updates(Dim,N_updates), Updates_index(N_updates)) + + !! Initialize S, S_inv, A and A_inv + S(1,1) = 1.0d0 + S(1,2) = 0.0d0 + S(1,3) = 1.0d0 + S(1,4) = -1.0d0 + S(2,1) = 0.0d0 + S(2,2) = 1.0d0 + S(2,3) = 1.0d0 + S(2,4) = 0.0d0 + S(3,1) = -1.0d0 + S(3,2) = 0.0d0 + S(3,3) = -1.0d0 + S(3,4) = 0.0d0 + S(4,1) = 1.0d0 + S(4,2) = 1.0d0 + S(4,3) = 1.0d0 + S(4,4) = 1.0d0 + + S_inv(1,1) = 1.0d0 + S_inv(1,2) = -1.0d0 + S_inv(1,3) = 1.0d0 + S_inv(1,4) = 1.0d0 + S_inv(2,1) = 1.0d0 + S_inv(2,2) = 0.0d0 + S_inv(2,3) = 2.0d0 + S_inv(2,4) = 1.0d0 + S_inv(3,1) = -1.0d0 + S_inv(3,2) = 1.0d0 + S_inv(3,3) = -2.0d0 + S_inv(3,4) = -1.0d0 + S_inv(4,1) = -1.0d0 + S_inv(4,2) = 0.0d0 + S_inv(4,3) = -1.0d0 + S_inv(4,4) = 0.0d0 - !! Initialize A with M=3 and fill acc. to Eq. (17) from paper A(1,1) = 1.0d0 - A(1,2) = 1.0d0 - A(1,3) = -1.0d0 - A(2,1) = 1.0d0 - A(2,2) = 1.0d0 - A(2,3) = 0.0d0 + A(1,2) = 0.0d0 + A(1,3) = 1.0d0 + A(1,4) = -1.0d0 + A(2,1) = 0.0d0 + A(2,2) = -1.0d0 + A(2,3) = 1.0d0 + A(2,4) = -1.0d0 A(3,1) = -1.0d0 A(3,2) = 0.0d0 A(3,3) = -1.0d0 + A(3,4) = 0.0d0 + A(4,1) = 1.0d0 + A(4,2) = 1.0d0 + A(4,3) = 1.0d0 + A(4,4) = 1.0d0 - !! Prepare the diagonal matrix A0 and the update matrix Ar - do i=1,dim - Ar_index(i) = i - do j=1,dim - if (i == j) then - A0(i,j) = A(i,j) - A0_inv(i,j) = 1.0d0 / A0(i,j) - else - A0(i,j) = 0.0d0 - A0_inv(i,j) = 0.0d0 - end if - Ar(i,j) = A(i,j) - A0(i,j) - end do - end do + A_inv(1,1) = 0.0d0 + A_inv(1,2) = -1.0d0 + A_inv(1,3) = -2.0d0 + A_inv(1,4) = -1.0d0 + A_inv(2,1) = 1.0d0 + A_inv(2,2) = 0.0d0 + A_inv(2,3) = 2.0d0 + A_inv(2,4) = 1.0d0 + A_inv(3,1) = 0.0d0 + A_inv(3,2) = 1.0d0 + A_inv(3,3) = 1.0d0 + A_inv(3,4) = 1.0d0 + A_inv(4,1) = -1.0d0 + A_inv(4,2) = 0.0d0 + A_inv(4,3) = -1.0d0 + A_inv(4,4) = 0.0d0 - call MaponiA3(A0, A0_inv, dim, n_updates, Ar, Ar_index) + !! Prepare Updates and Updates_index + Updates(1,1) = 0 + Updates(1,2) = 0 + Updates(2,1) = -2 + Updates(2,2) = -1 + Updates(3,1) = 0 + Updates(3,2) = 0 + Updates(4,1) = 0 + Updates(4,2) = 0 - do i=1,dim - do j=1,dim - write(*,"(F3.0,3X)", advance="no") A0_inv(i,j) + Updates_index(1) = 2 + Updates_index(2) = 4 + + write(*,*) + write(*,*) "Old S = " + do i=1,Dim + do j=1,Dim + write(*,"(F3.0,3X)", advance="no") S(i,j) end do write(*,*) end do - deallocate(Ar_index, A, A0, Ar, A0_inv) + write(*,*) + write(*,*) "Old S_inv = " + do i=1,Dim + do j=1,Dim + write(*,"(F3.0,3X)", advance="no") S_inv(i,j) + end do + write(*,*) + end do + + write(*,*) + write(*,*) "Updates = " + do i=1,Dim + do j=1,N_updates + write(*,"(F3.0,3X)", advance="no") Updates(i,j) + end do + write(*,*) + end do + + write(*,*) + write(*,*) "Updates_index = " + do i=1,N_updates + write(*,"(I1,3X)", advance="no") Updates_index(i) + end do + write(*,*) + + !! Update S + do i=1,N_updates + do j=1,Dim + col = Updates_index(i) + S(j,col) = S(j,col) + Updates(j,i) + end do + end do + + !! Update S_inv + call MaponiA3(S_inv, Dim, N_updates, Updates, Updates_index) + + write(*,*) + write(*,*) + write(*,*) "New computed S = " + do i=1,Dim + do j=1,Dim + write(*,"(F3.0,3X)", advance="no") S(i,j) + end do + write(*,*) + end do + + write(*,*) + write(*,*) "New actual S = " + do i=1,Dim + do j=1,Dim + write(*,"(F3.0,3X)", advance="no") A(i,j) + end do + write(*,*) + end do + + write(*,*) + write(*,*) + write(*,*) "New computed S_inv = " + do i=1,Dim + do j=1,Dim + write(*,"(F3.0,3X)", advance="no") S_inv(i,j) + end do + write(*,*) + end do + + write(*,*) + write(*,*) "New actual S_inv = " + do i=1,Dim + do j=1,Dim + write(*,"(F3.0,3X)", advance="no") A_inv(i,j) + end do + write(*,*) + end do + + deallocate(Updates_index, A, A_inv, S, Updates, S_inv) end program diff --git a/fMaponiA3_test.f90.org b/fMaponiA3_test.f90.org new file mode 100644 index 0000000..f63d88a --- /dev/null +++ b/fMaponiA3_test.f90.org @@ -0,0 +1,52 @@ +program Interface_test + use Sherman_Morrison, only : MaponiA3 + use, intrinsic :: iso_c_binding, only : c_int, c_double + implicit none + + integer i, j !! Iterators + integer(c_int) :: Dim, N_updates + integer(c_int), dimension(:), allocatable :: Updates_index + real(c_double), dimension(:,:), allocatable :: A, S, Updates + real(c_double), dimension(:,:), allocatable :: S_inv + + Dim = 3 + N_updates = 3 + allocate(Updates_index(Dim), A(Dim,Dim), S(Dim,Dim), Updates(Dim,Dim), S_inv(Dim,Dim)) + + !! Initialize A with M=3 and fill acc. to Eq. (17) from paper + A(1,1) = 1.0d0 + A(1,2) = 1.0d0 + A(1,3) = -1.0d0 + A(2,1) = 1.0d0 + A(2,2) = 1.0d0 + A(2,3) = 0.0d0 + A(3,1) = -1.0d0 + A(3,2) = 0.0d0 + A(3,3) = -1.0d0 + + !! Prepare the diagonal matrix S and the update matrix Updates + do i=1,Dim + Updates_index(i) = i + do j=1,Dim + if (i == j) then + S(i,j) = A(i,j) + S_inv(i,j) = 1.0d0 / S(i,j) + else + S(i,j) = 0.0d0 + S_inv(i,j) = 0.0d0 + end if + Updates(i,j) = A(i,j) - S(i,j) + end do + end do + + call MaponiA3(S_inv, Dim, N_updates, Updates, Updates_index) + + do i=1,Dim + do j=1,Dim + write(*,"(F3.0,3X)", advance="no") S_inv(i,j) + end do + write(*,*) + end do + + deallocate(Updates_index, A, S, Updates, S_inv) +end program