Wrong solution at phase 333 with schur complement option

[Niels_G]

Hello,

Here are information that you often request:
- OS details and hardware :
Description: CentOS Linux release 7.8.2003 (Core)
Release: 7.8.2003
Intel(R) Xeon(R) Gold 6252 CPU @ 2.10GHz (max RAM 252G)
- Version Intel MKL : 2020.1.217

I want to compute the schur complement of a symmetric sparse matrix, store it in CSR format and solve the linear system with phases 331,332 (with pardiso in this example but it will be done with a cutom subroutine in the future) and 333.
The backward solve phase (333) gives wrong result with a very small symmetric case (attached to this message). I have a python script that gives me the Schur matrix (which is correct) and the result of the linear system. Intermediate results from phase 331 and 332 are correct but phase 333 gives wrong result on the first componant of the soution. It should be [ 0.25 -0.25 1.5 ] and pardiso gives [ 0.3 -0.25 1.5 ]. I have tried many set of parameters (mtype, iparm) but the one that I sent you is the most relevant for me and the results are wrong so far. If you have any solution let me know!

Python Script :

import numpy as np
print("\n ")
print("\n ")
print("\n ")

M = np.matrix('5.0,1.0,0.0; 1.0,9.0,2.0; 0.0,2.0,1.0')

A = M[0:2,0:2]
B = M[0:2,2]
D = M[2,2]

invD = 1/D
S = A - np.dot(B,np.dot(invD,B.T))

print("\nSchur Complement S = A-B.inv(D).B^T =")
print(S)

b=np.array([1,1,1])
x = np.linalg.solve(M, b)
print("\n Solution Ax=b :")
print(x)

print("\n Residual Ax=b :")
print(M.dot(x)-b)

[ShanmukhS_Intel]

Hi Niels,

Thanks for posting on Intel Communities.

Thanks for sharing the source codes. We have tried running the same in our environment.

We could see the below results while running the same. We will get back to you soon with more information regarding the results wrt. shared Python code and Fortran schur 333 implementation. In addition, we recommend you use the latest version of oneAPI with supported OS environments.

Schur Complement S = A-B.inv(D).B^T =

[[5. 1.]

 [1. 5.]]

 Solution Ax=b :

[ 0.25 -0.25  1.5 ]

 Residual Ax=b :

[[0.00000000e+00 4.44089210e-16 2.22044605e-16]]

---------------------------------!

 ------- Output Restitution ------!

 ---------------------------------!

 SOL(      1 ) =  0.300000000000000   

 SOL(      2 ) = -0.250000000000000   

 SOL(      3 ) =  1.50000000000000

Best Regards,

Shanmukh.SS

[ShanmukhS_Intel]

Hi Niels,

We are working on your issue internally. We will get back to you soon with an update.

Best Regards,

Shanmukh.SS

[Niels_G]

Ok Thank you.

Niels

[ShanmukhS_Intel]

Hi Niels,

As we were able to reproduce the issue at our end, our team is working on it internally. We will get back to you with an update soon.

Best Regards,

Shanmukh.SS

[ShanmukhS_Intel]

Hi Niels,

Thanks for reporting this issue. We were able to reproduce it and we have informed the development team about it. We will get back to you soon with an update.

Best Regards,

Shanmukh.SS

[Niels_G]

Hi,

As it is no resolved for the moment, I calculated the Schur complement in dense format with iparm(36) = 3 and all phases 331, 332 and 333 are correct. I also tried to use pardiso_handle_store after phase=12 to store factorisation but I get the following error :

forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image PC Routine Line Source
a.out 000000000040694A Unknown Unknown Unknown
libpthread-2.17.s 00002B17D682C630 Unknown Unknown Unknown
libmkl_core.so 00002B17D360CDE0 mkl_pds_lp64_pard Unknown Unknown
a.out 00000000004055BA Unknown Unknown Unknown
a.out 00000000004048E2 Unknown Unknown Unknown
libc-2.17.so 00002B17D6A5B555 __libc_start_main Unknown Unknown
a.out 00000000004047E9 Unknown Unknown Unknown

My question is : Subroutines "pardiso_handle_store" and "...restore" are available if schur complement options are activated?

You will find attached to this message, an example of fortran code that points out this issue.

Best regards,

Niels