Problems solving sparse complex non-symmetric matrix with Padiso DSS interface routines

[mdajunge]

I encoutered some problems, when trying to solve a sparse non-symmetrix complex matrix with pardiso and the DSS interface routines. I receive a segmentation fault during the call "dss_factor_complex( handle, type, values );"

To reproduce the problem, I adapted the example in
/opt/intel/mkl/10.1.1.019/examples/solver/source/dss_sym_c.c
to solve the system with the (8x8) sparse complex non symmetric matrix which is provided in the example file
/opt/intel/mkl/10.1.1.019/examples/solver/source/pardiso_unsym_complex_c.c

I compiled and linked as specified in the makefile /opt2/intel/mkl/10.1.1.019/examples/solver/makefile
with gcc (GCC) 4.1.2 20061115 (prerelease) (SUSE Linux)
The following libraries are listed when calling ldd
/pardiso_test>ldd ./dss_usym_complex_c.out
libiomp5.so => /opt2/intel/mkl/10.1.1.019/lib/em64t/libiomp5.so (0x00002b61ebf46000)
libpthread.so.0 => /lib64/libpthread.so.0 (0x00002b61ec146000)
libm.so.6 => /lib64/libm.so.6 (0x00002b61ec361000)
libc.so.6 => /lib64/libc.so.6 (0x00002b61ec5b8000)
/lib64/ld-linux-x86-64.so.2 (0x00002b61ebf28000)
libdl.so.2 => /lib64/libdl.so.2 (0x00002b61ec8fa000)

If I modify the matrix to be real-valued and non-symmetric or complex and symmetric I get the correct results without any segmentation fault.

Below I attached the definition of matrix and the rhs vector. If authorized by Intel, I could also post the C-file, which is a modified version of the example /opt/intel/mkl/10.1.1.019/examples/solver/source/dss_sym_c.c

#define NROWS 8
#define NCOLS 8
#define NNONZEROS 20
#define NRHS 1
static _INTEGER_t rowIndex[NROWS+1] = { 1, 5, 8, 10, 12, 13, 16, 18, 21 };
static _INTEGER_t columns[NNONZEROS] = { 1,3, 6, 7, 2, 3, 5, 3, 8,4,7, 2, 3,6,8,2,7,3,7,8};
static _DOUBLE_COMPLEX_t values[NNONZEROS];
static _DOUBLE_COMPLEX_t rhs[NROWS];
in main:
values[0 ].r = 7.0, values[0 ].i = 1.0;
values[1 ].r = 1.0, values[1 ].i = 1.0;
values[2 ].r = 2.0, values[2 ].i = 1.0;
values[3 ].r = 7.0, values[3 ].i = 1.0;
values[4 ].r = -4.0, values[4 ].i = 0.0;
values[5 ].r = 8.0, values[5 ].i = 1.0;
values[6 ].r = 2.0, values[6 ].i = 1.0;
values[7 ].r = 1.0, values[7 ].i = 1.0;
values[8 ].r = 5.0, values[8 ].i = 1.0;
values[9 ].r = 7.0, values[9 ].i = 0.0;
values[10].r = 9.0, values[10].i = 1.0;
values[11].r = -4.0, values[11].i = 1.0;
values[12].r = 7.0, values[12].i = 1.0;
values[13].r = 3.0, values[13].i = 1.0;
values[14].r = 8.0, values[14].i = 0.0;
values[15].r = 1.0, values[15].i = 1.0;
values[16].r = 11.0, values[16].i = 1.0;
values[17].r = -3.0, values[17].i = 1.0;
values[18].r = 2.0, values[18].i = 1.0;
values[19].r = 5.0, values[19].i = 0.0;
rhs[0 ].r = 1.0; rhs[0 ].i = 1.0;
rhs[1 ].r = 1.0; rhs[1 ].i = 1.0;
rhs[2 ].r = 1.0; rhs[2 ].i = 1.0;
rhs[3 ].r = 1.0; rhs[3 ].i = 1.0;
rhs[4 ].r = 1.0; rhs[4 ].i = 1.0;
rhs[5 ].r = 1.0; rhs[5 ].i = 1.0;
rhs[6 ].r = 1.0; rhs[6 ].i = 1.0;
rhs[7 ].r = 1.0; rhs[7 ].i = 1.0;

Best regards
Michael Junge

[Gennady_F_Intel]

Is it 32 or 64 bit OS?
Which an additional options else did you use?
I mean for example: make < libem64t | dllem64t >
SUSE's version?
--Gennady

[mdajunge]

To compile the provided examples I used
make libem64t compiler=gnu

For my own demo program I took the command generate by make and adapted it. I then compiled from command line with:

