% version 1.2
n = 1000
maxval = 10
A = rand(n,n) * maxval;
for i = 1:n
	A(i,i) = rand*maxval + sum(A(i,:)) - A(i,i);
endfor
xexact = rand(n,1) * maxval - maxval/2 ;
b = A * xexact ;
D = diag(diag(A)) ;
B = D - A ;
H = inv(D) * B ; 
d = inv(D) * b ;
fdelta = max(abs(H * xexact + d - xexact)) ;
disp("Solution delta:"), disp(fdelta);
disp("Matrix norm of H: "), disp(norm(H));
disp("Condition number of A: "), disp(cond(A));
save "a.txt" A
save "b.txt" b
save "h.txt" H
save "d.txt" d
save "xexact.txt" xexact
