lu_blk_test.cc (flens/examples/lu_blk

      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

#include <flens/flens.cxx>

#include <iostream>

#ifdef WITH_OPERATORS

#   include <flens/examples/lu/lu_blk_with_operators.h>

#else

#   include <flens/examples/lu/lu_blk_with_flensblas.h>

#endif

#include <flens/examples/lu/apply_perm.h>

using namespace flens;

using namespace std;

int

main()

{

    const int m = 4000;

    const int n = 4000;

    GeMatrix<FullStorage<double, RowMajor> >  A(m,n), A_org;

    DenseVector<Array<int> >        p(m);

    Underscore<int>                 _;

    fillRandom(A);

    // store original A for checking the factorization later

    A_org = A;

    /// Uncomment for playing with small dimensions:

    //cout << "A = " << A << endl;

    if (lu_blk(A, p) != 0) {

        cout << "Matrix is singular or close to singular" << endl;

        return 1;

    }

#   ifdef CHECK

    /// and

    //cout << "L = " << A.lowerUnit() << endl;

    //cout << "U = " << A.upper() << endl;

    //cout << "p = " << p << endl;

    GeMatrix<FullStorage<double> >  LU(m,n);

    // blas::copy(NoTrans, A.upper(), LU.upper());

    LU.upper() = A.upper();

    // blas::mm(Left, NoTrans, 1.0, A(_(1,4),_(1,4)).lowerUnit(), LU);

    LU = A(_(1,m),_(1,m)).lowerUnit() * LU;

    // LU = P*LU

    apply_perm(p, LU);

    //cout << "P*L*U = " << LU << endl;

    // blas::axpy(NoTrans, -1.0, LU, A_org);

    A_org -= LU;

    //cout << "A - P*L*U = " << A_org << endl;

    cout << "|| A-P*L*U ||_1 = " << blas::asum(A_org) << endl;

    return 0;

#   endif

}