mirror of
https://github.com/TREX-CoE/Sherman-Morrison.git
synced 2024-11-03 20:54:08 +01:00
8d63dd1701
S = [1,0,1,-1; 0,1,1,0; -1,0,-1,0; 1,1,1,1]
S_inv = [1,-1,1,1; 1,0,2,1; -1,1,-2,-1; -1,0,-1,0]
u1 = [0,-2,0,0]
u2 = [0,-1,0,0]
upd_idx = [2,4]
To go from Maponi's examples where the number of updates is always equal
to the the dimension of the matrix, and the decomposition is always
diagonal, to cases with a non-diagonal decomposition and a number of
updates unequal to its size, the following changed needed to be made:
* in the calculation of the {y0k} an extra inner for-loop needs to be
added to make it a full matrix-vector multiplication due to the fact
that A0 is not a diagonal matrix
* in some places the use of the update-order vector p needs
the be replaced with that of upd_idx to make sure the correct
component of the ylk is selected and the proper rank-1 matrices are
constructed
* when a matrix is passed from Fortran to C++ with 2D adressing, it is
passed in colum-major order. The passed matrix needs to be transposed
before passing to C++. Doing this inside the algorithm will break
compatibility with called from C/C++.
59 lines
1.7 KiB
C++
59 lines
1.7 KiB
C++
// main.cpp
|
|
#include "SM_MaponiA3.hpp"
|
|
#include "Helpers.hpp"
|
|
#include <cstdlib>
|
|
#include <ctime>
|
|
|
|
int main() {
|
|
|
|
srand((unsigned) time(0));
|
|
unsigned int M = 3; // Dimension of the Slater-matrix
|
|
unsigned int i, j; // Indices for iterators
|
|
|
|
// Declare, allocate all vectors and matrices and fill them with zeros
|
|
unsigned int *Ar_index = new unsigned int[M];
|
|
double *A = new double[M*M]; // The matrix to be inverted
|
|
double *A0 = new double[M*M]; // A diagonal matrix with the digonal elements of A
|
|
double *Ar = new double[M*M]; // The update matrix
|
|
double *A0_inv = new double[M*M]; // The inverse
|
|
|
|
// Fill with zeros
|
|
for (i = 0; i < M; i++) {
|
|
for (j = 0; j < M; j++) {
|
|
A0[i*M+j] = 0;
|
|
Ar[i*M+j] = 0;
|
|
A0_inv[i*M+j] = 0;
|
|
}
|
|
}
|
|
|
|
// Initialize A with M=3 and fill acc. to Eq. (17) from paper
|
|
A[0] = 1; A[3] = 1; A[6] = -1;
|
|
A[1] = 1; A[4] = 1; A[7] = 0;
|
|
A[2] = -1; A[5] = 0; A[8] = -1;
|
|
|
|
showMatrix(A, M, "A");
|
|
|
|
// Initialize the diagonal matrix A0,
|
|
// the inverse of A0_inv of diagonal matrix A0_inv
|
|
// and the update matrix Ar
|
|
for (i = 0; i < M; i++) {
|
|
A0[i*M+i] = A[i*M+i];
|
|
A0_inv[i*M+i] = 1.0/A[i*M+i];
|
|
Ar_index[i] = i;
|
|
for (j = 0; j < M; j++) {
|
|
Ar[i*M+j] = A[i*M+j] - A0[i*M+j];
|
|
}
|
|
}
|
|
|
|
// Define pointers dim and n_updates to use in Sherman-Morrison(...) function call
|
|
unsigned int dim = M;
|
|
unsigned int n_updates = M;
|
|
MaponiA3(A0_inv, dim, n_updates, Ar, Ar_index);
|
|
showMatrix(A0_inv, M, "A0_inv");
|
|
|
|
// Deallocate all vectors and matrices
|
|
delete [] A, A0, A0_inv, Ar, Ar_index;
|
|
|
|
return 0;
|
|
}
|