Maponi algorithm 4 works for the modified Example 18 in the paper but not for general non-singlual NxN matrices yet.

.vscode/launch.json

@@ -0,0 +1,28 @@
+{
+    // Use IntelliSense to learn about possible attributes.
+    // Hover to view descriptions of existing attributes.
+    // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387
+    "version": "0.2.0",
+    "configurations": [
+        {
+            "name": "(GDb) Launch",
+            "type": "cppdbg",
+            "request": "launch",
+            "program": "${workspaceFolder}/SM-MaponiA3",
+            "args": [],
+            "stopAtEntry": false,
+            "cwd": "${workspaceFolder}",
+            "environment": [],
+            "externalConsole": false,
+            "MIMode": "gdb",
+            "miDebuggerPath": "/opt/intel/compilers_and_libraries/linux/bin/intel64/gdb-ia"
+            "setupCommands": [
+                {
+                    "description": "Enable pretty-printing for gdb",
+                    "text": "-enable-pretty-printing",
+                    "ignoreFailures": true
+                }
+            ]
+        }
+    ]
+}

.vscode/tasks.json

@@ -0,0 +1,27 @@
+{
+    "tasks": [
+        {
+            "type": "cppbuild",
+            "label": "C/C++: clang++ build active file",
+            "command": "/opt/intel/compilers_and_libraries_2020.4.304/linux/bin/clang++",
+            "args": [
+                "-g",
+                "${file}",
+                "-o",
+                "${fileDirname}/${fileBasenameNoExtension}"
+            ],
+            "options": {
+                "cwd": "${workspaceFolder}"
+            },
+            "problemMatcher": [
+                "$gcc"
+            ],
+            "group": {
+                "kind": "build",
+                "isDefault": true
+            },
+            "detail": "Task generated by Debugger."
+        }
+    ],
+    "version": "2.0.0"
+}

Makefile

@@ -1,6 +1,7 @@
 CC=icc
 CXX=icpc
-CFLAGS=
+CFLAGS=-xCORE-AVX2 -O0 -debug full
+CXXFLAGS=-xCORE-AVX2 -O0 -debug full
 
 DEPS = SM-MaponiA3.cpp
 OBJ = SM-MaponiA3.o 
@@ -10,3 +11,6 @@ OBJ = SM-MaponiA3.o
 
 SM-MaponiA3: $(OBJ)
 	$(CXX) -o $@ $^ $(CFLAGS)
+
+clean:
+	@rm -vf *.o

SMMaponiA4.m

@@ -8,15 +8,15 @@ clc ## Clear the screen
 %    1,1,0; ...
 %   -1,0,-1];
 
-## The modified example that gives all singular updates at some point
+%## The modified example that gives all singular updates at some point
-A=[1,1,1; ...
+%A=[1,1,1; ...
-   1,1,0; ...
+%   1,1,0; ...
-  -1,0,-1];
+%  -1,0,-1];
 
-%## A square uniform distributed random integer matrix with entries in [-1,1]  
+## A square uniform distributed random integer matrix with entries in [-1,1]  
-%do
+do
-%  A=randi([-5,5],4,4);
+  A=randi([-5,5],5,5);
-%until (det(A)!=0) ## We neec matrix non-singular
+until (det(A)!=0) ## We need matrix non-singular
 
 % ## A square uniform distributed random float matrix with entries in (0,1)  
 % do
