#ifndef HPC_ULMLAPACK_GETF2_H
#define HPC_ULMLAPACK_GETF2_H 1

#include <algorithm>
#include <hpc/ulmblas/ger.h>
#include <hpc/ulmblas/iamax.h>
#include <hpc/ulmblas/scal.h>
#include <hpc/ulmblas/swap.h>

namespace hpc { namespace ulmlapack {

template <typename Index, typename TA, typename TP>
Index
getf2(Index m, Index n, TA *A, Index incRowA, Index incColA, TP *p, Index incP)
{
    Index mn = std::min(m,n);

    for (Index j=0; j<mn; ++j) {
        Index i = j+ulmblas::iamax(m-j, &A[j*(incRowA+incColA)], incRowA);
        if (i!=j) {
            ulmblas::swap(n, &A[i*incRowA], incColA, &A[j*incRowA], incColA);
        }
        p[j*incP] = i;
        if (A[j*(incRowA+incColA)]==TA(0)) {
            return j;
        }
        ulmblas::scal(m-1-j, TA(1)/A[j*(incRowA+incColA)],
                      &A[j*(incRowA+incColA)+incRowA], incRowA);
        ulmblas::ger(m-1-j, n-1-j, TA(-1),
                     &A[j*(incRowA+incColA)+incRowA], incRowA,
                     &A[j*(incRowA+incColA)+incColA], incColA,
                     &A[(j+1)*(incRowA+incColA)], incRowA, incColA);
    }
    return -1;
}

} } // namespace ulmblas, hpc

#endif // HPC_ULMLAPACK_GETF2_H