MKL-DSS-DSS-Error, reordering problem

[Naveen_P_]

I am trying to solve a sparse system of equations but I am getting this error "MKL-DSS-DSS-Error, reordering problem" without any more specifics on what the error is. This is what my code looks like:

#include "mkl_dss.h"

#include "mkl_types.h"
/*
** Define the array and rhs vectors
*/
//#define NROWS         5
//#define NCOLS         5
//#define NNONZEROS     9
//#define NRHS         1
#if defined(MKL_ILP64)
#define MKL_INT long long
#else
#define MKL_INT int
#endif
//static const MKL_INT nRows =     NROWS     ;
//static const MKL_INT nCols =     NCOLS     ;
//static const MKL_INT nNonZeros = NNONZEROS  ;
//static const MKL_INT nRhs =      NRHS     ;
//static _INTEGER_t rowIndex[NROWS+1] = { 1, 6, 7, 8, 9, 10 };
//static _INTEGER_t columns[NNONZEROS] = { 1, 2, 3, 4, 5, 2, 3, 4, 5 };
//static _DOUBLE_PRECISION_t values[NNONZEROS] = { 9, 1.5, 6, .75, 3, 0.5, 12, .625, 16 };
//static _DOUBLE_PRECISION_t rhs[NCOLS] = { 1, 2, 3, 4, 5 };

void mkl_dss(double* K, int* col_ind, int* row_ptr, double* F, int nn, int nnz) {

    MKL_INT nRows = 2*nn;
    MKL_INT nCols = 2*nn;
    MKL_INT nNonZeros = nnz;
    MKL_INT nRhs = 1;
    MKL_INT i;

    _INTEGER_t *rowIndex, *columns;
    _DOUBLE_PRECISION_t *values, *rhs, *solValues;

    rowIndex = new _INTEGER_t[nRows + 1];
    columns = new _INTEGER_t[nnz];
    values = new _DOUBLE_PRECISION_t[nnz];
    rhs = new _DOUBLE_PRECISION_t[nCols];
    solValues = new _DOUBLE_PRECISION_t[nCols];

    for(i = 0; i < nnz; i++){
        values = K;
        columns = col_ind;
    }

    for(i = 0; i < nRows + 1; i++)
        rowIndex = row_ptr;
        
    for(i = 0; i < nCols; i++){
        rhs = F;
        solValues = 0.0;
    }
    
    /* Allocate storage for the solver handle and the right-hand side. */
//    _DOUBLE_PRECISION_t solValues[NROWS];
    _MKL_DSS_HANDLE_t handle;
    _INTEGER_t error;
    _CHARACTER_t statIn[] = "determinant";
    _DOUBLE_PRECISION_t statOut[5];
    MKL_INT opt = MKL_DSS_DEFAULTS;
    MKL_INT sym = MKL_DSS_SYMMETRIC;
    MKL_INT type = MKL_DSS_POSITIVE_DEFINITE;
/* --------------------- */
/* Initialize the solver */
/* --------------------- */
    error = dss_create(handle, opt );
    if ( error != MKL_DSS_SUCCESS ) goto printError;
/* ------------------------------------------- */
/* Define the non-zero structure of the matrix */
/* ------------------------------------------- */
    error = dss_define_structure(
        handle, sym, rowIndex, nRows, nCols,
        columns, nNonZeros );
    if ( error != MKL_DSS_SUCCESS ) goto printError;
/* ------------------ */
/* Reorder the matrix */
/* ------------------ */
    error = dss_reorder( handle, opt, 0);
    if ( error != MKL_DSS_SUCCESS ) goto printError;
/* ------------------ */
/* Factor the matrix */
/* ------------------ */
    error = dss_factor_real( handle, type, values );
    if ( error != MKL_DSS_SUCCESS ) goto printError;
/* ------------------------ */
/* Get the solution vector */
/* ------------------------ */
    error = dss_solve_real( handle, opt, rhs, nRhs, solValues );
    if ( error != MKL_DSS_SUCCESS ) goto printError;
/* ------------------------ */
/* Get the determinant (not for a diagonal matrix) */
/*--------------------------*/
    if ( nRows < nNonZeros ) {
        error = dss_statistics(handle, opt, statIn, statOut);
        if ( error != MKL_DSS_SUCCESS ) goto printError;
/*-------------------------*/
/* print determinant */
/*-------------------------*/
        printf(" determinant power is %g \n", statOut[0]);
        printf(" determinant base is %g \n", statOut[1]);
        printf(" Determinant is %g \n", (pow(10.0,statOut[0]))*statOut[1]);
    }
/* -------------------------- */
/* Deallocate solver storage */
/* -------------------------- */
    error = dss_delete( handle, opt );
    if ( error != MKL_DSS_SUCCESS ) goto printError;
/* ---------------------- */
/* Print solution vector */
/* ---------------------- */
    printf(" Solution array: ");
    for(i = 0; i< nCols; i++)
        printf(" %g", solValues );
    printf("\n");
    
printError:
    printf("Solver returned error code %d\n", error);

    for(i = 0; i < nCols; i++)
        F = solValues;

    delete[] rowIndex;
    delete[] columns;
    delete[] values; 
    delete[] rhs; 
    delete[] solValues;
}

