Task

Invert matrix   A   using Gauss-Jordan method.

A   being an   n × n   matrix.




/* Gauss-Jordan matrix inversion */

G_J: procedure options (main);       /* 4 November 2020 */

   declare t float;
   declare (i, j, k, n) fixed binary;

   open file (sysin) title ('/GAUSSJOR.DAT');

   get (n);                     /* Read in the order of the matrix. */
   put skip data (n);
   begin;
      declare a(n,n)float, aux(n,n) float;

      aux = 0;
      do i = 1 to n; aux(i,i) = 1; end;
      get (a);                  /* Read in the matrix. */
      put skip list ('The matrix to be inverted is:');
      put edit (a) ( skip, (n) F(10,4)); /* Print  the matrix. */

      do k = 1 to n;
         /* Divide row k by a(k,k) */
         t = a(k,k); a(k,*) = a(k,*) / t; aux(k,*) = aux(k,*) / t;

         do i = 1 to k-1, k+1 to n; /* Work down the rows. */
            t = a(i,k);
            a(i,*) = a(i,*) - t*a(k,*); aux(i,*) = aux(i,*)  - t*aux(k,*);
         end;
      end;
      put skip (2) list ('The inverse is:');
      put edit (aux) ( skip, (n) F(10,4));
   end; /* of BEGIN block */
end G_J;
