diff --git a/Matlab/Lectures/Lecture22/lecture22a.m b/Matlab/Lectures/Lecture22/lecture22a.m
index 8fdfa73910be4c6248fcb50b462a890cbadf5b5c..668e132bff69d068252cd116554d8ba920f160d6 100644
--- a/Matlab/Lectures/Lecture22/lecture22a.m
+++ b/Matlab/Lectures/Lecture22/lecture22a.m
@@ -15,28 +15,36 @@ cd Lectures/Lecture22;
close all
+example = 3;
+switch example
+ case 1 % sparse random matrix which we make diagonally dominant
+ rng(42);
+ n = 100;
+ A = sprand(n,n,0.1);
+ A = A +spdiags( 2*n*ones(n,1), 0,n,n);
+ disp('sparse random');
+ case 2 %from SuiteSparse, small
+ A = mmread('fs_183_6.mtx');
+ disp('fs_183_6');
+ case 3 %SuiteParse
+ A = mmread('Goodwin_040.mtx');
+ disp('Goodwin_040');
+ case 4
+ A = mmread('gre_1107.mtx');
+ disp('gre_1107');
+ otherwise
+ disp('Example does not exist')
+end
+
+ n = size(A,1);
+ fprintf('dimension is %d\n',n);
+ figure
+ subplot(1,2,1)
+ spy(A);
+ title('A')
-
-indx = [ 1 3 8 10];
-valx = [-1 2 1 -2];
-indy = [1 3 8];
-valy = [ 1 2 -3];
-%alpha = sp_inner_prod(indx,valx,indy,valy)
-
-%graph
-% 1 2 3 4 5 6 7 8 8 10
-% colindG = [1 7 1 5 1 4 5 3 6 7];
-% rowptrG= [1 3 5 8 11 11 11 11];
-% n = 7;
-% x = [ 1 4]
-% y = reach(rowptrG,colindG,x,n)
-
-rng(42);
-n = 100;
-A = sprand(n,n,0.1);
-A = A +spdiags( 2*n*ones(n,1), 0,n,n);
-spy(A);
+ pause
[rowindA,colindA,valA] = find(A);
[rowptrA,colindA,valA] = coord2csr(rowindA,colindA,valA,n);
@@ -49,7 +57,8 @@ L = sparse(rowindL,colindL,valL,n,n);
[rwoindU,colindU, valU] = csc2coord(rowindU,colptrU,valU,n);
U = sparse(rowindU,colindU,valU,n,n);
-figure
+subplot(1,2,2)
spy(L+U)
+title('L+U')
acc = norm(A-L*U,'fro')
\ No newline at end of file
diff --git a/Matlab/functions/__graph/reach.m b/Matlab/functions/__graph/reach.m
index 6e45c99a79e9795f3ea353885c469090181294ba..83ca29e09f78e7eeb17d9929b2b172ef1f8d7f93 100644
--- a/Matlab/functions/__graph/reach.m
+++ b/Matlab/functions/__graph/reach.m
@@ -1,4 +1,4 @@
-function y = reach(rowptrG,colindG,s,n);
+function y = reach(rowptrG,colindG,s,n)
% returns all reachable nodes y(i) from node set s in directed graph [rowptrG,colindG]
% does BFS from all nodes in x; keeps the marks so that we can abort early
% when using later nodes from x
diff --git a/Matlab/functions/__lu_symbolic/lu_sparse.asv b/Matlab/functions/__lu_symbolic/lu_sparse.asv
deleted file mode 100644
index df8d6bbde3425d61c1540b2056d897b78e5efdbb..0000000000000000000000000000000000000000
--- a/Matlab/functions/__lu_symbolic/lu_sparse.asv
+++ /dev/null
@@ -1,81 +0,0 @@
-function [rowindU,colptrU,valU,rowptrL,colindL,valL] = lu_sparse(rowptrA,colindA,valA,n)
-% computes the sparseLU factorization according to p. 96 in the book
-% no symbolic phase, we expand the structure for L and U as we go
-% Andreas Frommer, 07.01.2024
-
-% we need the csc and csr format of A, so we convert. And we make sure that
-% indices for A are sorted in both formats by converting again
-
-[rowindA,colptrA,valA_csc] = csr2csc(rowptrA,colindA,valA,n);
-[rowptrA,colindA,valA_csr] = csc2csr(rowindA,colptrA,valA_csc,n);
-
-% Initialize first diagonal elelement in L and U
-% U is csc, L is csr
-colptrU(1) = 1;
-colptrU(2) = 2;
-rowindU(1) = 1;
-valU(1) = valA_csc(1);
-
-rowptrL(1) = 1;
-rowptrL(2) = 2;
-colindL(1) = 1;
-valL(1) = 1;
-
-for k=2:n
- % compute row k of L
- % At this point, rowptrU, colindU, valU holds U(1:k-1,1:k-1)
- % We solve U^T(1:k-1,1:k-1)y = A^T(k,1:k-1)
- % Then y^T = L(k,1:k-1)
- % The csr format of U^T is [colptrU, rowindU, valU]
-
- % we need row k of A left to the diagonal (excluded)
- rowind = colindA(rowptrA(k):rowptrA(k+1)-1); % whole row
- l = 0;
- while (rowind(l+1) < k)
- l = l+1;
- end
- rowind = rowind(1:l); % take care of empty rows (l=0) later
- val = valA_csr(rowptrA(k):rowptrA(k)+l-1);
-
-
- [rowindy,valy] = forward_solve_sparse_sparse(colptrU,rowindU,valU,rowind,val,k-1);
- rowindy = [rowindy,k]; valy = [valy, 1]; % add diagonal
- rowptrL(k+1) = rowptrL(k) + length(rowindy);
- colindL(rowptrL(k):rowptrL(k+1)-1) = rowindy;
- valL(rowptrL(k):rowptrL(k+1)-1) = valy;
-
- % now compute column k of U
- % At this point, rowindL, colptrL, valL holds L(1:k,1:k), but we only need
- % L(1:k-1,1:k-1)
- % We solve L(1:k-1,1:k-1)y = A(1:k-1,k)
- % Then U(1:k-1,k) = y
-
- % we need column k of A above the diagonal (excluded)
- % % and the value of the diagonal element valA_kk for later use
-
- rowind = rowindA(colptrA(k):colptrA(k+1)-1); % whole column
- l = 0;
- while (rowind(l+1) < k)
- l = l+1;
- end
- rowind = rowind(1:l); % take care of empty rows (l=0) later
- val = valA_csc(colptrA(k):colptrA(k)+l-1);
- val_Akk = valA_csc(colptrA(k)+l);
-
- [rowindy,valy] = forward_solve_sparse_sparse(rowptrL,colindL,valL,rowind,val,k-1);
- rowindy = [rowindy,k]; valy = [valy, 1]; % add diagonal; 1 must be corrected later
- colptrU(k+1) = colptrU(k) + length(rowindy);
- rowindU(colptrU(k):colptrU(k+1)-1) = rowindy;
- valU(colptrU(k):colptrU(k+1)-1) = valy;
-
- % now correct the diagonal, we need the inner product of L(k,1:k-1)
- % with U(1:k-1,k). Both are sorted
- indu = rowindU(colptrU(k):colptrU(k+1)-2);
- valu = valU(colptrU(k):colptrU(k+1)-2);
- indl = colindL(rowptrL(k):rowptrL(k+1)-2);
- vall = valL(rowptrL(k):rowptrL(k+1)-2);
-
- valU(colptrU(k+1)-1) = val_Akk-sp_inner_prod(indu,valu,indl,vall);
-
-end
-end
\ No newline at end of file
diff --git a/Matlab/functions/__lu_symbolic/lu_sparse.m b/Matlab/functions/__lu_symbolic/lu_sparse.m
index e56c4fda35350b61e8deccc90b12c93779f2c9cf..18aabc0223e499a823a51ad2765f672f16cd06eb 100644
--- a/Matlab/functions/__lu_symbolic/lu_sparse.m
+++ b/Matlab/functions/__lu_symbolic/lu_sparse.m
@@ -29,15 +29,15 @@ for k=2:n
% The csr format of U^T is [colptrU, rowindU, valU]
% we need row k of A left to the diagonal (excluded)
- rowind = colindA(rowptrA(k):rowptrA(k+1)-1); % whole row
+ rowk = colindA(rowptrA(k):rowptrA(k+1)-1); % whole row
l = 0;
- while (rowind(l+1) < k)
+ while ( (rowk(l+1) < k) & (l+1 < length(rowk) ) )
l = l+1;
end
- rowind = rowind(1:l); % take care of empty rows (l=0) later
+ rowk = rowk(1:l); % take care of empty rows (l=0) later
val = valA_csr(rowptrA(k):rowptrA(k)+l-1);
- [rowindy,valy] = forward_solve_sparse_sparse(colptrU,rowindU,valU,rowind,val,k-1);
+ [rowindy,valy] = forward_solve_sparse_sparse(colptrU,rowindU,valU,rowk,val,k-1);
rowindy = [rowindy; k]; valy = [valy; 1]; % add diagonal
rowptrL(k+1) = rowptrL(k) + length(rowindy);
colindL(rowptrL(k):rowptrL(k+1)-1) = rowindy;
@@ -52,17 +52,17 @@ for k=2:n
% we need column k of A above the diagonal (excluded)
% % and the value of the diagonal element valA_kk for later use
- rowind = rowindA(colptrA(k):colptrA(k+1)-1); % whole column
+ rowk = rowindA(colptrA(k):colptrA(k+1)-1); % whole column
l = 0;
- while (rowind(l+1) < k)
+ while (( rowk(l+1) < k) & l+1 < length(rowk))
l = l+1;
end
- rowind = rowind(1:l); % take care of empty rows (l=0) later
+ rowk = rowk(1:l); % take care of empty rows (l=0) later
val = valA_csc(colptrA(k):colptrA(k)+l-1);
val_Akk = valA_csc(colptrA(k)+l);
% we work with L(1:k-1,1:k-1), note that k-th row has just been computed
- [rowindy,valy] = forward_solve_sparse_sparse(rowptrL(1:k),colindL(1:rowptrL(k)-1),valL,rowind,val,k-1);
+ [rowindy,valy] = forward_solve_sparse_sparse(rowptrL(1:k),colindL(1:rowptrL(k)-1),valL,rowk,val,k-1);
rowindy = [rowindy,; k]; valy = [valy; 1]; % add diagonal; 1 must be corrected later
colptrU(k+1) = colptrU(k) + length(rowindy);
rowindU(colptrU(k):colptrU(k+1)-1) = rowindy;
diff --git a/Matlab/functions/__matrix_data/Goodwin_040.mtx b/Matlab/functions/__matrix_data/Goodwin_040.mtx
new file mode 100644
index 0000000000000000000000000000000000000000..2305f6d0d7878b676edf34cab16c56523d93864e
--- /dev/null
+++ b/Matlab/functions/__matrix_data/Goodwin_040.mtx
@@ -0,0 +1,561695 @@
+%%MatrixMarket matrix coordinate real general
+%-------------------------------------------------------------------------------
+% SuiteSparse Matrix Collection, Tim Davis
+% https://sparse.tamu.edu/Goodwin/Goodwin_040
+% name: Goodwin/Goodwin_040
+% [Finite element, Navier-Stokes &