gcc -w -I/opt2/intel/mkl/10.1.1.019/include dss_usym_complex_c.c -L"/opt2/intel/mkl/10.1.1.019/lib/em64t" "/opt2/intel/mkl/10.1.1.019/lib/em64t"/libmkl_solver_lp64.a "/opt2/intel/mkl/10.1.1.019/lib/em64t"/libmkl_intel_lp64.a -Wl,--start-group "/opt2/intel/mkl/10.1.1.019/lib/em64t"/libmkl_intel_thread.a "/opt2/intel/mkl/10.1.1.019/lib/em64t"/libmkl_core.a -Wl,--end-group -L"/opt2/intel/mkl/10.1.1.019/lib/em64t" -liomp5 -lpthread -lm -o dss_usym_complex_c.out

I tested both the provided examples and the my demo on two machines:
1)
--Intel Xeon CPU 5160 @ 3.00GHz,
--SLES 10 SP1
--Linux 2.6.16.46-0.12-smp #1 SMP Thu May 17 14:00:09 UTC 2007 x86_64 x86_64 x86_64 GNU/Linux
--gcc (GCC) 4.1.2 20070115 (prerelease) (SUSE Linux)
2)
--Intel Core2 CPU 6300
--Linux 2.6.18.8-0.7-default #1 SMP Tue Oct 2 17:21:08 UTC 2007 x86_64 x86_64 x86_64 GNU/Linux
--openSUSE 10.2 (X86-64)
--gcc (GCC) 4.1.2 20061115 (prerelease) (SUSE Linux)

Best regards
Michael

[msolal]

Quoting - mdajunge

To compile the provided examples I used
make libem64t compiler=gnu

For my own demo program I took the command generate by make and adapted it. I then compiled from command line with:

gcc -w -I/opt2/intel/mkl/10.1.1.019/include dss_usym_complex_c.c -L"/opt2/intel/mkl/10.1.1.019/lib/em64t" "/opt2/intel/mkl/10.1.1.019/lib/em64t"/libmkl_solver_lp64.a "/opt2/intel/mkl/10.1.1.019/lib/em64t"/libmkl_intel_lp64.a -Wl,--start-group "/opt2/intel/mkl/10.1.1.019/lib/em64t"/libmkl_intel_thread.a "/opt2/intel/mkl/10.1.1.019/lib/em64t"/libmkl_core.a -Wl,--end-group -L"/opt2/intel/mkl/10.1.1.019/lib/em64t" -liomp5 -lpthread -lm -o dss_usym_complex_c.out

I tested both the provided examples and the my demo on two machines:
1)
--Intel Xeon CPU 5160 @ 3.00GHz,
--SLES 10 SP1
--Linux 2.6.16.46-0.12-smp #1 SMP Thu May 17 14:00:09 UTC 2007 x86_64 x86_64 x86_64 GNU/Linux
--gcc (GCC) 4.1.2 20070115 (prerelease) (SUSE Linux)
2)
--Intel Core2 CPU 6300
--Linux 2.6.18.8-0.7-default #1 SMP Tue Oct 2 17:21:08 UTC 2007 x86_64 x86_64 x86_64 GNU/Linux
--openSUSE 10.2 (X86-64)
--gcc (GCC) 4.1.2 20061115 (prerelease) (SUSE Linux)

Best regards
Michael

Hello,
I have a similar problem on windows. I am trying to use the Paradiso DSS interface for complex numbers. I also started from the example dss_sym_c.c and modified it to be complex (with 0 imaginary part).
In the case, I use a symmetric matrix, I get the same result as the one with real numbers and everything is fine. To check the non symmetric case, I took the same case but entering the full matrix. The program never get out of the factorization routine.
I think the problem is not in the linking, since the symmetric case is working. Do you have any working example of using this Paradiso DSS for non symmetric complex matrices?

Thanks,

Marc

Here is the code:


/*
********************************************************************************
* INTEL CONFIDENTIAL
* Copyright(C) 2004-2008 Intel Corporation. All Rights Reserved.
* The source code contained or described herein and all documents related to
* the source code ("Material") are owned by Intel Corporation or its suppliers
* or licensors. Title to the Material remains with Intel Corporation or its
* suppliers and licensors. The Material contains trade secrets and proprietary
* and confidential information of Intel or its suppliers and licensors. The
* Material is protected by worldwide copyright and trade secret laws and
* treaty provisions. No part of the Material may be used, copied, reproduced,
* modified, published, uploaded, posted, transmitted, distributed or disclosed
* in any way without Intel's prior express written permission.
* No license under any patent, copyright, trade secret or other intellectual
* property right is granted to or conferred upon you by disclosure or delivery
* of the Materials, either expressly, by implication, inducement, estoppel or
* otherwise. Any license under such intellectual property rights must be
* express and approved by Intel in writing.
*
********************************************************************************
* Content : MKL DSS C example
*
********************************************************************************
*/
/*
**
** Example program to solve symmetric positive definite system of equations.
** modified for non symmetric
*/
#include
#include
#include
#include "mkl_dss.h"
/*
** Define the array and rhs vectors
*/
#define NROWS 5
#define NCOLS 5
#define NNONZEROS 13
#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, 8, 10, 12, 14 };
static _INTEGER_t columns[NNONZEROS] = { 1, 2, 3, 4, 5,1, 2,1, 3,1, 4,1, 5 };
static _DOUBLE_PRECISION_t values[NNONZEROS] = { 9, 1.5, 6, .75, 3, 1.5, 0.5, 6, 12,.75, .625,3., 16 };
//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 };