The same piece of code worked for the default example problem, so I know the sparse solver compiles and links correctly. But it does not work when I use my own set of arrays. I have a hunch that opt = MKL_DSS_DEFAULTS is not working as it should. Should I be using something different ? I have used Fortran style indexing (starting from 1) for the columns and row index vectors. I have also tried zero based indexing and set the opt = MKL_DSS_MSG_LVL_WARNING + MKL_DSS_TERM_LVL_ERROR+MKL_DSS_ZERO_BASED_INDEXING which also did not work.

[mecej4]

If you write program based on a set of assumptions, and then run it with a new set data which may not satisfy the assumptions, what can you expect?

The example problem in dss_sym_c.c is positive-definite, and this can be verified. Is your K matrix known to be positive definite? Did you correctly fill in the arrays of K in the CSR representation? If not, blindly changing options and index conventions in an attempt to fix undiagnosed errors is probably going to fail.

[Naveen_P_]

mecej4 wrote:

If you write program based on a set of assumptions, and then run it with a new set data which may not satisfy the assumptions, what can you expect?

The example problem in dss_sym_c.c is positive-definite, and this can be verified. Is your K matrix known to be positive definite? Did you correctly fill in the arrays of K in the CSR representation? If not, blindly changing options and index conventions in an attempt to fix undiagnosed errors is probably going to fail.

Yes I checked everything you mentioned before I posted this. My matrix is symmetric positive definite and I checked my CSR representation. I also noticed that the example problem is 1 - indexed which is why my row index and column index vectors are also 1 - indexed.

[mecej4]

Some MKL routines can work with 0-based and 1-based array indices. Some work only with 1-based indices, and these can be noted by the absence of a subroutine argument to convey the indexing option that is desired.

You should perhaps put together a "reproducer" and post it here. 

[Naveen_P_]

This is the 4x4 matrix I used as an example:

10000 4.69541 -4.69541 -4.69541              
    4.69541 10000   -4.69541 -4.69541             
    -4.69541 -4.69541 10000 4.69541          
    -4.69541  -4.69541 4.69541  4.69541

The Rhs is Transpose{0 0 0 1000}

The arrays I used as entry into mkl_dss are:

K                          Columns

10000                    1
4.69541                 2
-4.69541                3
-4.69541                4
4.69541                 1
10000                    2
-4.69541                3
-4.69541                4
-4.69541                1
-4.69541                2
10000                    3
4.69541                 4
-4.69541                1
-4.69541                2
4.69541                 3
4.69541                 4

Row Index:

1
5
9
13
17

Number of non - zeros = 16

And it's easy to verify that K is positive definite and symmetric and not singular.

[mecej4]

Since you did not show the code, #5 is not yet a "reproducer", but I can see a hint of the problem. The MKL Reference Manual says, in the section Storage Formats for the Direct Sparse Solvers:

For symmetric matrices, it is necessary to store only the upper triangular half of the matrix (upper triangular format) or the lower triangular half of the matrix (lower triangular format).

The Intel MKL direct sparse solvers use a row-major upper triangular storage format: the matrix is compressed row-by-row and for symmetric matrices only non-zero elements in the upper triangular half of the matrix are stored.

See, for instance, https://software.intel.com/en-us/node/471374#9874F00E-C273-4838-AA94-7882DE6D749A .

You seem to be entering the full matrix, instead -- this I gather from your setting nnz=16 instead of 10, and the values that you displayed for the index and value arrays.

Please retry, with only the elements in the upper triangle of the symmetric matrix entered.

[Naveen_P_]

It seems that what you are referring to is the problem. When I select the matrix type to be non - symmetric, it is generating the correct solution. Thanks.

[mecej4]

Whether you (i) select "symmetric" and provide the upper triangle, or (ii) select "unsymmetric" and provide the full matrix, you should expect the same solution. However, the latter choice may double the run time, so it is worth making sure to select the former if running speed is important (particularly if the equation set will be solved many times in a single run of your program).