(subscriptrange, fofl, size):                         /* 2 Nov. 2013 */
LU_Decomposition: procedure options (main);
   declare a1(3,3) float (18) initial ( 1,  3,  5,
                                        2,  4,  7,
                                        1,  1,  0);
   declare a2(4,4) float (18) initial (11,  9, 24, 2,
                                        1,  5,  2, 6,
                                        3, 17, 18, 1,
                                        2,  5,  7, 1);
   call check(a1);
   call check(a2);


/* In-situ decomposition */
LU: procedure(a, p);
   declare a(*,*) float (18);
   declare p(*)   fixed binary;
   declare (maximum, rtemp) float (18);
   declare (n, i, j, k, ii, temp) fixed binary;

   n = hbound(a,1);
   do i = 1 to n; p(i) = i; end;

   do k = 1 to n-1;

      maximum = 0; ii = k;
      do i = k to n;
         if maximum < abs(a(p(i),k)) then
            do; maximum = abs(a(p(i),k)); ii = i; end;
      end;
      if ii ^= k then do; temp = p(k); p(k) = p(ii); p(ii) = temp; end;

      do i = k+1 to n; a(p(i),k) = a(p(i),k) / a(p(k),k); end;

      do j = k+1 to n;
         do i = k+1 to n;
            a(p(i),j) = a(p(i),j) - a(p(i),k) * a(p(k),j);
         end;
      end;

   end;
end LU;

CHECK: procedure(a);
   declare a(*,*) float (18) nonassignable;

   declare aa(hbound(a,1), hbound(a,2)) float (18);
   declare  L(hbound(a,1), hbound(a,2)) float (18);
   declare  U(hbound(a,1), hbound(a,2)) float (18);
   declare (p(hbound(a,1), hbound(a,2)), ipiv(hbound(a,1)) ) fixed binary;
   declare pp(hbound(a,1), hbound(a,2)) fixed binary;
   declare (i, j, n, temp(hbound(a,1))) fixed binary;

   n = hbound(a,1);
   aa = A; /* work with a copy */
   P = 0; L = 0; U = 0;
   do i = 1 to n;
      p(i,i) = 1; L(i,i) = 1; /* convert permutation vector to a matrix */
   end;

   call LU(aa, ipiv);

   do i = 1 to n;
      do j = 1 to i-1; L(i,j) = aa(ipiv(i),j); end;
      do j = i to n;   U(i,j) = aa(ipiv(i),j); end;
   end;

   pp = p;
   do i = 1 to n;
      p(ipiv(i), *) = pp(i,*);
   end;

   put skip list ('A');
   put edit (A) (skip, (n) f(10,5));

   put skip list ('P');
   put edit (P) (skip, (n) f(11));

   put skip list ('L');
   put edit (L) (skip, (n) f(10,5));

   put skip list ('U');
   put edit (U) (skip, (n) f(10,5));

end CHECK;

end LU_Decomposition;