MKL_INT main() {
MKL_INT i;
/* Allocate storage for the solver handle and the right-hand side. */
//_DOUBLE_PRECISION_t solValues[NROWS];
_DOUBLE_COMPLEX_t solValues[NROWS];
_DOUBLE_COMPLEX_t matvalues[NNONZEROS],rhsc[NCOLS];
_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;
MKL_INT sym=MKL_DSS_NON_SYMMETRIC;
MKL_INT type=MKL_DSS_INDEFINITE;
int ix[NROWS];
int ii;
for (ii=0;ii<=NNONZEROS-1;ii++){
matvalues[ii].r=values[ii];
matvalues[ii].i=0.0;
}
for (ii=0;ii<=NCOLS-1;ii++){
rhsc[ii].r=rhs[ii];
rhsc[ii].i=0.0;
}

/* --------------------- */
/* Initialize the solver */
/* --------------------- */
error = dss_create(handle, opt );
printf("create donen");


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 );
printf("define donen");


if ( error != MKL_DSS_SUCCESS ) goto printError;
/* ------------------ */
/* Reorder the matrix */
/* ------------------ */

error = dss_reorder( handle, opt, ix);
printf("reorder donen");


if ( error != MKL_DSS_SUCCESS ) goto printError;
/* ------------------ */
/* Factor the matrix */
/* ------------------ */
//error = dss_factor_real( handle, type, values );
error = dss_factor_complex( handle, type,matvalues );
printf("factor donen");


if ( error != MKL_DSS_SUCCESS ) goto printError;
/* ------------------------ */
/* Get the solution vector */
/* ------------------------ */
//error = dss_solve_real( handle, opt, rhs, nRhs, solValues );
error=dss_solve_complex(handle,opt,rhsc,nRhs,solValues);
printf("solve donen");


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.r );
printf("n");
exit(0);
printError:
printf("Solver returned error code %dn", error);
exit(1);
}

[Gennady_F_Intel]

Yes, we have a couple of such examples. Please look at pardiso_unsym_complex_c.c and pardiso_unsym_complex_f.f files. You can find these tests into examplessolversource directory. These examples to show the use of the "PARDISO" routine for complex unsymmetric linear systems.
--Gennady

[msolal]

Quoting - Gennady Fedorov (Intel)

Yes, we have a couple of such examples. Please look at pardiso_unsym_complex_c.c and pardiso_unsym_complex_f.f files. You can find these tests into examplessolversource directory. These examples to show the use of the "PARDISO" routine for complex unsymmetric linear systems.
--Gennady

I tried this unsym_complex_c.c example you give. I was able to compile it and to run it. But this is not your nice DSS interface. I think the DSS interface solver does not work for non symmetric complex matrices. I will have to use the PARDISO routine instead. It is ok for me, but we are two to have similar problems. You should check if we are right and in this case correct in a next version.

Thanks,

Marc Solal

[Gennady_F_Intel]

Marc,
Really, we have run time problem with DSS interface for non symmetric complex matrix.
The cause of the problem dealt with threading and we are going to fix this problem at the next release.

As a workaround, Id recommend you
or use sequential libraries:
for example for win32:
mkl_solver_sequential.lib mkl_intel_c.lib mkl_sequential.lib mkl_core.lib
or
before first MKL routine calling
mkl_set_num_threads( 1 );

I tried to use both of suggested methods and they work on win32:
Below the output:
++++++++++++++++++++++++++++++
define structure
reorder
factor_complex
solve complex
statistics
determinant power is 5
determinant base is 2.5872
Determinant is 258720
delete
Solution array: 0.107468 = i* 0.107468 -0.25 = i* -0.25 0.0568182 = i* 0.0568182 -0.00324675 = i*
-0.00324675 -0.227273 = i* -0.227273 -0.302273 = i* -0.302273 0.113636 = i* 0.113636 0.188636 = i*
0.188636
Press any key to continue . . .
++++++++++++++++++++++++++++++

I hope it helps you with this problem till the fix will availale.
--Gennady

[Gennady_F_Intel]

Mark,
The problem you reported was fixed in the version 10.2.
This version available for download from intel registration center: https://registrationcenter.intel.com/

In order to fix the problem you need to change the second parameter to dss_define_structure, namely:

int sym = MKL_DSS_NON_SYMMETRIC; should be changed to
int sym = MKL_DSS_NON_SYMMETRIC_COMPLEX;

--Gennady