Different results of calculation by DSS-interface on C# with real and complex data where complex = real + 0j

[vd19]

Hello.
I'm trying to use pardiso-based (DSS) routines on C#. To use original routines I've written 2 simple wrappers:
the fist wrapper written on C++ over original routines and built into dll, and the second written on C# over this dll.
So, if I solve a problem with real data - final result is proper, and if I solve a complex problem where all complex data have zero imaginary part (actually it's a real data, written in complex notation) - result is wrong. Correctness of calculated results compared with Mathcad and between itself (results must be equal).
I have inadequate results only in C#, analogous sample in C++ with using original and not wrapped routines returns correct results.

I can't understand where is difference between calling original routines from C++ and wrapped routines from C# considering that
the sample with real data gives a good result. May be in description of complex data-type? It described as structure with the first double field 'r' and the second - 'i'... I don't know.

Where I can made an error?

I don't know how to attach files using forum's possibilities. Four files (2 wrappers and 2 samples) are here: http://webfile.ru/4311060

[Gennady_F_Intel]

see here how attach file to forum post.

[vd19]

Thank you, Gennady.

I already reload files..

[vd19]

So.. I written the combined test of using pardiso and DSS-routines with complex data. They gives different results with using single data: pardiso gives adequate results, DSS - no.
I think it's very important. It may be error of DSS-interface routines.. Or may be this situation is not enough documented. Or may be I don't understand something. It is serious question to Intel's members.
Test code see below.
P.S. DLL with MKL routines was created similary description in forum's topic: http://software.intel.com/en-us/articles/using-intel-mkl-in-your-c-program/.

Results of computing using pardiso directly and via DSS:

The solution of the system solved by pardiso is:

x[0] = 0,107467532467532+0j
x[1] = -0,25+0j
x[2] = 0,0568181818181818+0j
x[3] = -0,00324675324675324+0j
x[4] = -0,227272727272727+0j
x[5] = -0,302272727272727+0j
x[6] = 0,113636363636364+0j
x[7] = 0,188636363636364+0j

The solution of the system solved via DSS is:
0,671346263984748-0,0769951948846796j
-0,734375000000091-0,187500000000026j
0,93572214008168+0,0333629269701237j
0,000780919725403221-0,0784913138809309j
-3,71036534829729-0,561772850460331j
-2,20701388789418+0,0933585467090734j
0,335450319308798+0,0542355371900888j
0,385110852606479-0,00710001241296985j

Main test code:

public static int Main(string[] args)
{
/* Matrix data. */
int[] ia = new int[] { 1, 5, 8, 10, 12, 13, 16, 18, 21 };
int[] ja = new int[] { 1, 3, 6, 7, 2, 3, 5, 3, 8, 4, 7, 2, 3, 6, 8, 2, 7, 3, 7, 8 };
int
n = 8,
nNonZeros = ja.Length;

Complex[] a = new Complex[]
{
new Complex( 7,0), new Complex( 1,0), new Complex( 2,0), new Complex( 7,0),
new Complex(-4,0), new Complex( 8,0), new Complex( 2,0),
new Complex( 1,0), new Complex( 5,0),
new Complex( 7,0), new Complex( 9,0),
new Complex(-4,0),
new Complex( 7,0), new Complex( 3,0), new Complex( 8,0),
new Complex( 1,0), new Complex(11,0),
new Complex(-3,0), new Complex( 2,0), new Complex( 5,0)
};
Complex[] b = new Complex[]
{
new Complex (1,0),
new Complex (1,0),
new Complex (1,0),
new Complex (1,0),
new Complex (1,0),
new Complex (1,0),
new Complex (1,0),
new Complex (1,0)
};
Complex[] x = new Complex;

int
mtype = 13, // Double complex unsymmetric matrix
maxfct = 1, // Maximum number of numerical factorizations.
mnum = 1, // Which factorization to use.
phase,
error = 0, // Initialize error flag
msglvl = 0, // Print statistical information in file
nrhs = 1, // Number of right-hand sides
opt = (int)Options.TerminationLevel.MKL_DSS_TERM_LVL_FATAL,
sym = (int)Options.Structure.MKL_DSS_NON_SYMMETRIC_COMPLEX,
type = (int)Options.Factorization.MKL_DSS_POSITIVE_DEFINITE;

IntPtr[] pt = new IntPtr[64];
IntPtr handle = new IntPtr();

int[] iparm = new int[64];

int i;
Complex[] ddum = new Complex[1];
int[] idum = new int[1];

for (i = 0; i < 64; i++)
iparm = 0;

iparm[0] = 1; /* No solver default */
iparm[1] = 2; /* Fill-in reordering from METIS */
/* Numbers of processors, value of OMP_NUM_THREADS */
iparm[2] = 1;
iparm[3] = 0; /* No iterative-direct algorithm */
iparm[4] = 0; /* No user fill-in reducing permutation */
iparm[5] = 0; /* Write solution into x */
iparm[6] = 0; /* Not in use */
iparm[7] = 2; /* Max numbers of iterative refinement steps */
iparm[8] = 0; /* Not in use */
iparm[9] = 13; /* Perturb the pivot elements with 1E-13 */
iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
iparm[11] = 0; /* Not in use */
iparm[12] = 0; /* Maximum weighted matching algorithm is switched-off
* (default for symmetric). Try iparm[12] = 1 in case of
* inappropriate accuracy */
iparm[13] = 0; /* Output: Number of perturbed pivots */
iparm[14] = 0; /* Not in use */
iparm[15] = 0; /* Not in use */
iparm[16] = 0; /* Not in use */
iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
iparm[18] = -1; /* Output: Mflops for LU factorization */
iparm[19] = 0; /* Output: Numbers of CG Iterations */

for (i = 0; i < 64; i++)
pt = IntPtr.Zero;

//**************************************************
// Using pardiso...
//**************************************************

// Reordering and Symbolic Factorization. This step also allocates all memory that is necessary for the factorization.
phase = 11;
DSS.pardiso(pt, ref maxfct, ref mnum, ref mtype, ref phase, ref n, a, ia, ja, idum, ref nrhs, iparm, ref msglvl, ddum, ddum, ref error);

// Numerical factorization.
phase = 22;
DSS.pardiso(pt, ref maxfct, ref mnum, ref mtype, ref phase, ref n, a, ia, ja, idum, ref nrhs, iparm, ref msglvl, ddum, ddum, ref error);

// Back substitution and iterative refinement.
phase = 33;
iparm[7] = 2; /* Max numbers of iterative refinement steps. */
DSS.pardiso(pt, ref maxfct, ref mnum, ref mtype, ref phase, ref n, a, ia, ja, idum, ref nrhs, iparm, ref msglvl, b, x, ref error);

// .. Termination and release of memory. */
phase = -1; /* Release internal memory. */
DSS.pardiso(pt, ref maxfct, ref mnum, ref mtype, ref phase, ref n, ddum, ia, ja, idum, ref nrhs, iparm, ref msglvl, ddum, ddum, ref error);

Console.WriteLine("The solution of the system solved by pardiso is:");
for (i = 0; i < n; i++)
Console.WriteLine("x[" + i + "] = " + x);


//**************************************************
// Using DSS...
//**************************************************
error = DSS.dss_create(ref handle, ref opt);
error = DSS.dss_define_structure(ref handle, ref sym, ia, ref n, ref n, ja, ref nNonZeros);
error = DSS.dss_reorder(ref handle, ref opt, new int[5]);
error = DSS.dss_factor_complex(ref handle, ref type, a);
error = DSS.dss_solve_complex(ref handle, ref opt, b, ref nrhs, x);
error = DSS.dss_delete(ref handle, ref opt);

Console.WriteLine("\nThe solution of the system solved via DSS is:");
for (i = 0; i < n; i++)
Console.WriteLine("{0} ", x.ToString());


return 0;
}

Wrapper and Complex definition code:

namespace MKL
{

public struct Complex
{
public double r, i;

public Complex(double r, double i)
{
this.r = r;
this.i = i;
}
public override string ToString()
{
string Res;
if (i < 0)
Res = this.r.ToString() + this.i.ToString() + "j";
else
Res = this.r.ToString() + "+" + this.i.ToString() + "j";
return Res;
}
}
public static class Options
{
public enum TerminationLevel
{
MKL_DSS_TERM_LVL_SUCCESS = 1073741832,
MKL_DSS_TERM_LVL_DEBUG = 1073741840,
MKL_DSS_TERM_LVL_INFO = 1073741848,
MKL_DSS_TERM_LVL_WARNING = 1073741856,
MKL_DSS_TERM_LVL_ERROR = 1073741864,
MKL_DSS_TERM_LVL_FATAL = 1073741872
}
public enum Structure
{
MKL_DSS_SYMMETRIC = 536870976,
MKL_DSS_SYMMETRIC_STRUCTURE = 536871040,
MKL_DSS_NON_SYMMETRIC = 536871104,
MKL_DSS_SYMMETRIC_COMPLEX = 536871168,
MKL_DSS_SYMMETRIC_STRUCTURE_COMPLEX = 536871232,
MKL_DSS_NON_SYMMETRIC_COMPLEX = 536871296
}

public enum Factorization
{
MKL_DSS_POSITIVE_DEFINITE = 134217792,
MKL_DSS_INDEFINITE = 134217856,
MKL_DSS_HERMITIAN_POSITIVE_DEFINITE = 134217920,
MKL_DSS_HERMITIAN_INDEFINITE = 134217984
}
}

[SuppressUnmanagedCodeSecurity]
public class DSS
{

const string MKL_DLL_Intel_Name = @"MKL_.dll";

[DllImport(MKL_DLL_Intel_Name, CallingConvention = CallingConvention.Cdecl)]
public static extern int dss_create(ref IntPtr handle, ref int opt);

[DllImport(MKL_DLL_Intel_Name, CallingConvention = CallingConvention.Cdecl)]
public static extern int dss_define_structure(ref IntPtr handle, ref int opt, int[] rowIndex, ref int nRows, ref int nCols, int[] columns, ref int nNonZeros);

[DllImport(MKL_DLL_Intel_Name, CallingConvention = CallingConvention.Cdecl)]
public static extern int dss_reorder(ref IntPtr handle, ref int opt, int[] perm);

[DllImport(MKL_DLL_Intel_Name, CallingConvention = CallingConvention.Cdecl)]
public static extern int dss_factor_real(ref IntPtr handle, ref int opt, double[] Values);

[DllImport(MKL_DLL_Intel_Name, CallingConvention = CallingConvention.Cdecl)]
public static extern int dss_factor_complex(ref IntPtr handle, ref int opt, Complex[] Values);

[DllImport(MKL_DLL_Intel_Name, CallingConvention = CallingConvention.Cdecl)]
public static extern int dss_solve_real(ref IntPtr handle, ref int opt, double[] RhsValues, ref int nRhs, double[] SolValues);

[DllImport(MKL_DLL_Intel_Name, CallingConvention = CallingConvention.Cdecl)]
public static extern int dss_solve_complex(ref IntPtr handle, ref int opt, Complex[] RhsValues, ref int nRhs, [Out] Complex[] SolValues);

[DllImport(MKL_DLL_Intel_Name, CallingConvention = CallingConvention.Cdecl)]
public static extern int dss_statistics(ref IntPtr handle, ref int opt, ref string RequestName, double[] retValues);

[DllImport(MKL_DLL_Intel_Name, CallingConvention = CallingConvention.Cdecl)]
public static extern int dss_delete(ref IntPtr handle, ref int opt);

[DllImport(MKL_DLL_Intel_Name, CallingConvention = CallingConvention.Cdecl)]
public static extern int pardiso(

[In, Out] IntPtr[] handle,
ref int maxfct,
ref int mnum,
ref int mtype,
ref int phase,
ref int n,
[In] Complex[] a,
[In] int[] ia,
[In] int[] ja,
[In] int[] perm,
ref int nrhs,
[In, Out] int[] iparm,
ref int msglvl,
[In, Out] Complex[] b,
[Out] Complex[] x,
ref int error);

}

}

[Vladimir_Koldakov__I]

Hello,

This is not a bugin DSS. This is effect of managed memory control. CLR guarantees pinning of thepassed arrays during the call of the unmanaged method only. It does not guarantee the same location of arrays between two consecutive unmanaged calls. Inthe lastcase you should explicitly pin arrays. Lets see, your original program prints the following:

The solution of the system solved by pardiso is:
x[0] = 0.107467532467532 +0j
x[1] = -0.25 +0j
x[2] = 0.0568181818181818 +0j
x[3] = -0.00324675324675324 +0j
x[4] = -0.227272727272727 +0j
x[5] = -0.302272727272727 +0j
x[6] = 0.113636363636364 +0j
x[7] = 0.188636363636364 +0j

The solution of the system solved via DSS is:
0.671346263984748 -0.0769951948846796j
-0.734375000000091 -0.187500000000026j
0.93572214008168 +0.0333629269701237j
0.000780919725403221 -0.0784913138809309j
-3.71036534829729 -0.561772850460331j
-2.20701388789418 +0.0933585467090734j
0.335450319308798 +0.0542355371900888j
0.385110852606479 -0.00710001241296985j

Ive a bit modified the program:

118a119,120
> GCHandle a_pin = GCHandle.Alloc(a,GCHandleType.Pinned);
> IntPtr a_ptr = a_pin.AddrOfPinnedObject();
126c128
< error = DSS.dss_factor_complex(ref handle, ref type, a);
---
> error = DSS.dss_factor_complex(ref handle, ref type, a_ptr);
129a132,133
> a_pin.Free();
>
211c215
< public static extern int dss_factor_complex(ref IntPtr handle, ref int opt, Complex[] Values);
---
> public static extern int dss_factor_complex(ref IntPtr handle, ref int opt, IntPtr Values);

The result now is the following:

The solution of the system solved by pardiso is:
x[0] = 0.107467532467532 +0j
x[1] = -0.25 +0j
x[2] = 0.0568181818181818 +0j
x[3] = -0.00324675324675324 +0j
x[4] = -0.227272727272727 +0j
x[5] = -0.302272727272727 +0j
x[6] = 0.113636363636364 +0j
x[7] = 0.188636363636364 +0j

The solution of the system solved via DSS is:
0.107467532467532 +0j
-0.25 +0j
0.0568181818181818 +0j
-0.00324675324675324 +0j
-0.227272727272727 +0j
-0.302272727272727 +0j
0.113636363636364 +0j
0.188636363636364 +0j

Note, its better to pin ia and ja as well because mklman.pdf declares:

WARNING. Pointers to the rowIndex and columns are stored in the handle. These arrays are used when the solver runs. Do not free or delete them before calling dss_solve_real or dss_solver_complex.

Thanks,
Vladimir