gausselimscaledpivot.m

A very simple Gaussian elimination code in matlab

function [x] = gaussElimScaledPivot(A,b, verbose)
  % File gaussElim.m
  %   This subroutine will perform Gaussian elmination
  %   on the matrix that you pass to it.
  %   i.e., given A and b it can be used to find x,
  %        Ax = b
  %
  %   To run this file you will need to specify several
  %   things:
  %   A - matrix for the left hand side.
    %   b - vector for the right hand side
      %
      %   The routine will return the vector x.
      %   ex: [x] = gaussElim(A,b)
      %     this will perform Gaussian elminiation to find x.
      %
      %


      N = max(size(A));
      L = zeros(N,N);
      P = eye(N,N);
      A0= A;
      % Perform Gaussian Elimination
      if verbose != 0
	A
	b
      end

      for k=1:N,
	s(k) = max(abs(A(k,:)));
      end
      
      s=s';
	

      for j=2:N,
	maxVal = A(j-1, j-1)/s(j-1);
	index = j-1;
	for k=j-1:N,
	  checking = abs(A(k, j-1))/s(k);
	  if checking > maxVal
	    maxVal = checking;
	    index  = k;
	  end
	end
	if index != j-1
	  A([index,j-1],:) = A([j-1,index],:);
	  b([index,j-1],:) = b([j-1,index],:);
	  s([index,j-1],:) = s([j-1,index],:);
	  L([index,j-1],:) = L([j-1,index],:);
	  P([index,j-1],:) = P([j-1,index],:);

	  if verbose != 0
	    printf("***** Swapping %d with %d val=%f******\n", index, j-1, maxVal);
	    A
	    b
	  end
	else
	  if verbose != 0
	    printf("***** No Swap ******\n");
	  end
	end
	
	for i=j:N,
	  m = myRound(myRound(A(i,j-1))/myRound(A(j-1,j-1)));
	  L(i, j-1) = m;
	  for k=1:N,
	    A(i,k) = myRound(myRound(A(i,k)) - myRound(myRound(A(j-1,k))*myRound(m)));
	  end
	  b(i) = myRound(myRound(b(i)) - myRound(m)*myRound(b(j-1)));
	  if verbose != 0
	    printf("****** Eliminating %d, %d ******\n", i, j-1);
	    A
	    b
	  end
	end  
      end

      % Perform back substitution

      x = zeros(N,1);
      x(N) = myRound(myRound(b(N))/myRound(A(N,N)));
      
      for j=N-1:-1:1,
	x(j) = myRound(myRound((myRound(b(j))-myRound(A(j,j+1:N)*x(j+1:N))))/myRound(A(j,j)));
      end
      

      A0;
      L=L+eye(N,N);
      U=A;
      P   ;  
      LU = L*U;
      PA = P*A0;

      
      % A = [.6667 .2857 .2000; .3333 .1429 -.5000; .2000 -.4286 .4000]
      % b = [2.867 .8333 -2.400]'

      % A = [2 3 1; 4 1 4; 3 4 6]
      % b = [-4 9 0]'