POSV and GETRS

[Danesh_D]

Hi all,

I have a symmetric matrix "A" and a non-symmetric matrix "B". I want to solve "AX=B" where "X" will be "X=(A^-1)*B". Now, when I use:

PDPOSV('L', A, B)

and

GETRF(A, pvtA)
GETRS('N', A, pvtA, B)

I expect that "B" which holds the value of "X" in "AX=B" should be same, but I get two different values for "B". Why "PDPOSV" doesn't return same values as "GETRS" while "A" is symmetric?

Regards,

Danesh

[Julien Langou]

Can you check the INFO message from DPOSV?

For DPOSV to work, A needs to Symmetric, yes, but as important it needs to be Positive
Definite. If A is Symmetric but not Positive Definite then DPOSV will fail, this is expected.
The INFO flag will not be ZERO. The initial matrix is destroyed but the right-hand side B
is still B in output (an not X).

(DPOSV tries to perform a Cholesky factorization of your initial matrix A and this is
possible if and only if A is Symmetric Positive Definite.)

So, basically, if all you know is that your matrix is Symmetric, and want to use DPOSV,
you really need to check the INFO flag to be sure the answer is correct.

In the other hand, you can try DSYSV, this one is for any Symmetric matrix.

Julien

[Danesh_D]

I can not find "DSYSV" in ScaLAPACK and PBLAS and thus can not write wrapper for it in C. "POSV" was already in ScaLAPACK but not "DSYSV". Am I right?

[Julien Langou]

Hello,

sorry I forgot that you were using ScaLAPACK. In that case you are right. There is
currently no SYSV in ScaLAPACK (while there is one in LAPACK). So if you have
a symmetric indefinite matrix, the best you can do is to call GESV. (This doubles
the FLOPS count and double the memory.)

SYSV is one important feature that is missing in ScaLAPACK, we know ....
Best is to send request email or forum post to let us know. At some point,
some one will implement the routine.

Julien.

[Danesh_D]

Hi Julien,

I checked returned value from "PDPOSV" and it was zero which means that the operation should be done successfully, but when I check the result of "PDPOSV" with "PDGETRF+PDGETRS" they are different with same inputs! I also checked the resuls of "PDPOTRF+PDPOTRS" and the results are as same as "PDPOSV". It is very weired to me. Why these symmetric functions return difference values than "PDGETRF+PDGETRS"?

Danesh

[Julien Langou]

Hello Danesh,
so in this case, there is probably a bug in your code.
Standard errors are:
- do not provide the same input (is b the same as PDPOSV and for PDGESV)?
- do not provide the same matrix
- mess up in the transpose
... good luck ...
a good idea is to check the residual (b-Ax) in both case
Julien.

[Danesh_D]

Thanks Julien.

Well, the only difference is that "A" matrix in "POSV" and "POTRF+POTRF" is filled only in this lower triangle and these functions are also called with "L" as "uplo" variable to reference it to lower triangle, but "GETRF+GETRS" uses complete "A" matrix which has both upper and lower triangles filled with data and they are same since "A" is symmetric. Should I uses full filled "A" for symmetric functions too, however I reference only "A"'s lower trinagle?

Regards,

Danesh

[Julien Langou]

Hello Danesh,

can you try symmetrize your matrix using PDGEADD ?
I can not believe there is a problem in PDGESV or PDPOSV.

Julien.

[Danesh_D]

Actually I did and when my matrix is symmetric, for example POSV works fine. But, when I only fill lower triangle (to save time) and thus my matrix is not symmetric anymore and only reference lower triangle by calling POSV, the result changes from when my matrix was symmetric. I mean should my matrix be symmetric even if I reference lower triangle of it when I call POSV?

Thanks,

Danesh

[Julien Langou]

I think you are messing up with what is a matrix for ScaLAPACK.
Here is an example, consider the 7-by-7 matrix

Code: Select all

( 0 1 2 3 4 5 6 )
( 1 7 8 9 A B C )
( 2 8 D E F G H )
( 3 9 E I J K L )
( 4 A F J M N O )
( 5 B G K N P Q )
( 6 C H L O Q R )


we are going to map it on a 2-by-2 process grid with 2-by-2 blocks, this gives
the following blocking

Code: Select all

( 0 1 | 2 3 | 4 5 | 6 )
( 1 7 | 8 9 | A B | C )
-----------------------
( 2 8 | D E | F G | H )
( 3 9 | E I | J K | L )
-----------------------
( 4 A | F J | M N | O )
( 5 B | G K | N P | Q )
-----------------------
( 6 C | H L | O Q | R )


and get partitionned as:

Code: Select all

Process (0,0)    Process (0,1)    Process (1,0)    Process (1,1)
( 0 1 4 5 )      ( 2 8 F G )     ( 2 3 6 )         ( D E H )
( 1 7 A B )      ( 3 9 J K )     ( 8 9 C )         ( E I L )
( 4 A M N )      ( 6 C O Q )     ( F J O )         ( H L R )
( 5 B N P )                      ( G K Q )


So that when you only want to reference the lower part you end up with

Code: Select all

Process (0,0)    Process (0,1)    Process (1,0)    Process (1,1)
( 0       )      ( 2 8     )     (       )         ( D     )
( 1 7     )      ( 3 9     )     (       )         ( E I   )
( 4 A M   )      ( 6 C O Q )     ( F J   )         ( H L R )
( 5 B N P )                      ( G K   )


We agree on this?

Start again with your full symmetric matrix, give a try to PDLASET, to set ZERO in
the upper part, then try PDPOSV (this will work) and check what form the input
matrix as.

Good luck,
Julien.