PARDISO: Suspected bug in partial / Schur-complement solve

[rbunger]

Hi,

I have a complex symmetric system (mtype = 6). I want to compute a partial solution, specifically just a small upper part of the solution vector. I have tested two different matrices:

1. For the first one (ID = 746704) both the partial solve mode (iparm[30] =3) and the sparse Schur-complement mode (iparm[35] = -2) deliver the same correct result.
2. The second matrix is special: It contains two individual sub-matrices so that the corresponding systems could be solved individually, like a block-diagonal matrix (the matrix is a FEM matrix where two domains are not connected (they are connected only via a hybridization with an integral equation)). The computed reduced solutions are wrong.

Each matrix is about 25 MB in size (compressed), so I can't upload them.

Below, I just print the first six coefficients of the solution vector. "partial_solve = 0": Schur-complement mode; "partial_solve = 1": Partial solve mode. 

Thanks & regards,

Rainer

---

ID = 746704
partial_solve = 0
N_schur = 12341

Reference solution:
x(0) = (0,26.0314)
x(1) = (0,-25.8455)
x(2) = (0,41.7267)
x(3) = (0,-67.5421)
x(4) = (0,16.0706)
x(5) = (0,-16.0299)

Reduced solution:
x(0) = (0,26.0314)
x(1) = (0,-25.8455)
x(2) = (0,41.7267)
x(3) = (0,-67.5421)
x(4) = (0,16.0706)
x(5) = (0,-16.0299)

Full solution:
x(0) = (0,26.0314)
x(1) = (0,-25.8455)
x(2) = (0,41.7267)
x(3) = (0,-67.5421)
x(4) = (0,16.0706)
x(5) = (0,-16.0299)

---

ID = 746704
partial_solve = 1

Reference solution:
x(0) = (0,26.0314)
x(1) = (0,-25.8455)
x(2) = (0,41.7267)
x(3) = (0,-67.5421)
x(4) = (0,16.0706)
x(5) = (0,-16.0299)

Reduced solution:
x(0) = (0,26.0314)
x(1) = (0,-25.8455)
x(2) = (0,41.7267)
x(3) = (0,-67.5421)
x(4) = (0,16.0706)
x(5) = (0,-16.0299)

======

ID = 202443
partial_solve = 0
N_schur = 6508

Reference solution:
x(0) = (0,-2.94291)
x(1) = (0,4.24624)
x(2) = (0,-7.5356)
x(3) = (0,8.64307)
x(4) = (0,-2.44929)
x(5) = (0,-3.66346)

Reduced solution:
x(0) = (0,0.124542)
x(1) = (0,0.13367)
x(2) = (0,0.133519)
x(3) = (0,0.151462)
x(4) = (0,0.152197)
x(5) = (0,0.148365)

Full solution:
x(0) = (0,0.365116)
x(1) = (0,0.589676)
x(2) = (0,-1.5622)
x(3) = (0,2.30311)
x(4) = (0,-0.417001)
x(5) = (0,0.473597)

---

ID = 202443
partial_solve = 1

Reference solution:
x(0) = (0,-2.94291)
x(1) = (0,4.24624)
x(2) = (0,-7.5356)
x(3) = (0,8.64307)
x(4) = (0,-2.44929)
x(5) = (0,-3.66346)

Reduced solution:
x(0) = (0,0.124542)
x(1) = (0,0.13367)
x(2) = (0,0.133519)
x(3) = (0,0.151462)
x(4) = (0,0.152197)
x(5) = (0,0.148365)

[rbunger]

I wanted to attach my test code but I got the following error: "The attachment's pardiso_schur.cpp content type (text/x-c++src) does not match its file extension and has been removed."

[rbunger]

oneAPI 2026.0.

[Aleksandra_K]

Hi, Rainer

Thank you for reaching out. Could you please try uploading the code again using a .c or .cpp file? If needed, try compressing it or upload it as a text file.

Regards,

Alex

[rbunger]

Hi Alex, many thanks for your interest. I will upload the code next week, but not before Tuesday. It is more or less a combination of existing demo codes. Since it works for one of the matrices, it should be ok. I think PARDISO analyses the matrix for special features, independent block solve for example, and these special cases are possibly not properly combined with partial / Schur complement solve...

Thanks, Rainer

[rbunger]

Hi,

I have attached the code as an archive. This is the version I have used to generate the above output. You can't run the code because I can't upload a necessary sparse matrix class I use to generate the reference. I will try to strip down the code and upload it today or tomorrow.

Thanks

Rainer

[rbunger]

Hi,

here is a code version which just depends on the Eigen library...

Thanks for your efforts,

Rainer

[rbunger]

Hi again,

here is a code version without any external requirement (apart from the STL).

Rainer

[rbunger]

Above, I have tested the "partial solve" feature and the Schur-complement feature. Apart from the problem with the wrong result for one of the two matrices, I have compared the PARDISO outputs concerning used resources: They were identical, which is why I think that the partial solve feature is essentially the Schur-complement with a sparse Schur matrix.

I have tested the partial solve feature for another bigger case which delivers the correct result (mtype=6):

• Partial solve factorization step is orders of magnitude slower than the standard factorization (can't give a factor because I measure in seconds and the standard factorization is < 1 sec.)
• The factorization for the partial solve needs about 21 times more RAM
• I call the "solve" step frequently inside an outer iterative solver (another Schur-complement formulation). With partial solve, the outer iterative solver is about 39 times slower. Here my aim was to get a speedup because the outer solver needs just a "surface solution" from PARDISO instead of the full 3D solution

Summary: The partial solve / Schur-complement feature is much slower and much more memory intensive than the std. factorization and solve.

I'm not an expert but it is understandable since a number of (sparse) matrix-matrix products are involved and an efficient implementation seems to be very complicated. This is also explained in a PARDISO manual by O. Schenk.

Regards,

Rainer