Public paste
sdsd
By: sdsd | Date: Jan 6 2012 20:43 | Format: None | Expires: never | Size: 1.74 KB | Hits: 786

  1. //Stefan Domański
  2. function [singular, x] = gaussian(A,b)
  3.     M = A;
  4.     w = b;
  5.     osob = 0;
  6.     perm = 1;
  7.    
  8.     for i=2:length(b)
  9.         perm(i)=i;
  10.     end
  11.    
  12.     for i=1:length(b)
  13.         max=0;
  14.         xm=0;
  15.         ym=0;
  16.         for x=i:length(b)
  17.             for y=i:length(b)
  18.                 if abs(M(x, y))>max then
  19.                     max = abs(M(x, y))
  20.                     xm = x
  21.                     ym = y
  22.                 end
  23.             end
  24.         end
  25.         if max<1E-1 then
  26.             osob = 1;
  27.         end
  28.         for j=1:length(b)
  29.             temp = M(j, ym);
  30.             M(j, ym) = M(j, i);
  31.             M(j, i) = temp;
  32.         end
  33.         temp = w(xm);
  34.         w(xm)= w(i);
  35.         w(i) = temp;
  36.        
  37.         temp = perm(xm);
  38.         perm(xm)= perm(i);
  39.         perm(i) = temp;
  40.         for j=1:length(b)
  41.             temp = M(xm, j);
  42.             M(xm, j) = M(i, j);
  43.             M(i, j) = temp;
  44.         end
  45.        
  46.         for j= i+1:length(b)
  47.             ratio = M(j,i)/M(i,i);
  48.             for k=i:length(b)
  49.                 M(j,k) = M(j, k) - ratio*M(i, k);
  50.             end
  51.             w(j) = w(j) - ratio*w(i);
  52.         end
  53.     end
  54.     for i=1:length(b)
  55.         j = length(b)-i+1;
  56.         r = 0;
  57.         for k=j+1:length(b)
  58.             r = r + M(j, k)*w(k);
  59.         end
  60.         w(j)=w(j)/M(j, j);
  61.     end
  62.     for i=1:length(b)
  63.         temp = perm(i);
  64.         while temp <> i
  65.             temp1 = w(i);
  66.             w(i)=w(perm(i));
  67.             w(perm(i))=temp1;
  68.            
  69.             temp1 = perm(i);
  70.             perm(i)=perm(temp);
  71.             perm(temp)=temp1;
  72.            
  73.             temp=perm(i);
  74.         end
  75.     end
  76.    
  77.     singular = osob;
  78.     x = w;
  79. endfunction