BLR compression solver gives wrong result

[dimitargslavchev]

Hello PIeter,

I am trying to solve a system of linear algebraic equations with the BLR solver. I am using testBLR_seq as an example. I have successfully solved the same problem with the HSS solver with some accuracy. The matrix is not symmetric.

However the results that I receive are very different from the results obtained from a direct LU solver (dgesv from MKL). In the MKL result all values are ~10^{1} or below, while BLR gives me results that are between ~10^{6} to ~10^{12}.

I tried changing the rel_tol to a smaller value of 10^{-8} or even set all tiles to be non-compressed (i.e. adm.fill(false);). I always receive the wrong answer.

Here is the code of the solver part that I use:

// STRUMPACK headers
#include "dense/DenseMatrix.hpp"
#include "BLR/BLRMatrix.hpp"

using namespace std;
using namespace strumpack;
using namespace strumpack::BLR;

// my_size is matrix size; D is the matrix (col major order) and sv is the right hand side arrays
void my_solve(int my_size, double *D, double *sv)
{
    //////////////////////////////////////////////////////////
    // S T R U M P A C K    W O R K S
    //    HSSOptions<double> blr_opts;
    BLROptions<double> blr_opts;
    blr_opts.set_verbose(false);
    blr_opts.set_max_rank(my_size);

    blr_opts.set_rel_tol(1e-2);
    //blr_opts.set_abs_tol(1e-8);

    // OR manually define the tree
    strumpack::structured::ClusterTree tree(my_size);
    tree.refine(blr_opts.leaf_size());

    // Get a vector of tiles
    auto tiles=tree.template leaf_sizes<std::size_t>();

    //ADMISSIBILITY -- weak
    std::size_t nt = tiles.size();
    DenseMatrix<bool> adm(nt, nt);
    adm.fill(true);
    for (std::size_t t=0; t<nt; t++) {
        adm(t, t) = false;
    }
    std::vector<int> piv;

    // Initialize the matrices
    // rows = N, cols = n, Array = D, Lead Dimention(== cols) = N
    DenseMatrix<double> A = DenseMatrix<double>(my_size,my_size,D,my_size);
    DenseMatrix<double> B = DenseMatrix<double>(my_size,1,sv,my_size);

    // Compress the matrix into BLR. directly from constructor
    BLRMatrix<double> H(A, tiles, adm, piv, blr_opts);

    H.solve(piv, B);
    copy(B, sv, my_size);

}

This is the experiment matrix, right hand side and the result from MKL and BLR with $rel_tol = 10^{-2}$ BLR_test_1000.zip

[pghysels]

Dear Dimitar,

My apologies for the delay. I tried your code, with the matrix you sent. And indeed the results seem wrong.

I believe the problem is a numerical issue. MKL dgesv performs partial pivoting, it looks for the largest element in each column, and uses that element as the pivot. But in BLR, we only look for the largest element in the column within the diagonal block. I think for this matrix, this restricted pivoting might not be sufficient. We could implement another variant of the BLR compression and factorization algorithm that is more numerically stable, I will discuss this with my collaborators. Another thing you could try is to add a diagonal shift to the matrix, then compress and factor the shifted matrix and use that as a preconditioner.

We will try to fix the StructuredMatrix interface soon.

Pieter

[dimitargslavchev]

Thank you for the quick response.

Best Regards, Dimitar Slavchev