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
      87
      88
      89
      90
      91
      92
      93
      94
      95
      96
      97
      98
      99
     100
     101
     102
     103
     104
     105
     106
     107
     108
     109
     110
     111
     112
     113
     114
     115
     116
     117
     118
     119
     120
     121
     122
     123
     124
     125
     126
     127
     128
     129

#include <cassert>
#include <random>
#include <hpc/aux/walltime.h>
#include <hpc/matvec/copy.h>
#include <hpc/matvec/gematrix.h>
#include <hpc/matvec/mm.h>
#include <hpc/matvec/print.h>
#include <hpc/ulmlapack/getf2.h>

using namespace hpc::matvec;

template <typename T, typename I>
double
checkLU(const GeMatrixView<T, I> &A,
        GeMatrixView<T, I> &A_,
        const GeMatrixView<I, I> &p)
{
    assert(A.numRows==A_.numRows);
    assert(A.numCols==A_.numCols);
    assert(A.numCols>=A.numRows);
    assert(A.numRows==p.numRows);

    typedef typename GeMatrix<T, I>::Index  Index;

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

    // Copy L form A_ and set A_ to U
    GeMatrix<T, I> L(m, m);
    for (Index j=0; j<m; ++j) {
        for (Index i=0; i<m; ++i) {
            if (i==j) {
                L(i,i) = 1;
            } else if (i>j) {
                L(i,j)  = A_(i,j);
                A_(i,j) = 0;
            } else if (j>i) {
                L(i,j)  = 0;
            }
        }
    }

    // Compute LU = L*U
    GeMatrix<T, I> LU(m, n);
    for (Index j=0; j<n; ++j) {
        for (Index i=0; i<m; ++i) {
            LU(i,j) = 0;
        }
    }
    hpc::matvec::mm(T(1), L, A_, T(0), LU);

    // Apply P
    for (Index i=m; i>0; --i) {
        if (i-1!=p(i-1,0)) {
            hpc::ulmblas::swap(n,
                               &LU(i-1,0), LU.incCol,
                               &LU(p(i-1,0),0), LU.incCol);
        }
    }

    double diff = 0;
    for (Index j=0; j<n; ++j) {
        for (Index i=0; i<m; ++i) {
            diff += std::abs(A(i,j)-LU(i,j));
        }
    }
    return diff;
}

//
//  Random initializer for general matrices: real and complex valued
//
template <typename Index, typename T>
void
randomInit(Index m, Index n, T *A, Index incRowA, Index incColA)
{
    std::random_device                  random;
    std::default_random_engine          mt(random());
    std::uniform_real_distribution<T>   uniform(-100,100);

    for (Index i=0; i<m; ++i) {
        for (Index j=0; j<n; ++j) {
            A[i*incRowA+j*incColA] = uniform(mt);
        }
    }
}

int
main()
{
    const std::ptrdiff_t   m_min = 10;
    const std::ptrdiff_t   n_min = 10;

    const std::ptrdiff_t   m_max = 2000;
    const std::ptrdiff_t   n_max = 2000;

    const std::ptrdiff_t   m_inc = 10;
    const std::ptrdiff_t   n_inc = 10;

    GeMatrix<double, std::ptrdiff_t>      A(m_max, n_max);
    GeMatrix<double, std::ptrdiff_t>      A_org(m_max, n_max);
    GeMatrix<std::ptrdiff_t, std::ptrdiff_t>        p(m_max, 1);
    hpc::aux::WallTime<double>  timer;

    for (std::ptrdiff_t m=m_min, n=n_min; m<=m_max && n<=n_max; m+=m_inc, n+=n_inc) {

        randomInit(m, n, A.data, A.incRow, A.incCol);
        hpc::ulmblas::gecopy(m, n,
                             A.data, A.incRow, A.incCol,
                             A_org.data, A_org.incRow, A_org.incCol);
        //print(A_org, "A_org");

        timer.tic();
        std::ptrdiff_t info = hpc::ulmlapack::getf2(m, n,
                                                    A.data, A.incRow, A.incCol,
                                                    p.data, p.incRow);
        double t = timer.toc();
        //print(A, "A");
        //print(p, "p");

        auto A_ = A_org(0, 0, m, n);
        auto LU = A(0, 0, m, n);
        auto p_ = p(0, 0, m, 1);

        double diff = checkLU(A_, LU, p_);
        //print(A, "LU");
        std::printf("%4ld %4ld %4ld %16.3e %16.5lf\n", m, n, info, diff, t);
    }
}





































































































Home