Matlab bindings for pardiso

[Zohar]

Are there convenient bindings for matlab to call pardiso, e.g. like the ones in pardiso-project.org that are based on Peter Carbonetto's code?

[mecej4]

The situations of the Pardiso Project's version and the MKL version of their Pardiso solvers are quite different. For the former, the Pardiso solver is the entirety of their product. On the other hand, Pardiso is a tiny part of the very large set of capabilities that MKL offers. Note also that the Matlab interfaces available through a link from pardiso-project-org were prepared by users at a university. 

Depending on how many people wish to directly call Pardiso from Matlab, someone could volunteer to write the appropriate wrappers.

Should that be Intel? I do not know. If "yes", why not develop wrappers for many of the other excellent solvers in MKL, such as TRNLSP, for example? Why not for all of MKL? What would be the cost of building all those wrappers and keeping them updated to work with various versions of Matlab, MKL and platforms?

[Zohar]

I'll take it as a no. Here, I modified their bindings to work with mkl's version:

https://drive.google.com/file/d/1ACLbguq1Il-wYKaWJZ2tnXbxefa1uwDv/view

Can you please take a look and see why examplesym.m fails with (the unhelpful) "error type -1" (line 32, pardisoreorder)?

[mecej4]

Re "I'll take it as a no": I am just another user, and so I don't speak for Intel. I am not a C++ programmer, either.

Here are a couple of comments about your modified files. The argument dparm in  the Basel/Lugano Pardiso is not present in the MKL version. You have to remove it wherever it is used as an argument. The two versions of Pardiso may differ as to the meanings of the large number of arguments that are passed through the iparm array.

If you look at the various phases in the Pardiso calls in the C or Fortran examples in the MKL documentation, you can see that in many places dummy arguments are passed. It is not easy to understand from the documentation when a dummy argument is acceptable (because it will not be used in that phase) and when it is needed and needs to be passed to a later phase.

Of the four example m-files provided by Carbonetto, the real-symmetric example has one feature that the other three m-files do not: using a user-provided permutation vector. That part is the one that causes the failure that you reported when the MKL version of Pardiso is used. Remove that request, and this example also works.  The other three example m-files run without such a modification.

I did my testing using Matlab 2013a, which is the latest version that I have.

[Zohar]

Strange, when I run exampleunsym.m, I get:

Error using pardisoreorder
PARDISO reports an error of type -1
Error in exampleunsym (line 28)
info = pardisoreorder(A,info,verbose);

