高斯消元法解線性方程組 (高斯-塞得爾算法解方程組)

{輸入方程個數n,再輸入n個方程}
{例如方程3x+2y=1,2x+7y=0,則輸入:
2
3 2 1
2 7 0
}
{使用了高斯消元法}
program p1052;{解線性方程組}
const maxn=100;
type float=real;
     int=integer;
var a:array[1..maxn,1..maxn]of float;
    b:array[1..maxn]of float;
    d,t,sum:float;
    k,l,i,j,m,n,ii,jj,loop:int;
    ch:char;
begin
     readln(n);
     m:=n;
     for ii:=1 to n do begin
         for jj:=1 to m do read(a[ii][jj]);
         read(b[ii]);
     end;
     k:=1;
     d:=0;
     l:=0;
     while(k< =n)do begin
                   d:=a[k][k];
                   l:=k;
                   for i:=k+1 to n do begin
                       if(abs(a[i][k])>abs(d)) then begin
                                               d:=a[i][k];
                                               l:=i;
                       end;
                   end;
                   if(l<>k)then begin
                              for j:=k to n do begin
                                  t:=a[l,j];
                                  a[l,j]:=a[k,j];
                                  a[k,j]:=t;
                              end;
                              t:=b[k];
                              b[k]:=b[l];
                              b[l]:=t;
                   end;
                   for j:=k+1 to n do a[k,j]:=a[k][j]/a[k][k];
                   b[k]:=b[k]/a[k,k];
                   for i:=k+1 to n do begin
                       for j:=k+1 to n do a[i,j]:=a[i,j]-a[i,k]*a[k,j];
                       j:=1;
                       b[i]:=b[i]-a[i,k]*b[k];
                   end;
                   inc(k);
     end;
     for i:=n-1 downto 1 do begin
         sum:=0;
         for j:=i+1 to n do sum:=sum+a[i,j]*b[j];
         b[i]:=b[i]-sum;
     end;
     for loop:=1 to n do begin
         if loop=n then ch:=chr(13) else ch:=' ';
         write(round(b[loop]),ch);
     end;
     readln;
     readln;
end.

#include <cmath>
#include <cstring>
#include <algorithm>
/* Gaussian Elimination for Solving Linear Equations Set */
template <typename T>
inline T * gauss(T * * A, T * b, int m, T eps=1e-15)
{
   // [cof con']=[A B']
   int p, i, j;
   T * * cof=new T *[m];
   for (i=0; i<m ; i++) {
      cof[i]=new T[m];
      memcpy(cof[i], A[i], m*sizeof(T));
   }
   T * con=new T[m];
   memcpy(con, b, m*sizeof(T));
   T * & x=con;
   for (p=0; p<m-1; p++) {
      // find major element
      int max=p;
      for (j=p+1; j<m; j++) {
         if (fabs(cof[max][p])<fabs(cof[j][p])) {
            max=j;
         }
      }
      if (fabs(cof[max][p])<eps) {
         throw "Coefficient matrix is odd";
      }
      // apply row switching when necessary
      if (max!=p) {
         std::swap(cof[max], cof[p]);
         std::swap(con[max], con[p]);
      }
      // elimination
      for (i=p+1; i<m; i++) {
         T rate=cof[i][p]/cof[p][p];
         for (j=p; j<m; j++) {
            cof[i][j]-=rate*cof[p][j];
         }
         con[i]-=rate*con[p];
      }
   }
   // back-tracing the unknowns
   for (p=m-1; p>0; p--) {
      x[p]=con[p]/cof[p][p];
      for (j=1; j< =p; j++) {
         con[p-j]-=cof[p-j][p]*x[p];
         cof[p-j][p]=0.0;
      }
   }
   x[0]=con[0]/cof[0][0];
   delete [] cof;
   return x;
}