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
     130
     131
     132
     133
     134
     135
     136
     137
     138
     139
     140
     141
     142
     143
     144
     145
     146
     147
     148
     149
     150
     151
     152
     153

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

using namespace hpc::matvec;

template <typename MA, typename MA_, typename VP>
typename std::enable_if<hpc::matvec::IsGeMatrix<MA>::value
                     && hpc::matvec::IsGeMatrix<MA_>::value
                     && hpc::matvec::IsDenseVector<VP>::value,
         double>::type
checkLU(const MA &A, MA_ &A_, const VP &p)
{
    assert(A.numRows==A_.numRows);
    assert(A.numCols==A_.numCols);
    assert(A.numCols>=A.numRows);
    assert(A.numRows==p.length);

    typedef typename MA::ElementType  T;
    typedef typename MA::Index        I;
    typedef typename MA::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)) {
            hpc::ulmblas::swap(n,
                               &LU(i-1,0), LU.incCol,
                               &LU(p(i-1),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);
        }
    }
}

template <typename MA>
typename std::enable_if<hpc::matvec::IsGeMatrix<MA>::value,
                        void>::type
randomInit(MA &A)
{
    randomInit(A.numRows, A.numCols, A.data, A.incRow, A.incCol);
}


int
main()
{
    using namespace hpc::matvec;

    typedef double       T;
    typedef std::size_t  Index;

    const Index   m_min = 10;
    const Index   n_min = 10;

    const Index   m_max = 2000;
    const Index   n_max = 2000;

    const Index   m_inc = 10;
    const Index   n_inc = 10;

    GeMatrix<T, Index>            A(m_max, n_max);
    GeMatrix<T, Index>            A_org(m_max, n_max);
    DenseVector<Index, Index>     p(m_max);

    hpc::aux::WallTime<double>  timer;

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

        randomInit(A);
        copy(A, A_org);
        //print(A_org, "A_org");

        timer.tic();
        /*
        Index info = hpc::ulmlapack::getf2(m, n,
                                           A.data, A.incRow, A.incCol,
                                           p.data, p.inc);
        */
        Index info = hpc::ulmlapack::tf2(A(0,0,m,n), p(0,m));
        double t = timer.toc();
        print(A_org(0,0,m,n), "A_org");
        print(A(0,0,m,n), "A");
        print(p(0,m), "p");

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

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





































































































Home