Can you confirm that you built the mex and used my code as is (and you didn't use pardiso-prj by mistake, etc.)?

[mecej4]

I did not use your code. I took Carbonetto's  C++ sources, modified them by removing DPARM in the Pardiso calls, and changed the Pardisoinit calls to match MKL. I also removed the provision for passing a permutation vector in the C++ code, and removed the corresponding code in examplesym.m.

[Zohar]

I'm not passing dparam, else it won't compile.

Also, I refer to exampleunsym.m, which doesn't use a permutation, as you pointed out.

Can you please try my code or release yours so I'll have something working?

[mecej4]

Attached is a Zip file containing my modifications to Carbonetto's package for calling Pardiso from Matlab, to enable the MKL version of Pardiso to be called from Matlab. Extract the files to a working directory. Open a 64-bit OneAPI command window, change to the working directory, set the environment to build Mex files, and run the batch file buildall.bat . Five Mexw64 files should be built. You can then run Matlab and execute the four example m-files.

As I stated earlier, I am not a C++ user and I am not prepared to support users of this combination of MKL and Matlab. I have not tested the package with any versions of Matlab other than 2013a-64 bit on Windows. Since my version of Matlab is old, its mex build facility does not work, so I use batch files as noted above.

[Zhongmin]

Dear mecej4:

I have install oneAPI, and  I am using oneAPI2022,Vs2019 and matlab 2018a, 

I build the files as your recommended steps, and got an error,

{ pardisoreorder.cpp(1): catastrophic error: cannot open source file "mex.h"
#include "mex.h" }

I think I miss the step "set the environment to build Mex files",

But I am a new in oneAPI, colud you tell me how to set the environment ?

And due to use the newer vision of one API and Vs., should I update my MATLAB to 2020 or later?

Thanks so much

[Zohar]

I don't see significant difference between your code and mine. I built and ran it anyway, and got the same error message.

I'm using matlab2019b.

I tried to build the mex using 3 different tools:

- matlab's mex 

- vs2019 cl.exe

- intel oneApi 2022 icx.exe

Can anyone try and verify this?

I haven't tried a stand-alone c++ app that does something similar (but uses its own matrices). Do you have, perhaps, a quick example to do a sanity check?

[Zohar]

Sanity check failed: I just tried to build and run 'pardiso_sym_c.c', and I got the same

ERROR during symbolic factorization: -1

At least, I reduced the problem.

[Zohar]

I solved it. Linking with mkl_intel_ilp64.lib gives the error, linking with mkl_intel_lp64.lib runs fine.

"The Intel(R) Math Kernel Library (Intel(R) MKL) ILP64 libraries use the 64-bit integer type (necessary for indexing large arrays, with more than 2³¹-1 elements), whereas the LP64 libraries index arrays with the 32-bit integer type."

So, of course, I needed the latter lib. However, I recalled that it crashed matlab while the other one didn't. Perhaps I changed things since then and forgot about this option.

Thanks @mecej4 

[mecej4]

The batch files that I provided would have linked to the MKL LP64 libraries. As this thread shows, there are so many ways that building Mex files can go wrong, and there are no simple methods to untangle those problems.

[Zohar]

Yes, I had to modify your batch files, and the LP64 libs were one of them (until I changed it back):

bld_pre.bat

icx /O2 /c -I"C:\Program Files\matlab\R2019b\extern\include" *.cpp

bld.bat

icx /LD %1.obj common.obj matlabmatrix.obj sparsematrix.obj pardisoinfo.obj libmx.lib libmex.lib /Fe%1.mexw64 /link /LIBPATH:"C:\Program Files\matlab\R2019b\extern\lib\win64\microsoft" /LIBPATH:"C:\Program Files (x86)\Intel\oneAPI\mkl\2022.0.0\lib\intel64" /export:mexFunction mkl_intel_lp64_dll.lib mkl_core_dll.lib

bldall.bat

for %%f in (free,factor,init,solve,reorder) do call bld pardiso%%f

I don't think it's matlab's fault here or anything to do with mex, as my sanity check discovered. Initially, if you recall, I said "returns (unhelpful) type -1." Somewhere in the docs, there should be a warning about it.

[Zohar]

There's a problem with the symmetric solver.

When I run examplesym.m, it's fine:

The factors have 13 nonzero entries.
The matrix has 3 positive and 1 negative eigenvalues.
PARDISO performed 0 iterative refinement steps.
The maximum residual for the solution X is 4.46e-14.

When I load my matrix instead, adding the following lines at line 22:

load ab.mat;
A = A + eps*speye( size(A) );
B = full( B );

https://drive.google.com/file/d/1AzG573417J7X6vJqHLnqEqf2PSaqa5OR/view

I get:

The factors have 495 nonzero entries.
The matrix has 28 positive and 16 negative eigenvalues.
PARDISO performed 0 iterative refinement steps.
The maximum residual for the solution X is 6.04e+03.

If I take a smaller eps, e.g. 1e-5, the residual decreases.

[Zohar]

An easier way to reproduce the problem is to modify A in examplesym.m:

A = A - diag(A) + eps*speye( size(A) );

and just use

B = ones( n, 1 );

Output:

The factors have 17 nonzero entries.
The matrix has 2 positive and 2 negative eigenvalues.
PARDISO performed 0 iterative refinement steps.
The maximum residual for the solution X is 0.5.

[VidyalathaB_Intel]

Hi Zohar,

Thanks for reaching out to us.

>>Are there convenient bindings for matlab to call pardiso

We don’t plan to provide any kind of bindings for any external tool/software including Matlab.

>>I solved it.

Glad to know that you got the solution to the issue with your own implementation.

The thread is closing and we will no longer respond to this thread. If you require additional assistance from Intel, please start a new thread. Any further interaction in this thread will be considered community only.

Regards,

Vidya.

[Zohar]

I did a sanity check and tested the matrix with Eigen bindings. I used my generic cmake, so remove unnecessary stuff that you don't have.

https://drive.google.com/file/d/1BqfWXPF2qP1W9oxqbyRF4-BflVOC6XO6/view

Now, need to combine this with the mex code to pinpoint the issue.

[Zohar]

In the above example, I meant

A = A - diag(diag(A)) + eps*speye( size(A) );

[Zohar]

Note to self:

My conclusion was that I need to write a new matlab wrapper from scratch using mlx and Eigen bindings (when I need this again).

https://au.mathworks.com/matlabcentral/answers/1731480-mlx-pass-a-sparse-matrix