Task

Solve   Ax=b   using Gaussian elimination then backwards substitution.

A   being an   n by n   matrix.

Also,   x and b   are   n by 1   vectors.

To improve accuracy, please use partial pivoting and scaling.



See also

   the Wikipedia entry:   Gaussian elimination


Solve: procedure options (main);    /* 11 January 2014 */
 
   declare n fixed binary;
   put ('Program to solve n simultaneous equations of the form Ax = b. Please type n:' );
   get (n);
 
begin;
   declare (A(n, n), b(n), x(n)) float(18);
   declare (SA(n,n), Sb(n)) float (18);
   declare i fixed binary;
 
   put skip list ('Please type A:');
   get (a);
   put skip list ('Please type the right-hand sides, b:');
   get (b);
 
   SA = A; Sb = b;
 
   put skip list ('The equations are:');
   do i = 1 to n;
      put skip edit (A(i,*), b(i)) (f(5), x(1));
   end;
 
   call Gauss_elimination (A, b);
 
   call Backward_substitution (A, b, x);
 
   put skip list ('Solutions:'); put skip data (x);
 
   /* Check solutions: */
   put skip list ('Residuals:');
   do i = 1 to n;
      put skip list (sum(SA(i,*) * x(*)) - Sb(i));
   end;
end;
 
Gauss_elimination: procedure (A, b) options (reorder); /* Triangularise */
   declare (A(*,*), b(*)) float(18);
   declare n fixed binary initial (hbound(A, 1));
   declare (i, j, k) fixed binary;
   declare t float(18);
 
   do j = 1 to n;
      do i = j+1 to n; /* For each of the rows beneath the current (pivot) row. */
         t = A(j,j) / A(i,j);
         do k = j+1 to n; /* Subtract a multiple of row i from row j. */
            A(i,k) = A(j,k) - t*A(i,k);
         end;
         b(i) = b(j) - t*b(i); /* ... and the right-hand side. */
      end;
   end;
end Gauss_elimination;
 
Backward_substitution: procedure (A, b, x) options (reorder);
   declare (A(*,*), b(*), x(*)) float(18);
   declare t float(18);
   declare n fixed binary initial (hbound(A, 1));
   declare (i, j) fixed binary;
 
   x(n) = b(n) / a(n,n);
 
   do j = n-1 to 1 by -1;
      t = 0;
      do i = j+1 to n;
         t = t + a(j,i)*x(i);
      end;
      x(j) = (b(j) - t) / a(j,j);
   end;
end Backward_substitution;
 
end Solve;