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
     66
     67
     68
     69
     70
     71
     72
     73
     74
     75
     76
     77
     78
     79
     80
     81
     82
     83
     84
     85
     86

#ifndef HPC_ULMLAPACK_GETF2_H
#define HPC_ULMLAPACK_GETF2_H 1

#include <algorithm>
#include <hpc/matvec/r.h>
#include <hpc/matvec/iamax.h>
#include <hpc/matvec/scal.h>
#include <hpc/matvec/swap.h>

#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;
}

template <typename MA, typename VP>
typename std::enable_if<hpc::matvec::IsGeMatrix<MA>::value
                     && hpc::matvec::IsDenseVector<VP>::value,
         typename std::remove_reference<MA>::type::Index>::type
tf2(MA &&A, VP &&p)
{
    typedef typename std::remove_reference<MA>::type MatrixA;

    typedef typename MatrixA::ElementType    T;
    typedef typename MatrixA::Index          Index;

    using namespace hpc::matvec;

    Index m = A.numRows;
    Index n = A.numCols;


    for (Index j=0; m>0 && n>0; ++j, --m, --n) {
        auto Ajj = A(j, j, m, n);

        Index i = j+iamax(Ajj.col(0));
        if (i!=j) {
            swap(A.row(i), A.row(j));
        }
        p(j) = i;
        if (A(j,j)==T(0)) {
            return j;
        }
        if (m==1 || n==1) {
            break;
        }

        auto a01 = Ajj.row(0)(1,n-1);
        auto a10 = Ajj.col(0)(1,m-1);
        auto A11 = Ajj(1, 1, m-1, n-1);

        scal(T(1)/A(j,j), a10);
        r(T(-1), a10, a01, A11);
    }
    return -1;
}

} } // namespace ulmblas, hpc

#endif // HPC_ULMLAPACK_GETF2_H





































































































Home