
randn('state',sum(100*clock));

disp(' ')
disp(' ')
n = 2000;
disp(['Generate random data for A,x and b=Ax, with n=  ',num2str(n)])
%A = randn(n);   
 A = sprandn(n,n,.001) + diag(rand(n,1));% or use a sparse example
xhat = randn(n,1);
b = A*xhat;

%condest(A)


disp(' ')
disp('Solve by x = inv(A)*b:');
pause(.1)

tic, x = inv(A)*b; toc
res = norm(A*x-b); disp(['residual ',num2str(res)]);
err = norm(x-xhat);disp(['error ',num2str(err)]);


disp(' ')
disp('Solve by x = A\b:');
pause(.1)

tic, x = A\b; toc
res = norm(A*x-b);  disp(['residual ',num2str(res)]);
err = norm(x-xhat); disp(['error ',num2str(err)]);


disp(' ')
disp('Solve by [L,U,P]=lu(A); x=U\(L\b):');
pause(.1)

tic, [L,U,p] = lu(A,'vector'); x=U\(L\b(p)); toc
res = norm(A*x-b);   disp(['residual ',num2str(res)]);
err = norm(x-xhat);  disp(['error ',num2str(err)]);


disp(' ')
disp('linsolve with [L,U,P]=lu(A):')
pause(.1)

tic
[L,U,p] = lu(A,'vector');
if issparse(L), L=full(L); U=full(U); end
y = linsolve(L,b(p),struct('LT',true));
x = linsolve(U,y,struct('UT',true));
toc
res = norm(A*x-b);    disp(['residual ',num2str(res)]);
err = norm(x-xhat);   disp(['error ',num2str(err)]);

