%     Block Gauss-Jordan (KJI-morfh) 
%     function x = gausswf(A,b);
%     A=[3 1 1; 2 1 -1;5 2 2]  
%     b=[-2 -3 -3]
%     A=[3 -2 3 4; 3 2 -4 1; -2 3 4 1; 6 -1 2 3];
%     b=[8 2 6 10];
      A=rand(8,8);
      disp(A);
      x=[-1 2 2 3 -1 2 2 3]';
      disp(x);
      b=A*x;

      disp(b); 

      tv1 = cputime

      flops(0)

      [m,n] = size(A);
      if m ~=n
        disp (' matrix A is not square ');
        return
      end;
      r=4;
      q=n/r;  
      for k = 1:q
           ik1=r*(k-1)+1;
           ik2=ik1+r-1
           A1=inv(A([ik1:ik2],[ik1:ik2]))
        for j = k:q
           ic1=r*(j-1)+1;
           ic2=ic1+r-1;
           A([ik1:ik2],[ic1:ic2])= A1*A([ik1:ik2],[ic1:ic2]);
           for i = 1:q
              if i ~= k 
                il1=r*(i-1)+1;
                il2=il1+r-1;
                A2=A([il1:il2],[ik1:ik2])*A([ik1:ik2],[ic1:ic2]);
                A3=A([il1:il2],[ic1:ic2])-A2;
                A([il1:il2],[ic1:ic2])=A3; 
              end;
           end;
        end;
      end;
      disp(A);     
     
      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;
 


   