Is the U factor in LU decomposition for rectangular matrices always in row echelon form? I have come across the following rectangular 5 x 10 matrix and carried out a LU decomposition of it, in the form PA = LU. The following matrices were obtained by function scipy.linalg.lu from module scipy to Python language. At first I thought was a computational issue, but after a contributor obtained the same result in both matlab and R, I decided to post it here (See: https://stackoverflow.com/questions/25186560/is-the-upper-triangular-matrix-in-function-scipy-linalg-lu-always-in-row-echelon)
\begin{equation} 
    A = \begin{bmatrix}
        1 & 1 & 1 & 1 & 0 & 1 & 1 & 1 & 1 & 0 \\
        1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 1 & 1 \\
        1 & 1 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 0 \\
        0 & 1 & 0 & 1 & 1 & 0 & 1 & 0 & 0 & 1 \\
        1 & 1 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 \end{bmatrix}
\end{equation}
\begin{equation}
    U = \begin{bmatrix}
        1 &  1 &  1 &  1 &  0 &  1 &  1 &  1 &  1 & 0 \\
        0 &  1 &  0 &  1 &  1 &  0 &  1 &  0 &  0 &  1 \\
        0 &  0 & -1 & -1 &  0 &  0 &  0 & -1 & -1 &  0 \\
        0 &  0 &  0 &  0 &  1 & -1 &  0 &  0 &  1 &  1 \\
        0 &  0 &  0 &  0 &  1 &  0 &  0 &  1 &  1 &  1
      \end{bmatrix}
\end{equation}
Notice that U is upper triangular, though it is not in row echelon form. Element U[5,5] = 1, but it should be zero, since the pivot for row 4 is immediately above it. According to Strang (1988, p.72), the U factor of a LU decomposition of a rectangular matrix should be in row echelon form, with the following characteristics:
"(i) The nonzero rows come first - otherwise there would have been row exchanges - and pivots are the first nonzero entries in those rows. (ii) Below each pivot is a column of zeros, obtained by elimination. (iii) Each pivot lies to the right of the pivot in the row above. This produces the staircase pattern."
STRANG, G. Linear Algebra and its applications. 3rd Ed., Thomson, 1988.
The L factor and permutation matrix P are below. Notice in P that Gaussian elimination swapped rows 2 and 4.
\begin{equation}
    L = \begin{bmatrix}
         1 &  0 &  0 &  0 &  0 \\ 
         0 &  1 &  0 &  0 &  0 \\
         1 &  0 &  1 &  0 &  0 \\
         1 &  0 &  1 &  1 &  0 \\
         1 &  0 &  1 &  0 &  1
    \end{bmatrix}
\end{equation}
\begin{equation}
    P = \begin{bmatrix}
         1 &  0 &  0 &  0 &  0 \\
         0 &  0 &  0 &  1 &  0 \\
         0 &  0 &  1 &  0 &  0 \\
         0 &  1 &  0 &  0 &  0 \\
         0 &  0 &  0 &  0 &  1 \\
    \end{bmatrix}
\end{equation}
Say one tries to cure the problem by cancelling out element U[5,5] = 1, by subtracting row 5 from row 4, and finally producing the row echelon form. L will not be modified, since L[5,5] will remain equals 1, and $PA \neq LU$. I suspect that Gaussian elimination does not always produce the row echelon form.
 A: 
According to Strang (1988, p.72), the U factor of a LU decomposition of a rectangular matrix should be in row echelon form,

What you can conclude from this is that in the 1988 book by Strang, the LU decomposition is defined so that $U$ is in row echelon form. By picking up another linear algebra book, you may find LU described differently. Also, none of those books constitute the documentation of numerical algebra routines. 
The LAPACK routine dgetrf (which is probably used by the aforementioned software) promises the following output: 

The factorization has the form $A = P  L  U$  where $P$ is a permutation matrix, $L$ is lower triangular with unit  diagonal elements (lower trapezoidal if $m > n$), and $U$ is upper  triangular (upper trapezoidal if $m < n$).

No promise is made to produce row-echelon form. A simpler example than yours is 
$$
A = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 \end{pmatrix}
$$
for which the aforementioned software returns $U=A$.
This may happen for square matrices, too: e.g., $U=A$ is also returned for 
$$
A = \begin{pmatrix} 1 & 0 & 0  \\ 0 & 0 & 1  \\ 0 & 0 & 1  \end{pmatrix}
$$
