open Array;

     fun pivot(M,m,n) =
             if m>=n then 0.0
             else if m=n-1 then sub(sub(M,m),m)
             else (* 0<m<n-1 *) (
                 let (* normalize  row m to right of diagonal *)
                     val i = ref (m+1);
                     val a = sub(sub(M,m),m)
                 in
                     while !i<n do (
                         update(sub(M,m), !i, sub(sub(M,m),!i)/a);
                         i := !i + 1
                     )
                 end;
                 let (* subtract M[i,m]*M[m,j] from M[i,j]
                        for all i and j such that m<i,j<n *)
                     val i = ref (m+1);
                     val j = ref (m+1)
                 in
                     while !i<n do (
                         while !j<n do (
                             update(sub(M,!i),!j,sub(sub(M,!i),!j)-
                                 sub(sub(M,!i),m)*sub(sub(M,m),!j));
                             j := !j + 1
                         );
                         j := m + 1;
                         i := !i + 1
                     )
                 end;
                 (* result is M[m,m] times recursive call on matrix
                    starting at row and column m+1 *)
                 sub(sub(M,m),m)*pivot(M,m+1,n)
             );