Program is either killed or produces segfault with pardiso inverse

I'm hoping to get some resolution on an interesting issue I am running into. Specifically, I am running into issues when attempting to invert a matrix using pardiso. In my case, the process running my code is either killed or produces a segfault.

In conjunction with NGSolve, I am running the following:
-Python 3.6.10, compiled with gcc-9.2.0
-CMake version 3.16.2, compiled with gcc-9.2.0
-NGSolve-6.2.1910-28-g130a6d7, compiled with gcc-9.2.0
-NGSolve is compiled so that pardiso and UMFPACK are enabled via -DUSE_MKL and -DUSE_UMFPACK

The script resolvent_test.py is the script I am running to reproduce the issue, and contains a number of parameters in the
__main__ section at the bottom of the script. The parameters I am concerned with are as follows:
-inverse: A string that specifies the driver for the computation of the inverse
-hermitian_conjugate: A flag that indicates whether to compute the hermitian conjugate of the
resolvents (z*B - A), where A and B are complex square matrices.

To run the script, run python3 resolvent_test.py with inverse = '' or inverse = 'pardiso'. This will result in the
process stalling when inside of the apply_resolvents() method, after which I get one of two single-line messages: 'Killed' or 'Segmentation fault'.
-This happens regardless of whether hermitian_conjugate is True or False. When setting
inverse = 'umfpack' or 'sparsecholesky', however, everything works regardless of whether
hermitian_conjugate is True or False.

I'm wondering if there is something I am missing in the compilation step of NGSolve. Either way, it seems like a strange issue to run into. Any help on the matter is greatly appreciated.

-Ben

It seems that the pardiso inverse functionality began to fail in this context in the commit 32d05dc on December 27. I tested multiple builds, based on different commits:

for these three versions:
bfd5ae04275efb4acbc847f7507d1c7d458105cb January 13
68c745993bf678de680eda74978e8891788f7f40  January 3
32d05dc5d245ab80d69b3201234b9c9fb018f6e6 December 27
Benjamin's script works with umfpack and fails with pardiso.

for these versions:
97bcafd2a2b585f6c990c6259dc31d9471c7b0eb from December 27
869decbc7ade6c070ab069ad99502cc41185f894 from December 27
8382ce7ffd28b0a58f2a4e1e600e97439db5cc16 from December 26
73b3a455c44559a47d1f1a002d897bb1b3911b9a on December 23
Benjamin's script works with both umfpack and pardiso.

Best,
Dow

Hello Ben and Dow,

Thanks for the report. The issue, as you spotted out correctly, started with commit
32d05dc5d245ab80d69b3201234b9c9fb018f6e6

The problem is that PardisoInverse implements neither MultTransAdd() nor MultTrans(). The default implementation of either function calls the other, which results in a stack overflow. I will let you know, when this is fixed.

Best,
Matthias

In the meantime, because your matrices are actually symmetric, as a workaround, in line 74, instead of you could use Kind of disgusting, but it works (for this case!)

Best,
Lukas

Matthias,

I appreciate the reply. I'll plan to move to a previous commit in the meantime.

-Ben

Hey Lukas,

Thank you for the heads up, much appreciated. I initially included the .H since there will be situations where the matrices I am working with are non-hermitian. I just wanted to make sure I had my bases covered before submitting a new topic.

-Ben