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++.
