Re: Re:Wrong results with DSS and the ATandT/pre2 matrix from Matrix Marke

cpmech · ‎09-05-2023

Hello,

I'm solving the A*x=rhs problem with DSS and the ATandT/pre2 matrix from Matrix Market. The right-hand side is set to unity. Unfortunately, DSS is computing incorrect results. It is worth noting that I've also run the same problem with MUMPS, and the results are correct.

The difficult matrix is the ATandT/pre2 from here: https://sparse.tamu.edu/ATandT/pre2

By the way, Intel DSS yields correct results with other matrices from Matrix Market (e.g., the bfwb62.mtx, included in the zip file). Also, Valgrind reports no leaks.

I'm wondering if there is an option for DSS to warn about such problematic matrices.

The complete code is listed below and included in the zip file. The main function starts on line 216 with most of the code being the reading of the MatrixMarket file and converting the COO (coordinates) matrix format to CSR (compressed sparse row) format.

#include "mkl.h"
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <memory>
#include <string>
#include <vector>

struct CooMatrix {
    MKL_INT m;
    MKL_INT nnz;
    MKL_INT nnz_max;
    bool symmetric;
    std::vector<MKL_INT> indices_i;
    std::vector<MKL_INT> indices_j;
    std::vector<double> values_aij;

    inline static std::unique_ptr<CooMatrix> make_new(MKL_INT m, MKL_INT nnz_max) {
        return std::unique_ptr<CooMatrix>{new CooMatrix{
            m,
            0,
            nnz_max,
            false,
            std::vector<MKL_INT>(nnz_max, 0),
            std::vector<MKL_INT>(nnz_max, 0),
            std::vector<double>(nnz_max, 0.0),
        }};
    }

    void put(MKL_INT i, MKL_INT j, double aij) {
        if (i < 0 || i >= static_cast<MKL_INT>(this->m)) {
            throw "CooMatrix::put: index of row is outside range";
        }
        if (j < 0 || j >= static_cast<MKL_INT>(this->m)) {
            throw "CooMatrix::put: index of column is outside range";
        }
        if (this->nnz >= this->nnz_max) {
            throw "CooMatrix::put: max number of items has been exceeded";
        }
        this->indices_i[this->nnz] = i;
        this->indices_j[this->nnz] = j;
        this->values_aij[this->nnz] = aij;
        this->nnz++;
    }
};

std::unique_ptr<CooMatrix> read_matrix_market(const std::string &filename) {
    FILE *f = fopen(filename.c_str(), "r");
    if (f == NULL) {
        std::cout << "filename = '" << filename << "'" << std::endl;
        throw "read_matrix_market: cannot open file";
    }

    const int line_max = 2048;
    char line[line_max];

    if (fgets(line, line_max, f) == NULL) {
        fclose(f);
        throw "read_matrix_market: cannot read any line in the file";
    }

    char mm[24], opt[24], fmt[24], kind[24], sym[24];
    int nread = sscanf(line, "%24s %24s %24s %24s %24s", mm, opt, fmt, kind, sym);
    if (nread != 5) {
        fclose(f);
        throw "read_matrix_market: number of tokens in the header is incorrect";
    }
    if (strncmp(mm, "%%MatrixMarket", 14) != 0) {
        fclose(f);
        throw "read_matrix_market: header must start with %%MatrixMarket";
    }
    if (strncmp(opt, "matrix", 6) != 0) {
        fclose(f);
        throw "read_matrix_market: option must be \"matrix\"";
    }
    if (strncmp(fmt, "coordinate", 10) != 0) {
        fclose(f);
        throw "read_matrix_market: type must be \"coordinate\"";
    }
    if (strncmp(kind, "real", 4) != 0) {
        fclose(f);
        throw "read_matrix_market: number kind must be \"real\"";
    }
    if (strncmp(sym, "general", 7) != 0 && strncmp(sym, "symmetric", 9) != 0) {
        fclose(f);
        throw "read_matrix_market: matrix must be \"general\" or \"symmetric\"";
    }

    std::unique_ptr<CooMatrix> coo;
    bool symmetric = strncmp(sym, "symmetric", 9) == 0;

    bool initialized = false;
    size_t m, n, nnz, i, j;
    double x;

    while (fgets(line, line_max, f) != NULL) {
        if (!initialized) {
            if (line[0] == '%') {
                continue;
            }
            nread = sscanf(line, "%zu %zu %zu", &m, &n, &nnz);
            if (nread != 3) {
                fclose(f);
                throw "read_matrix_market: cannot parse the dimensions (m,n,nnz)";
            }
            coo = CooMatrix::make_new(m, nnz);
            coo->symmetric = symmetric;
            initialized = true;
        } else {
            nread = sscanf(line, "%zu %zu %lg", &i, &j, &x);
            if (nread != 3) {
                fclose(f);
                throw "read_matrix_market: cannot parse the values (i,j,x)";
            }
            if (symmetric) {
                // must swap lower triangle to upper triangle
                // because MatrixMarket is lower triangle and
                // intel DSS requires upper triangle
                coo->put(j - 1, i - 1, x);
            } else {
                coo->put(i - 1, j - 1, x);
            }
        }
    }

    fclose(f);
    return coo;
}

