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

#include <cassert>
#include <random>
#include <hpc/matvec/copy.h>
#include <hpc/matvec/gematrix.h>
#include <hpc/matvec/mm.h>
#include <hpc/matvec/print.h>
#include <hpc/ulmblas/trlsm.h>
#include <hpc/ulmblas/trusm.h>
#include <hpc/ulmlapack/getf2.h>
#include <hpc/ulmlapack/swap.h>


//
//  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()
{
    using namespace hpc::matvec;

    typedef double       T;
    typedef std::size_t  Index;

    const Index   m = 600;
    const Index   n = 800;

    GeMatrix<T>      A(m, m);
    GeMatrix<T>      X(m, n);
    GeMatrix<T>      B(m, n);
    GeMatrix<Index>  p(m, 1);

    randomInit(m, m, A.data, A.incRow, A.incCol);
    randomInit(m, n, X.data, X.incRow, X.incCol);

    //
    //  A*X = B
    //
    mm(T(1), A, X, T(0), B);

    //
    //  Factorize A=P*L*U
    //
    std::size_t info = hpc::ulmlapack::getf2(A.numRows, A.numCols,
                                             A.data, A.incRow, A.incCol,
                                             p.data, A.incRow);
    std::printf("getf2 returned: info = %ld\n", info);

    //
    //  Solve P*L*U*X = B
    //
    hpc::ulmlapack::swap(B.numRows, B.numCols, B.data, B.incRow, B.incCol,
                         Index(0), B.numRows,
                         p.data, p.incRow);
    hpc::ulmblas::trlsm(m, n, T(1), true,
                        A.data, A.incRow, A.incCol,
                        B.data, B.incRow, B.incCol);
    hpc::ulmblas::trusm(m, n, T(1), false,
                        A.data, A.incRow, A.incCol,
                        B.data, B.incRow, B.incCol);

    //
    //  Compute residual
    //
    double res = 0;
    for (Index j=0; j<n; ++j) {
        for (Index i=0; i<m; ++i) {
            res += std::abs(B(i,j)-X(i,j));
        }
    }

    std::printf("res=%e\n", res);
}





































































































Home