inertia computation from magmablas_dsiinertia returns wrong values

[Nai-Yuan Chiang]

Magma version: 2.6.1

After using magma_dsytrf_gpu to factorize a 42x42 matrix, function magmablas_dsiinertia returns wrong inertia computation.

I saved the matrix in a matlab file (attached), in which I compute the eigenvalues from Matlab eig() funciton, and compute the inertia from the factorization factors obtained from magma_dsytrf_gpu. Both of them return the expected value (18,24,0).

However, routine magmablas_dsiinertia returns (17,22,3).

‌

[Stan Tomov]

Thanks for this bug report. I managed to reproduce it.

At first I made the matrix symmetric in matlab, before I saved it. With that matrix magma_dsytrf_gpu produces the correct result (18,24,0) for either lower or upper triangular matrix specified.

If I don’t make the matrix explicitly symmetric, i.e., pass just the lower triangle and use “L” option, I get indeed (17,22,3). However, if I store the matrix as upper triangular and compute with magma_dsytrf_gpu with “U” option, I get the correct result (18,24,0).

At this point I am thinking that there is a bug somewhere in “L” version, e.g., reading from the upper triangular part. However, using LAPACK also gives me the “U” correct and the “L” wrong. The MAGMA code is hybrid and in the “L” case uses LAPACK for panels, so if there is indeed a bug in LAPACK, MAGMA will also be wrong. This will require some more testing to confirm if the bug comes from LAPACK or is something in MAGMA.

[Stan Tomov]

It looks like the bug comes from the LINPACK dsidi routine. This routine works with factorization coming from dsifa, which assumes that the matrix is given only by the upper triangular part, so LINPACK did not support lower triangle storage. Now that we support lower storage in MAGMA (and LAPACK) the bug comes from dsidi because it still assumes upper triangular storage (so “U” case and making the matrix symmetric was OK, but not the “L” case). This is easy to fix and will be updated in the magma master shortly.

[Cosmin Petra]

interesting, thanks!

[Nai-Yuan Chiang]

@Stan Thank a lot for your quick fix! Cheers!

[Stan Tomov]

Hi Cosmin,

We added uplo flag to dsidi and the dsiinertia routines in MAGMA master to support symmetric matrices given in the lower triangular part (see commits f7053d1, a9a6fd7, 4098c17, and 9d9db1e). Before support was just for upper case through the original dsidi for CPUs, and similarly in dsiinertia for GPUs. As discussed, the problem was happening when dsytrf is called with “L” option, followed by dsidi (or dsiinertia) that were not supporting “L” option.

Please check if this fixes the issue.

Thanks,
Stan

[Nai-Yuan Chiang]

Hi Stan,

I tested the latest master branch with some examples and they all work fine.

Thanks a lot for your help!

Cheers,
NY

‌

‌

‌

[Stan Tomov]

Hi NY,

That’s great to hear everything works now!
Thank you for checking and testing more.
In this case I will close the issue shortly.

Thank you also again for discovering and reporting this bug.
Best regards,
Stan

[Stan Tomov]

Resolved by extending the functionality and support for symmetric matrices stored in the lower triangular part, and changes tested and verified by NY Chiang.

[Nai-Yuan Chiang]

Hi Stan,

Do you have an approximate time for the next release, in which I assume this fix will be included?

Cheers,

NY

[Stan Tomov]

Hi NY,

We plan next release to be around SC21, i.e., mid November.

Best regards,

Stan

[Nai-Yuan Chiang]

Thank you!

NY

[Cosmin Petra]

thank you Stan!