#ifndef HPC_ULMLAPACK_GETRF_H #define HPC_ULMLAPACK_GETRF_H 1 #include <algorithm> #include <hpc/ulmblas/gemm.h> #include <hpc/ulmblas/trlsm.h> #include <hpc/ulmlapack/getf2.h> #include <hpc/ulmlapack/swap.h> namespace hpc { namespace ulmlapack { template <typename Index, typename TA, typename TP> Index getrf(Index m, Index n, TA *A, Index incRowA, Index incColA, TP *p, Index incP) { Index mn = std::min(m,n); Index bs = 64; Index info = -1; if (bs<=1 || bs>mn) { return getf2(m, n, A, incRowA, incColA, p, incP); } else { for (Index j=0; j<mn; j+=bs) { Index jb = std::min(mn-j, bs); Index info_ = getf2(m-j, jb, &A[j*(incRowA+incColA)], incRowA, incColA, &p[j*incP], incP); if (info==-1 && info_>-1) { info = j + info_; } for (Index k=j; k<j+jb; ++k) { p[k*incP] += j; } swap(m, j, A, incRowA, incColA, j, j+jb, p, incP); if (j+jb<n) { swap(m, n-(j+jb), &A[(j+jb)*incColA], incRowA, incColA, j, j+jb, p, incP); ulmblas::trlsm(jb, n-(j+jb), TA(1), true, &A[j*(incRowA+incColA)], incRowA, incColA, &A[j*incRowA+(j+jb)*incColA], incRowA, incColA); } if (j+jb<n && j+jb<m) { ulmblas::gemm(m-(j+jb), n-(j+jb), jb, TA(-1), &A[(j+jb)*incRowA+j*incColA], incRowA, incColA, &A[j*incRowA+(j+jb)*incColA], incRowA, incColA, TA(1), &A[(j+jb)*(incRowA+incColA)], incRowA, incColA); } } } return info; } } } // namespace ulmblas, hpc #endif // HPC_ULMLAPACK_GETRF_H |