#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