new; cls;
k=5;   /* k x k matrix */
x=1000*rndu(k,k);
print "x=" x;
print "x_1= " inv(x); /* by built-in function inv(x) */
print "x_1'=" sinv(x);


proc sinv(x);
   local k,i,j,l,y,a;
   k=rows(x);
   y=eye(k);
   j=1;
   do while j<=k;
      i=1;
      do while i<=k;
         if i==j;
            a=1/x[i,j];
            l=1;
            do while l<=k;
               x[i,l]=a*x[i,l];
               y[i,l]=a*y[i,l];
               l=l+1;
            endo;
         endif;
         i=i+1;
      endo;
      i=1;
      do while i<=k;
         if i/=j;
            a=-x[i,j]/x[j,j];
            l=1;
            do while l<=k;
               x[i,l]=a*x[j,l]+x[i,l];
               y[i,l]=a*y[j,l]+y[i,l];  
               l=l+1;
            endo;
         endif;
         i=i+1;
      endo;
      j=j+1;
   endo;
   retp(y);
endp;