Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

Matlab bindings for pardiso

Zohar
Novice
1,993 Views

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?

0 Kudos
19 Replies
mecej4
Honored Contributor III
1,975 Views

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?

0 Kudos
Zohar
Novice
1,967 Views

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)?

0 Kudos
mecej4
Honored Contributor III
1,937 Views

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.

0 Kudos
Zohar
Novice
1,928 Views

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.)?

 

0 Kudos
mecej4
Honored Contributor III
1,917 Views

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.

0 Kudos
Zohar
Novice
1,902 Views

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?

0 Kudos
mecej4
Honored Contributor III
1,884 Views

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.

0 Kudos
Zhongmin
Beginner
1,509 Views

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

0 Kudos
Zohar
Novice
1,869 Views

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?

0 Kudos
Zohar
Novice
1,866 Views

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.

0 Kudos
Zohar
Novice
1,863 Views

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 231-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
Honored Contributor III
1,841 Views

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.

0 Kudos
Zohar
Novice
1,834 Views

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.

0 Kudos
Zohar
Novice
1,828 Views

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.

0 Kudos
Zohar
Novice
1,824 Views

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.

0 Kudos
VidyalathaB_Intel
Moderator
1,804 Views

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.

 

0 Kudos
Zohar
Novice
1,788 Views

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.

0 Kudos
Zohar
Novice
1,767 Views

In the above example, I meant

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

0 Kudos
Zohar
Novice
1,441 Views

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

0 Kudos
Reply