PARDISO, reuse symbolic factorization

I installed the matlab package from pardiso-project.org, and I built the .mex linking with libpardiso600-WIN-X86-64.lib.

I'm solving in a loop a series of problems Ax=b, where A changes, but still has the same sparsity pattern. To exploit that, in the initialization, I only perform symbolic factoriazation once (which supposes to allocate the data structures), and I free it only once; this is the general idea:

pardiso_info = []
bSolverInitialized = 0;

for i = 1:n
    A = As{i};
    rhs = rhses{i};
    
    A = A + eps*speye(size(A)); % diagonal must be full
    A = tril(A);

    % initialize
    if ~bSolverInitialized
        bSolverInitialized = 1;
        pardiso_info = pardisoinit(-2,0);

        % Analyze the matrix and compute a symbolic factorization.
        pardiso_info = pardisoreorder(A, pardiso_info, false);
    end

    % Compute the numeric factorization
    pardiso_info = pardisofactor(A, pardiso_info, false);

    % Compute the solutions using the symbolic factorization
    [sol, obj.pardiso_info] = pardisosolve(A, rhs, pardiso_info, false);
end

pardisofree(pardiso_info);
clear pardiso_info;

I execute this code many times (in an outer loop). Matlab memory usage increases until it hangs. I suspect a memory leak, and pardiso is my prime suspect. I'll need to do further testing to try and pinpoint the problem, but I thought to verify first that I'm using it correctly since the manual didn't mention reusing symbolic factorization.

I'm not sure if it's the right place for the question. pardiso-project.org doesn't seem to have a forum, but Intel may support it as part of MKL.

[Yiyeung Lau]

Hi Zohar, could you please elaborate how you solve the memory issues? I'm using the pardiso's MATLAB interface too (on a Ubuntu system), and I encounter fatal error from MATLAB whenever I call the pardisofactor function and when A is not diagonal (the function works fine if A is diagonal). The error message reads:

malloc.c:4036: _int_malloc: Assertion `(unsigned long) (size) >= (unsigned long) (nb)' failed

Do you have any idea how to deal with this? Thank you very much!

[Zohar]

Sorry, no.

Answer 1

[Zohar]

I played with it a bit. I tried tracking the memory allocations using 'profile -memory on' and 'memory'. Unfortunately, the memory management in matlab is a fickle thing:

• When an array is allocated, some of the memory can be reserved.
• When you fill an array, allocation may occur.
• There's the garbage collection who's reponsible to free the memory, and it's unpredicatable.

Pardiso had nothing to do with the leak, and what I did seems fine. Pure luck led me to my logger. Through a socket, I send messages to a local app (since I'm calling matlab from cpp, and I want to see output). I used tcpip for that. I didn't see anywhere mentioned that delete should be called in the end:

https://au.mathworks.com/help/instrument/tcpip.html

Also, adding delete didn't help.

Now, I'm using instead tcpclient, and for now, the leak seems to be resolved.

[Jeremy Hu]

Hi Zohar, glad to see you got this working. I am also trying to use this version of pardiso in matlab.

I was wondering, when you compiled in mex, were you able to use gfortran in the mex makefile as in the example provided on the website, or did you go with some version of Intel ifort?

Best, Jeremy

[Zohar]

I used the supplied 'libpardiso600-WIN-X86-64.lib', and it didn't need anything else.

[Jeremy Hu]

Thanks Zohar. Does this mean that for pardisoinit you compiled like this:

mex -L'C:...\pardiso-library'...
    -llibpardiso600-WIN-x86-64 -output pardisoinit...
    common.cpp matlabmatrix.cpp sparsematrix.cpp pardisoinfo.cpp pardisoinit.cpp

and similarly for the other functions (pardisosolve, pardisoreorder, pardisofactor, etc.)

The reason I'm asking is because it seems like we need to link with a fortran compiler such as intel or gfortran for it to work.

[Zohar]

My mex is set up for vs2019, and I built it with:

names = {'pardisoinit' 'pardisoreorder' 'pardisofactor' 'pardisosolve' 'pardisofree'};
for i = 1:length(names)
    name = names{i};
    mex('-llibpardiso600-WIN-X86-64', '-output', name, 'common.cpp', 'matlabmatrix.cpp', 'sparsematrix.cpp', 'pardisoinfo.cpp', [name '.cpp'])
end

[Jeremy Hu]

Ah ok cool, thanks again Zohar!