struct CsrMatrix {
    MKL_INT m;
    MKL_INT nnz;
    MKL_INT *row_pointers = NULL;
    MKL_INT *column_indices = NULL;
    double *values = NULL;
    sparse_matrix_t handle = NULL;

    static std::unique_ptr<CsrMatrix> from(std::unique_ptr<CooMatrix> &coo) {
        // access triplet data
        auto ai = coo->indices_i.data();
        auto aj = coo->indices_j.data();
        auto ax = coo->values_aij.data();

        // create COO MKL matrix
        auto m = coo->m;
        auto nnz = coo->nnz;
        sparse_matrix_t coo_mkl = NULL;
        auto status = mkl_sparse_d_create_coo(&coo_mkl, SPARSE_INDEX_BASE_ZERO, m, m, nnz, ai, aj, ax);
        if (status != SPARSE_STATUS_SUCCESS) {
            throw "Intel MKL failed to create COO matrix";
        }

        // convert COO to CSR
        sparse_matrix_t csr = NULL;
        status = mkl_sparse_convert_csr(coo_mkl, SPARSE_OPERATION_NON_TRANSPOSE, &csr);
        if (status != SPARSE_STATUS_SUCCESS) {
            throw "Intel MKL failed to convert COO matrix to CSR matrix";
        }

        // destroy handle to COO
        mkl_sparse_destroy(coo_mkl);

        // get access to internal arrays
        sparse_index_base_t indexing;
        MKL_INT *pointer_b = NULL; // size = dimension
        MKL_INT *pointer_e = NULL; // size = dimension, last entry == nnz (when zero-based)
        MKL_INT *columns = NULL;   // size = nnz
        double *values = NULL;     // size = nnz

        // NOTE about mkl_sparse_d_export_csr
        // The routine returns pointers to the internal representation
        // and DOES NOT ALLOCATE additional memory.
        // https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-2/mkl-sparse-export-csr.html

        status = mkl_sparse_d_export_csr(csr, &indexing, &m, &m, &pointer_b, &pointer_e, &columns, &values);
        if (status != SPARSE_STATUS_SUCCESS) {
            throw "Intel MKL failed to export csr matrix";
        }

        // compute the row indices vector because we need a CSR3 matrix
        MKL_INT *row_pointers = (MKL_INT *)malloc((m + 1) * sizeof(MKL_INT));
        if (row_pointers == NULL) {
            throw "cannot allocate row_pointers";
        }
        for (MKL_INT i = 0; i < m; i++) {
            row_pointers[i] = pointer_b[i];
        }
        row_pointers[m] = pointer_e[m - 1];

        // final number of non-zero values
        MKL_INT final_nnz = row_pointers[m];

        // results
        return std::unique_ptr<CsrMatrix>{new CsrMatrix{
            m,
            final_nnz,
            row_pointers,
            columns,
            values,
            csr,
        }};
    }

    ~CsrMatrix() {
        if (this->row_pointers != NULL) {
            free(this->row_pointers);
        }
        if (this->handle != NULL) {
            mkl_sparse_destroy(this->handle);
        }
    }
};

