%     Gauss me meriki odigisi(xwris fysiki antimetathesi) 
      function x = gausswf(A,b);
%     A=[3 1 1; 2 1 -1;5 2 2]  
%     b=[-2 -3 -3]
      A=rand(10,10)
      x=[-1 2 2 -1 2 2 -1 2 2 3]';
      b=A*x;
%     A=[3 -2 3 4; 3 2 -4 1; -2 3 4 1; 6 -1 2 3];
      clc;
      disp(A);
%     b=[8 2 6 10];
      disp(b); 
      tv1 = cputime

      flops(0)

      [m,n] = size(A);
      if m ~=n
        disp (' matrix A is not square ');
        return
        end;
      p=[1:n];
      U = zeros(n,n);  
      for k = 1:n-1
         disp(p);
         M1 = eye(n,n);  
         d=p(k)
         s=A(d,k);

        for i = k+1:n
           if abs(A(p(i),k)) > abs(s)
              s = A(p(i),k);
              d = p(i);
           end;
        end;
        if s==0
            disp('the algorithm has encountered a zero pivot');
            x=[];
        return
        end;
        if d ~= k
        temp=p(k);
        p(k)=p(d);
        p(d)=temp;
    %   [p(k),p(d)] = intch(p(k),p(d));
        end;
            'inter'
        for i = k+1:n
           M1(i,k) = -A(i,k)/A(k,k);
           A(i,k) = M1(i,k);
        end;
        for j = k+1:n
          for i = k+1:n
               A(i,j) = A(i,j) + M1(i,k)*A(k,j);
          end
        end;
        for i = k+1:n
            b(i) = b(i) + M1(i,k)*b(k);
        end;
      end;
     U=triu(A);
     x=bcksub(U,b);
 
     tv2 = cputime - tv1;
     disp('the cpu time in seconds taken for this algorithm');

     tv2
     disp('the flop count for this algorithm is')
     flps = flops;

     end;

%____________________________________________________________
 
%    function [v1,v2] = intch(v1,v2)
%      temp=v1;
%      v1=v2;
%      v2=temp;
%    end;

%____________________________________________________________     
 
     function y = bcksub(T,b);
     [m,n]=size(T);
     if m~=n
     disp('matrix T is not square');
     return
     end;
     y = zeros(n,1);
      for i = n:-1:1
         sum = 0;
         for j = i+1:n
            sum = sum + T(i,j)*y(j);
         end;
         y(i) = (b(i) - sum) / T(i,i);
      end;
     end;       
 
 


   