%
% testlu_cp.m
%
%

example = input(' Example 1, 2, 3 (enter 0 to exit )? ')

if( example == 1 )
    A      = [ 2 4 2
               4 6 8
               1 3 2 ];
    b      = [  4; 12; 3 ];
    Xexact = [ 1; 0; 1 ];
elseif( example == 2 )

    A      = [ 1 2 2;
               4 4 2;
               4 6 4];
    b      = [ 3; 6 ; 10];
    Xexact = [-1; 3; -1 ];

elseif( example == 3 )

    % Now the same with random matrices
    n = 5;
    A=rand(n,n);       %  A is "almost surely" nonsingular

    Xexact=rand(n,1);  % Generate an "exact solution"
    b=A*Xexact;        %  and the right-hand side that will produce it

else
    break;
end


disp('  The system matrix A ')
disp( A )
disp( ' and the right hand side vector b ')
disp( b )

[A, ipivtr, ipivtc, iflag] = lu_cp( A );
disp( ' The output A of lu_cp ')
disp( A  )
disp( ' ipivtr ')
disp( ipivtr  )
disp( ' ipivtc ')
disp( ipivtc  )

[b, iflag] = lu_cp_sl( A, b, ipivtr, ipivtc );
disp( ' The solution of A x = b computed using lu_cp and lu_cp_sl')
disp( b )

fprintf(1,'   Error between exact and computed solution = %12.8e \n', ...
          norm(b-Xexact))