int main(int argc, char **argv) try {
    // read the matrix
    auto home = std::string(std::getenv("HOME"));
    // auto matrix = home + "/Downloads/matrix-market/bfwb62.mtx";
    auto matrix = home + "/Downloads/matrix-market/pre2.mtx";

    // convert to CSR
    auto coo = read_matrix_market(matrix);
    auto csr = CsrMatrix::from(coo);

    // set some constants for DSS
    MKL_INT opt = MKL_DSS_MSG_LVL_WARNING + MKL_DSS_TERM_LVL_ERROR + MKL_DSS_ZERO_BASED_INDEXING;
    MKL_INT sym = MKL_DSS_NON_SYMMETRIC;
    if (coo->symmetric) {
        sym = MKL_DSS_SYMMETRIC;
    }
    MKL_INT type = MKL_DSS_INDEFINITE;

    // initialize the solver
    _MKL_DSS_HANDLE_t handle;
    auto error = dss_create(handle, opt);
    if (error != MKL_DSS_SUCCESS) {
        throw "DSS failed when creating the solver handle";
    }

    // define the non-zero structure of the matrix
    error = dss_define_structure(handle, sym, csr->row_pointers, csr->m, csr->m, csr->column_indices, csr->nnz);
    if (error != MKL_DSS_SUCCESS) {
        throw "DSS failed when defining the non-zero structure of the matrix";
    }

    // reorder the matrix
    error = dss_reorder(handle, opt, 0);
    if (error != MKL_DSS_SUCCESS) {
        throw "DSS failed when reordering the matrix";
    }

    // factor the matrix
    error = dss_factor_real(handle, type, csr->values);
    if (error != MKL_DSS_SUCCESS) {
        throw "DSS failed when factorizing the matrix";
    }

    // get the solution vector
    auto rhs = std::vector<double>(csr->m, 1.0);
    auto x = std::vector<double>(csr->m, 0.0);
    MKL_INT n_rhs = 1;
    error = dss_solve_real(handle, opt, rhs.data(), n_rhs, x.data());
    if (error != MKL_DSS_SUCCESS) {
        throw "DSS failed when solving the system";
    }

    // clear dss data
    dss_delete(handle, opt);

    // check the solution
    auto op = SPARSE_OPERATION_NON_TRANSPOSE;
    auto desc = matrix_descr{
        coo->symmetric ? SPARSE_MATRIX_TYPE_SYMMETRIC : SPARSE_MATRIX_TYPE_GENERAL,
        SPARSE_FILL_MODE_UPPER,
        SPARSE_DIAG_NON_UNIT,
    };
    auto rhs_new = std::vector<double>(csr->m, 0.0);
    mkl_sparse_d_mv(op, 1.0, csr->handle, desc, x.data(), 0.0, rhs_new.data());
    for (MKL_INT k = 0; k < coo->m; k++) {
        auto diff = fabs(rhs[k] - rhs_new[k]);
        if (diff > 1e-10) {
            std::cout << "diff[" << k << "] = " << diff << " is too high" << std::endl;
            return 1;
        }
    }
    return 0;

} catch (std::exception &e) {
    std::cout << e.what() << std::endl;
} catch (std::string &msg) {
    std::cout << msg.c_str() << std::endl;
} catch (const char *msg) {
    std::cout << msg << std::endl;
} catch (...) {
    std::cout << "some unknown exception happened" << std::endl;
}

Any assistance is appreciated.

Kind regards.

Doriv4.

ShanmukhS_Intel · ‎09-08-2023

Hi Dorival,

Thanks for posting in Intel Communities.

Could you please get back to us with the MKL version and environmental details?

We have tried running the shared files. However, we are facing an error as pre2.mtx seems to be missing among the shared files. Could you please share the same?

Best Regards,

Shanmukh.SS

cpmech · ‎09-08-2023

Hi Shanmukh,

Thanks for looking into this.

I've used the default oneAPI MKL version that comes with Ubuntu 22.04 and Ubuntu 23.04 (I've tried them both). Specifically, the following version: 2023.2.0.