@@ -26,6 +26,8 @@ A=[1,1,1; ...
 
 nCols=columns(A); ## The number of coluns of A (M in accompanying PDF)
 Id=eye(nCols);
+% d=norm(A);
+% d=max(eig(A));
 d=0.1
 A0=d*eye(nCols);
 Ar=A-A0;           ## The remainder of A
@@ -34,6 +36,8 @@ Ainv=zeros(nCols,nCols);
 ylk=zeros(nCols,nCols,nCols);
 breakdown=zeros(nCols,1);
 cutOff=1e-10;
+p=zeros(nCols,1);
+P=zeros(nCols);
 
 ## Calculate the inverse of A0 and populate p-vector
 for i=1:nCols
@@ -59,24 +63,29 @@ endfor
 
 
 ## Calculate all the ylk from the y0k calculated previously
-for l=2:3#nCols
+for l=2:nCols
   el=Id(:,l-1);
   ## Calculate break-down conditions and put in a vector
   for j=l-1:nCols
-    #breakdown(j) = abs( el(j) + ylk(j,l-1,l-1) )
+    %breakdown(j) = abs( el(j) + ylk(j,l-1,l-1) ) ## Condition from Maponi. Probably not correct!
+    l
+    ylk(j,j,l-1);
     breakdown(j) = abs( el(j) + ylk(j,j,l-1) )
     % printf("|el(%d) + ylk(%d,%d,%d)|\n", j, j, l-1, l-1);
   endfor
-  [val, s] = max(breakdown); ## Find the index of the max value
+  [val, s] = max(breakdown) ## Find the index of the max value
-  printf("l = %d\ns = %d\n",l-1,s);
+  % printf("l = %d\ns = %d\n",l-1,s);
   breakdown=zeros(nCols,1); ## Reset the entries to zero for next l-round
   if (s!=l-1) ## Apply partial pivoting
   ## Swap yl-1,k(r) and yl-1,k(s) for all k=l,l+1,...,M
-    r=l-1
+    r=l-1;
+    tmp=p(r);
+    p(r)=p(s);
+    p(s)=tmp;
     for k=l-1:nCols
-      tmp=ylk(r,k,l-1)
+      tmp=ylk(r,k,l-1);
       ylk(r,k,l-1)=ylk(s,k,l-1);
-      printf("ylk(%d,%d,%d)=ylk(%d,%d,%d)\n",r,k,l-1,s,k,l-1);
+      % printf("ylk(%d,%d,%d)=ylk(%d,%d,%d)\n",r,k,l-1,s,k,l-1);
       ylk(s,k,l-1)=tmp;
     endfor
     ## Modify yl-1,r and yl-1,s
@@ -89,7 +98,7 @@ for l=2:3#nCols
   for k=l:nCols
     for i=1:nCols
       ylk(i,k,l) = ylk(i,k,l-1) - ylk(l-1,k,l-1) / (1 + ylk(l-1,l-1,l-1)) * ylk(i,l-1,l-1);
-      printf("ylk(%d,%d,%d) = ylk(%d,%d,%d) - (ylk(%d,%d,%d) / (1+ylk(%d,%d,%d) * (ylk(%d,%d,%d)\n",i,k,l,i,k,l-1,l-1,k,l-1,l-1,l-1,l-1,i,l-1,l-1);
+      % printf("ylk(%d,%d,%d) = ylk(%d,%d,%d) - (ylk(%d,%d,%d) / (1+ylk(%d,%d,%d) * (ylk(%d,%d,%d)\n",i,k,l,i,k,l-1,l-1,k,l-1,l-1,l-1,l-1,i,l-1,l-1);
     endfor
   endfor
 endfor
@@ -100,13 +109,15 @@ endfor
 Ainv=A0inv;
 for l=1:nCols
   Ainv=(Id - ylk(:,l,l) * transpose(Id(:,l)) / (1 + ylk(l,l,l))) * Ainv;
+  P(:,l)=Id(:,p(l)); ## Construct permutation matrix to swap columns of Ainv
   % printf("Ainv=(Id - ylk(:,%d,%d) * transpose(Id(:,%d)) / (1 + ylk(%d,%d,%d))) * Ainv\n",l,l,l,l,l,l);
 endfor
+Ainv=transpose(P*transpose(Ainv)); ## Swap the 
 
 
 
 ## Test if the inverse found is really an inverse (does not work if values are floats)
-IdTest=A*Ainv
+IdTest=A*Ainv;
 if (IdTest==eye(nCols))
    printf("\n");
    printf("Inverse of A^{-1} FOUND!\n");