lu_blk.h (flens/examples/lu

      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

#include <cxxstd/algorithm.h>

#include <flens/flens.h>

#include <flens/examples/lu_unblk.h>

#ifndef BS

#define BS 128

#endif

namespace flens {

template <typename MA>

typename GeMatrix<MA>::IndexType

lu_blk(GeMatrix<MA> &A)

{

    using flens::min;

    typedef typename GeMatrix<MA>::ElementType  ElementType;

    typedef typename GeMatrix<MA>::IndexType    IndexType;

    const IndexType  m = A.numRows();

    const IndexType  n = A.numCols();

    const IndexType  mn = min(m, n);

    const ElementType One(1);

    const Underscore<IndexType>  _;

    IndexType info;

    for (IndexType j=1; j<=mn; j+=BS) {

        const IndexType bs = min(BS, m-j, n-j);

        auto A11 = A(_(j,j+bs-1), _(j,j+bs-1));

        auto A12 = A(_(j,j+bs-1), _(j+bs,n));

        auto A21 = A(_(j+bs,m), _(j,j+bs-1));

        auto A22 = A(_(j+bs,m), _(j+bs,n));

        info = lu_unblk(A11);

        if (info) {

            return info+j-1;

        }

        // Compute A21*inv(A11.upper())

        blas::sm(Right, NoTrans, One, A11.upper(), A21);

        // Compute inv(A11.lowerUnit())*A12

        blas::sm(Left, NoTrans, One, A11.lowerUnit(), A12);

        A22 -= A21*A12;

    }

    return 0;

}

} // namespace flens