Hi!

I wonder if anyone can help with or verify real problem with the following, which has me stumped and questioning my sanity. Question is about what could be (hope not)an error in the Pardiso and DSS forward solve routine, seen e.g. in the following problem:

Consider system Ax=b, where the numerical values are

3 1 0 1

1 2 2 x = 2

0 2 3 3

This system has solution (0.66667, -1, 1.66667), which is correctly found by Pardiso and DSS.

However, something funny appears to happen with the forward solve mode: According to the MKL 10.2 manual, Pardiso and DSS do for sparse symmetric positive definite systems an A=LL^T Cholesky factorization. The correct solution for L (verified by Mathematica, and which can be verified manually as well), is

sqrt(3) 0 0

L = 1/sqrt(3) sqrt(5/3) 0

0 2*sqrt(3/5) sqrt(3/5)

The pardiso and DSS step 331 forward substitution should solve Ly=b. Fromabove one sees by substitution that sqrt(3)x1=1 => x1 = 1/sqrt(3), and so on, giving solution for y, y=(0.57735, 1.29099, 1.29099).

Instead of giving this result, Pardiso and DSS both give (0.57735, -0.57735, 1.73205).

From several things we tried, I should note the following:

- Iterative refinement was turned off

-IF the upper right hand corner element is entered as zero (logically present making it a full upper triangle, but with one element with value zero), both pardiso and DSS return the correct result. This strongly suggests that DSS operates wrong with the sparse structure.

- Result did not seem to depend on whether full solve has been done before fwd solve or not.

- It occurs also with Pardiso default parameters

- It occurred for another researcher who independently implemented the problem using DSS, suggesting that at the very least documentation is misleading

- The same error occurred with larger matrixes as well, which first prompted us to try verify operation with the small matrix above. (The full solve works fine.)

- The computer is pretty typical 2-processor Dell laptop, running 64bit Windows 7, MKL 10.2, and the Intel C++ compiler under QT creator.

I have pasted the relevant bit of code below. It is effectively the symmetric positive definite matrix sample code from ref manual section C, with the matrix replaced by the matrix above, and solve put to fwd solve without iterative refinement.

A solution to this problem or even verification of error condition would be greatly appreciated. We really need both forward and backward solutions working in our application, and it would be great to be able to proceed. We cannot work with dense Cholesky as the matrix dimension may go up to millions.

Additionally, I'd add a vote to earlier mails asking for Cholesky factors being made visible from Pardiso/DSS. The decomposition has several uses in computation. At the very least, Pardiso could report the determinant, which can be computed from the diagonal of the L matrix (DSS does this fine, = 3 for the matrix above).

Kind regards - thanks for any help - and sorry for the slightly wordy description above!

Atte Moilanen, University of Helsinki, Finland

------------------ the relevant function -------------------------

#include

#include

#include

#include

#include "mkl_pardiso.h"

#include "mkl_types.h"

#include "mkl_dss.h"

/***************************************************************/

int MKL_DSS_tst() {

#define NROWS 3

#define NCOLS 3

#define NNONZEROS 5

#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, 3, 5, 6};

static _INTEGER_t columns[NNONZEROS] = { 1, 2, 2, 3, 3 };

static _DOUBLE_PRECISION_t values[NNONZEROS] = { 3,1,2,2,3 };

static _DOUBLE_PRECISION_t rhs[NCOLS] = { 1, 2, 3};

MKL_INT i;

/* 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_STR_t statIn[] = "determinant";

_DOUBLE_PRECISION_t statOut[5];

MKL_INT opt = MKL_DSS_DEFAULTS;

MKL_INT sol_opt = MKL_DSS_FORWARD_SOLVE+MKL_DSS_REFINEMENT_OFF;

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, sol_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[i] );

printf("\\n");

exit(0);

printError:

printf("Solver returned error code %d\\n", error);

exit(1);

}