LU разложение неквадратной матрицы

На википедии LU разложение определено как разложение всякой матрицы в произведение двух квадратных матриц. Но как тогда сделать LU разложение неквадратной матрицы, и есть ли в этом смысл?

У меня есть следующий алгоритм, который работает для квадратных матриц, а также для матриц количество строк которых N < кол-ва столбцов M. Как можно его модифицировать, чтобы он работал для всех матриц?

for (size_t k = 0; k < N - 1; ++k) {
    for (size_t i = k + 1; i < N; ++i)
        if (std::abs(matrix[i][k]) > std::abs(matrix[k][k]))
            std::swap(matrix[i], matrix[k]);
    for (size_t i = k + 1; i < N; ++i) {
        matrix[i][k] /= matrix[k][k];
        for (size_t j = k + 1; j < M; ++j)
            matrix[i][j] -= matrix[i][k] * matrix[k][j];
    }
}