My system is based on a 13th Gen Intel(R) Core(TM) i9-13900KF

The pre2.mtx file is not included in the zip file because it's too large (129 Mb). However, it can be easily downloaded from https://suitesparse-collection-website.herokuapp.com/ATandT/pre2 (select Matrix Market version at the bottom of the page)

I have also uploaded my project to GitHub (see here https://github.com/cpmech/intel-matrix-market), where you may use a script to download the matrix market files (in the scripts folder).

Alternatively, you may use Docker. I've built the project into a Docker image that can be easily downloaded with:

docker pull cpmech/intel-matrix-market

You may "shell into Docker" with the following command:

docker run --rm -it cpmech/intel-matrix-market:latest /bin/bash

Then, download the Matrix Market files:

bash scripts/download-from-matrix-market.bash

And run the code:

bash all.bash

Which will output the following:

-- The C compiler identification is GNU 11.4.0
-- The CXX compiler identification is GNU 11.4.0
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Check for working C compiler: /usr/bin/cc - skipped
-- Detecting C compile features
-- Detecting C compile features - done
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Check for working CXX compiler: /usr/bin/c++ - skipped
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- Upgrade to CMake version 3.25.2 or later for native Intel(R) oneAPI DPC++/C++ Compiler support. You are running version 3.22.1.
-- MKL_ARCH: None, set to ` intel64` by default
-- MKL_ROOT /opt/intel/oneapi/mkl/2023.2.0
-- MKL_LINK: None, set to ` dynamic` by default
-- MKL_INTERFACE_FULL: None, set to ` intel_ilp64` by default
-- MKL_THREADING: None, set to ` intel_thread` by default
-- Found MKL: /opt/intel/oneapi/mkl/2023.2.0  
-- Found MKL: /opt/intel/oneapi/mkl/2023.2.0/lib/intel64/libmkl_intel_ilp64.so  
-- Found MKL: /opt/intel/oneapi/mkl/2023.2.0/lib/intel64/libmkl_intel_thread.so  
-- Found MKL: /opt/intel/oneapi/mkl/2023.2.0/lib/intel64/libmkl_core.so  
-- Found MKL: /opt/intel/oneapi/compiler/latest/linux/compiler/lib/intel64_lin/libiomp5.so  
-- Configuring done
-- Generating done
-- Build files have been written to: /tmp/build-intel-matrix-market
[ 50%] Building CXX object CMakeFiles/solve_matrix_market.dir/solve_matrix_market.cpp.o
[100%] Linking CXX executable solve_matrix_market
[100%] Built target solve_matrix_market
diff[0] = 5.43927e+67 is too high

Note the last line: diff[0] = 5.43927e+67 is too high

Thanks again.

Dorival

ShanmukhS_Intel · ‎09-17-2023

Hi,

Thanks for sharing the details. We are looking into your issue. We will get back to you soon with an update.

Best Regards,

Shanmukh.SS

ShanmukhS_Intel · ‎09-25-2023

Hi,

Thanks for the details.

We are able to reproduce your issue, we are working on your issue. We will get back to you soon.

Thanks & Regards,

Shanmukh.SS

cpmech · ‎09-29-2023

Hello,

Thanks again for looking into this.

For reference, I've got more information about the ATandT/pre2 matrix (calculated by MUMPS):

normA = 2.95543e+10
normX = 1.05307e+13
resid = 3.96148e-28
omega1 = 1.07564e-16
omega2 = 4.81081e-21
delta = 1.21214e-08
cond1 = 2
cond2 = 2.51961e+12

where omega1 and omega2 are the backward error estimates according to the Arioli, Demmel, Duff 1989 paper and according to MUMPS' manual page 14. cond1 is the condition number for the matrix, and cond2 is the condition number for the linear system. delta is the ratio between small changes in the solution vector norm, divided by the solution vector norm (MUMPS manual page 15).

Yes, it's a "bad" matrix and corresponding linear system :-D. However, MUMPS solves the problem (see log below):

{
  "platform": "Russell",
  "blas_lib": "OpenBLAS",
  "solver_name": "MUMPS-local",
  "matrix_name": "pre2",
  "symmetry": "None",
  "nrow": 659033,
  "ncol": 659033,
  "nnz": 5959282,
  "time_read_matrix_market_nanosecond": 475828606,
  "time_read_matrix_market_human": "475.828606ms",
  "time_factorize_nanosecond": 7037532439,
  "time_factorize_human": "7.037532439s",
  "time_solve_nanosecond": 972491236,
  "time_solve_human": "972.491236ms",
  "time_total_nanosecond": 8010023675,
  "time_total_human": "8.010023675s",
  "requested_ordering": "Auto",
  "requested_scaling": "Auto",
  "requested_openmp_num_threads": 1,
  "effective_ordering": "Metis",
  "effective_scaling": "Unknown",
  "effective_strategy": "Unknown",
  "effective_openmp_num_threads": 32,
  "verify_max_abs_a": 14777170000.0,
  "verify_max_abs_a_times_x": 1.0000829040072858,
  "verify_relative_error": 8.34342487157907e-15,
  "verify_time_nanosecond": 7836414,
  "verify_time_human": "7.836414ms",
  "compute_determinant": false,
  "determinant_mantissa": 0.0,
  "determinant_base": 2.0,
  "determinant_exponent": 0.0
}

Below is the output from my code with Intel DSS:

{
  "platform": "Russell",
  "blas_lib": "OpenBLAS",
  "solver_name": "IntelDSS",
  "matrix_name": "pre2",
  "symmetry": "None",
  "nrow": 659033,
  "ncol": 659033,
  "nnz": 5959282,
  "time_read_matrix_market_nanosecond": 470561468,
  "time_read_matrix_market_human": "470.561468ms",
  "time_factorize_nanosecond": 5338142985,
  "time_factorize_human": "5.338142985s",
  "time_solve_nanosecond": 126911638,
  "time_solve_human": "126.911638ms",
  "time_total_nanosecond": 5465054623,
  "time_total_human": "5.465054623s",
  "requested_ordering": "Auto",
  "requested_scaling": "Auto",
  "requested_openmp_num_threads": 1,
  "effective_ordering": "Unknown",
  "effective_scaling": "Unknown",
  "effective_strategy": "Unknown",
  "effective_openmp_num_threads": 32,
  "verify_max_abs_a": 14777170000.0,
  "verify_max_abs_a_times_x": 1.6086144329336344e122,
  "verify_relative_error": 1.0885808533195302e112,
  "verify_time_nanosecond": 8562183,
  "verify_time_human": "8.562183ms",
  "compute_determinant": false,
  "determinant_mantissa": 0.0,
  "determinant_base": 10.0,
  "determinant_exponent": 0.0
}

In particular, note the line: "verify_relative_error": 1.0885808533195302e112

For reference, below is the header of the Matrix Market file:

%%MatrixMarket matrix coordinate real general
%-------------------------------------------------------------------------------
% UF Sparse Matrix Collection, Tim Davis
% http://www.cise.ufl.edu/research/sparse/matrices/ATandT/pre2
% name: ATandT/pre2
% [AT&T,harmonic balance method, large example]
% id: 285
% date: 2001
% author: R. Melville, D. Long
% ed: T. Davis
% fields: title A Zeros name id date author ed kind
% kind: frequency-domain circuit simulation problem
%-------------------------------------------------------------------------------

This matrix gives me nightmares LOL.

Cheers.

Dorival

ShanmukhS_Intel · ‎10-09-2023

Hi,

Regarding any option in DSS to receive warnings or notifications about problematic matrices.
>> Unfortunately, It is not possible to receive warnings when using the DSS interface.

Could you please try using the MKL Pardiso Interface instead of DSS one?

Best Regards,
Shanmukh.SS

ShanmukhS_Intel · ‎10-16-2023

Hi,

A gentle reminder:

Has the information provided helped? Could you please let us know if you need any other information?

Best Regards,

Shanmukh.SS

ShanmukhS_Intel · ‎11-13-2023

Hi,

A gentle reminder:

Has the information provided helped? Could you please let us know if you need any other information?

Best Regards,

Shanmukh.SS