lu_unblk-test.cc (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

#include <cxxstd/algorithm.h>

#include <flens/flens.cxx>

#include <flens/examples/lu_unblk.h>

using namespace flens;

using namespace std;

int

main()

{

    typedef GeMatrix<FullStorage<double> >   RealGeMatrix;

    const int m = 2000;

    const int n = 2000;

    const int mn = std::min(m, n);

///

/// Setup test case and make a copy of the input

///

    RealGeMatrix  A(m, n), A_check;

    fillRandom(A);

    A_check = A;

///

///  Factorize $A$.

///

    if (lu_unblk(A)) {

        cerr << "A is singular" << endl;

        return 1;

    }

///

/// Compute $|A - L*U|_1$ by brute force.  Note that $A$ is not required to

/// be square.

///

    Underscore<int> _;

    RealGeMatrix L = A(_,_(1,mn)).lowerUnit();

    RealGeMatrix U = A(_(1,mn),_).upper();

    A_check -= L*U;

    cout << "asum(A-L*U) = " << blas::asum(A_check) << endl;

    return 0;

}