[Documentation]: Wrong decomposition in rocsolver_sgeblttrf_npvt

[amontoison]

Description of errors

The LU decomposition presented in the documentation is only valid for blocks of size 1x1.
Otherwise it is impossible to have L_1 I = B_1.

I tried to interface it in the Julia interface AMDGPU.jl and unable to recover the correct results with a block size > 2.
JuliaGPU/AMDGPU.jl#746

[EdDAzevedo]

Perhaps the notation of using "L_k" on the diagonal blocks is misleading since "L_k" is not lower triangular, but it holds the LU factorization ( "(L-eye) + U" ) of the (updated) diagonal blocks.

Maybe a simple fix is to use a different notation, say, "E_k" instead of "L_k" on the main diagonal blocks.

[amontoison]

Thanks for your answer @EdDAzevedo but I still don't understand how you recover E_1 * I = B_1. In all my tests B_1 != E_1 (except for blocks of size 1).

[EdDAzevedo]

I think the convention for LU factorization is something like L * U = B, where lower triangular matrix L has 1's on the main diagonal, U is upper triangular. The matrix B will be over-written by "L + U - eye", where "eye" is the identity matrix.

Since L has 1's on the main diagonal, there is no need to store those values.

So take for example a 2x2 matrix

B = [ B11 B12
B21 B22 ]
we want

L * U = B
or

[1 0 ] * [ U11 U12 ] = [B11 B12 ]
[ L21 1 ] [ 0 U22 ] [B21 B22 ]

so by expanding the above, we have

U11 = B11
L21 * U11 = B21, => L21 = B21/U11
U12 = B12
L21 * U12 + U22 = B22, => U22 = (B22 - L21 * U12 )

thus we can compute the 2x2 LU factorization by the above expressions.

By the convention, matrix B will be over-written as "L + U - eye" or as

[ U11 U12 ]
[ L21 U22 ]

A special case of B is a scalar or 1x1 matrix is "L" is just scalar 1 (so not stored) , U11 is just scalar B, so no change.

Note that knowing the L U factorization, we can form L and U (extracted from the over-written matrix) then multiply together to recover the original matrix B.

Storing the LU factorization is convenient when we need to perform the forward block triangular solve.

[amontoison]

I was probably not clear but my issue is not the storage, I just wanted to say that by products of the blocks in your documentation. You must have L_1 * I = B_1 and never found a way to assert it numerically.
I did a copy before B_1 is overwritted but never found valid factors.
Just by doing the matrix-matrix products between L and U by blocks.
I tried different storage that you could have used (block Doolittle or Crout versions) and none of them give the correct result except for a block size 1.

Can you confirm that you can recover a correct result numerically with cases of block size > 1?

I may do a C++ reproducer if it works on your side.

[amontoison]

Off-topic: I wanted to compare this solver with what we developed for SPD block tridiagonal solvers on GPU but I never recovered the correct result with the routines of rocSOLVER.
-> https://www.researchgate.net/publication/395206381_Harnessing_Batched_BLASLAPACK_Kernels_on_GPUs_for_Parallel_Solutions_of_Block_Tridiagonal_Systems

[EdDAzevedo]

What is the number of blocks "n" and block size "nb" and number of rhs vectors "nrhs"? If you need the batched version to solve multiple independent linear systems, what is the batch count?

[EdDAzevedo]

I suspect the rocsolver algorithm is the straight-forward thomas algorithm (not odd/even reduction). Moreover, it is not assuming symmetry. However, if you are not obtaining the correct results using rocsolver, would you kindly consider providing a reproducer or provide the problem sizes so we can investigate the issue?

Thanks.

[assistant-librarian]

This issue has been migrated to: ROCm/rocm-libraries#1681

[ammallya]

Imported to ROCm/rocm-libraries