Home Browse Files 
      1
      2
      3
      4
      5
      6
      7
      8
      9
     10
     11
     12
     13
     14
     15
     16
     17
     18
     19
     20
     21
     22
     23
     24
     25
     26
     27
     28
     29
     30
     31
     32
     33
     34
     35
     36
     37
     38
     39
     40
     41
     42
     43
     44
     45
     46
     47
     48
     49
     50
     51
     52
     53
     54
     55
     56
     57
     58
     59
     60
     61
     62
     63
     64
     65

#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





































































































Home