accuracy on calculating Schur complement with MKL PARDISO

[Jay_9]

Hi all,

I am using MKL PARDISO 2018.

When I calculating the Schur complement with MKL PARDISO 2018, the nonzero pattern of the Schur complement, which is calculated with MKL PARDISO, is different from the reference result, which is calculated with simply PARDISO LU ( you can see the attached file).

Additional I set iparm[35] to -2 and I extracted from CSR3 format.

The result is attached to *.PNG.  the left figure represents the reference result and the right figure represents the result resulting from PARDISO iparm[35]= -2.

I have looked for the bug lists(MKLD-4069 schur complement method accuracy is not enough) in your site(https://www.intel.com/content/www/us/en/developer/articles/troubleshooting/intel-math-kernel-library-2020-bug-fixes.html).

Since the fixed bug might be related to my problem, I would like to know contents of MKLD-4069.

Could you please tell me what the bug is exactly ?

[VidyalathaB_Intel]

Hi,

Thanks for reaching out to us.

>>Could you please tell me what the bug is exactly?

Sorry for the inconvenience caused. Unfortunately, we could not get enough details to share regarding the issue on which the bug was reported, will update you if we get full details. 

More or less the issue is related to accuracy it seems and as mentioned in the article the issue is fixed in MKL 2020.0 release.

Meanwhile, you can give it a try with the latest oneMKL 2022 which is available for download. You can get that by downloading oneAPI Base Toolkit (https://www.intel.com/content/www/us/en/developer/tools/oneapi/base-toolkit-download.html)

You can also share your reproducer code and steps to reproduce the issue (along with your OS details) so that we could test it with the latest MKL version on our side and see if the results are still different from the reference results that you have provided.

Regards,

Vidya.

[Jay_9]

Hi Vidya,

Thanks for your reply.  As you recommended, I downloaded the latest oneMKL 2022 that you mentioned and calculated Schur complement.

However, the result that are computed with the latest Version is different from the reference. 

I don't know what setting I did wrong. Is it possible to check it if I provide the input matrix?  If you can check it, should I provide the input matrix with CSR matrix format and *.cpp file that I wrote?

Best regards,

Jay  

[VidyalathaB_Intel]

Hi,

Thanks for getting back to us.

>>setting I did 

Could you please check the value of iparm[12]? we can try setting iparm[12]=1 in case of inappropriate accuracy.

Please refer to the below link for more details that might help you

https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-iface/pardiso-iparm-parameter.html#pardiso-iparm-parameter_IPARM37

You can also check the example codes for calculating Schur complement which is present inside the installed MKL directory as shown "oneAPI\mkl\2022.0.0\examples\examples_core_c\c\sparse_directsolvers\source"

>>Is it possible to check it if I provide the input matrix? If you can check it, should I provide the input matrix with CSR matrix format and *.cpp file that I wrote?

If you still get any issues, feel free to provide us with the .cpp file and the input matrix along with the compilation and execution commands you have used so that we could check it on our end and assist you accordingly.

Regards,

Vidya.

[Jay_9]

Hi Vidya,

I set iparm[12] = 1 as you said. But the calculated Schur complement is still inaccurate.

I have attached a zipfile. 

In azip file, CSR.mat represents the input matrix. you can see the file format in schur_examp.xpp. You can just read  it with the file that I haved offered,

schur_reference.mat represents the Schur complement  that is calculated without using iparm[35] = 2. It simply makes the Schur complement direct by factorizing A11.

schur.mat is the calculated Schur complement with iparm[35] =2.

schur_examp.cpp is simple code to calculate the Schur complement.  In this code, input matrix is read to calculate the Schur complement.

If you look into schur_reference.mat and schur.mat, you can find that the non-zero pattern of Schur complement calculated with iparm[35] =2 is different from the non-zero pattern of Schur complement calculated without iparm[35] =2

Could you tell me why there is difference between them? 

Best regards,

Jay

[Jay_9]

Hi Vidya,

Have you checked the file that I prepared?   Did I set the iparm wrong? 

When I solved a matrix system with PARDISO iparm[35]=2, the solution is inaccurate.

I have been debugging to find the cause of inaccurate. I checked the input matrix and iparm. I can not find what is wrong.

Could you tell me which one is wrong?

I am sorry to bother you. 

Best regards, 

Jay

[VidyalathaB_Intel]

Hi,

>>schur_reference.mat represents the Schur complement that is calculated without using iparm[35] = 2.

Since you are comparing the Schur.mat values with schur_refernce.mat values and calculating the difference between them,

Could you please help us in interpreting the statement regarding Without using iparm[35] = 2?

For more details regarding iparm[35] parameter, please refer to the below link

https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-iface/pardiso-iparm-parameter.html

Regards,

Vidya.

[Jay_9]

Hi,

>>> Could you please help us in interpreting the statement regarding Without using iparm[35] = 2?

  -  As I mentioned the previous (It simply makes the Schur complement direct by factorizing A11) to you, I meant that "Without using iparm[35] = 2" directly computes the Schur complement. 

     When a 2 by 2 block matrix ( A = [ A11 A12 ]

                                                                  [ A21 A22] )  is given, 

     I directly computed the Schur complement as follows

     S = A22 - A21 * inv(A11) * A12 , 

     where inv(A11) represents the inversion of A11

    I also attached the sample code in order to compute the reference data.  you can run the sample code with octave. 

    the Schur complement computed with octave is different from one computed with iparm[35] =2 in PARDISO

    I think the Schur complement is inaccurate as the Schur complement is computed with the iparm[35]=2. 

    If you directly compute the Schur complement ( A22 - A21 *inv(A11)*A12 ) by factorizing A11, you can check that the Schur

    complement, which is computed with iparm[35]=2, is different from the result computed from the other methods.

    you also can run *.cpp, which I provided to you, with visual studio and gcc.   

   Could you please check why there is difference?    

Regards,

Jay

[VidyalathaB_Intel]

Hi Jay,

Thank you for your patience in sharing the details and I apologize for the delay in response from my end.

We are able to reproduce your results from our end as well with iparm[35]=2 by running the schur_exam.cpp with the latest oneMKL 2022.

Before passing this issue to the concerned development team we would like to get some more clarifications.

>>result computed from the other methods

As part of that, we just need the details for 3rd scenario which you have given in that screenshot regarding Schur complement computed directly with pardiso by factorizing A11.

Could you please share the code for the above 3rd scenario so that we can reproduce it and report it to our experts?

Regards,

vidya.

[Jay_9]

Hi,

Thanks for your reply.  

 >> Could you please share the code for the above 3rd scenario

I have attached a zip file.  You can see 5 files in the zip file.  CSR.mat and schur_reference.mat are the files that I offered you previously.

In shur_examp.cpp, I added the sample code that you requested. 

If you run shur_examp.cpp, 2 files are created.  "direct_pardiso_schur.mat" is the result that is directly computed with pardiso by factorizing A11.  "direct_pardiso_schur.mat" is created in Calculate_Schur_Complement function.

"pardiso_schur.mat" is the result that is computed with iparm[35]=2 in PARDISO

If you look into two files(direct_pardiso_schur.mat and pardiso_schur.mat), you can check the results in the screenshot 

Let me know if you have any questions.

Regards,

Jay

[VidyalathaB_Intel]

Hi Jay,

Thanks for providing the details.

We have forwarded the issue to our experts and will get back to you soon with updates regarding the issue.

Thanks for your patience.

Regards,

Vidya.

[Khang_N_Intel]

Hi Jay,

Just want to let you know that we are working on this issue.

Best,

Khang

[Jay_9]

Hi Khang, 

Have you ever found the cause of the issue?

In my application, I need to calculate Schur complement using block LU for speed up.

Could you please let me know in which cases the inaccuracy occurs?

Thanks,

Jay

[Khang_N_Intel]

Hi Jay,

The mkl team is still working on this issue.

I will let you know as soon as I hear any news.

Best regards,

Khang

[Ruqiu_C_Intel]

Hi Jay,

Thank you again for reaching out to us.

The differences between the reference result for Schur complement and Pardiso output can be attributed to expected absence of zeros in Pardso output. 

Best regards,

Ruqiu