pdpotrf() fails to identify SPD matrix

[Georgios_S_]

I am trying to perform a Cholesky decomposition via pdpotrf() of MKL-Intel's library, which uses ScaLAPACK. I am reading the whole matrix in the master node and then distribute it like in this example. Everything works fine when the dimension of the SPD matrix is even. However, when it's odd, `pdpotrf()` thinks that the matrix is not positive definite.

Could it be because the submatrices are not SPD? I am working with this matrix and the submatrices are (with 4 processes and blocks of size 2x2):

    A_loc on node 0
      4   1   2
      1 0.5   0
      2   0  16
    
    nrows = 3, ncols = 2
    A_loc on node 1
      2 0.5
      0   0
      0   0
    
    nrows = 2, ncols = 3
    A_loc on node 2
      2   0   0
    0.5   0   0
    
    nrows = 2, ncols = 2
    A_loc on node 3
      3   0
      0 0.625

Here, every submatrix is not SPD, however, the overall matrix is SPD (have checked with running with 1 process). What should I do? Or there is nothing I can do and `pdpotrf()` does not work with matrices of odd size?

[Ying_H_Intel]

Hi Georgios, 

Not sure how  the input parameters  look like.  Could you please attach the cpp code so we can see what is the problem? 

From document, 

(local)

REAL for pspotrf

DOUBLE PRECISON for pdpotrf

COMPLEX for pcpotrf

DOUBLE COMPLEX for pzpotrf.

Pointer into the local memory to an array of dimension (lld_a, LOCc(ja+n-1)).

On entry, this array contains the local pieces of the n-by-n symmetric/Hermitian distributed matrix sub(A) to be factored.

Depending on uplo, the array a contains either the upper or the lower triangular part of the matrix sub(A) (see uplo).  Not the entire SPD matrix. 

Best Regards,
Ying 
 

[Georgios_S_]

Hi Ying,

  I tried to do so, by editing my first post. However, it required a captcha code to verify that I am a human being. I think I made a mistake at my first attempt. When I tried again for about 7 times, I am pretty sure I typed the correct symbols. However, I triggered the spam filter. I did sent a feedback, hope this will help.

  I did make a minor example, which behaves similarly to my real code. I have posted it at Stackoverflow. If needed, I can post the code here as well.

Thanks,

George

[Ying_H_Intel · Accepted solution]

Hi George, 

I modify the code and  add the call of pdoptrf. It looks similiar as yours, please try it and let us know if any findings. 

Best Regards,

Ying

Here is the compile command and run result. 

[yhu5@prc-mic01 scalapack]$ mpiicc -o pdpotrf pdpotrf.cpp -I. -I/opt/intel/mkl/include  /opt/intel/mkl/lib/intel64/libmkl_scalapack_lp64.a -L/opt/intel/mkl/lib/intel64 -Wl,--start-group /opt/intel/mkl/lib/intel64/libmkl_intel_lp64.a /opt/intel/mkl/lib/intel64/libmkl_core.a /opt/intel/mkl/lib/intel64/libmkl_sequential.a /opt/intel/mkl/lib/intel64/libmkl_blacs_intelmpi_lp64.a -Wl,--end-group -lpthread -lm -ldl
[yhu5@prc-mic01 scalapack]$ mpirun -n 4 ./pdpotrf a.txt 5 5 2 2
Matrix A:
  4   1   2 0.5   2
  1 0.5   0   0   0
  2   0   3   0   0
0.5   0   0 0.625   0
  2   0   0   0  16

Processes grid pattern:
0 1
2 3
A_loc on node 0
  4   1   2
  1 0.5   0
  2   0  16

A_loc on node 1
  2 0.5
  0   0
  0   0

A_loc on node 2
  2   0   0
0.5   0   0

A_loc on node 3
  3   0
  0 0.625

pdpotrf 00
pdpotrf 20
pdpotrf 10
pdpotrf 30
Matrix A test:
  2   1   2 0.5   2
0.5 0.5   0   0   0
  1  -1   1   0   0
0.25 -0.25 -0.5 0.5   0
  1  -1  -2  -3   1

[Georgios_S_]

Hi Ying,

 this seems to be OK. I commented #include <pblas_examples.h>, since it produced a fattal error and moreover seems non-needed. Moreover, I can see that you are based on the example, rather than the minor example, which is OK. However, the example has the values of procrows and proccols hardcoded. I did this dynamic, by letting the user set a value, just like for N, M, Nb and Mb. Is that OK, or you think I am on the wrong track? Recall that procrows * proccols = number of processes (they are the dimensions of the grid of processes in other words).

  As stated here, the problem came from this line:

MPI_Bcast(&lda, 1, MPI_INT, 0, MPI_COMM_WORLD);

, where lda is the leading dimension of the local matrix. I thought lda was the leading dimension of the global matrix.

  Oh, also notice that in your example, you use two global matrices. I had done some experiments and did check that it will work fine with one global matrix as well. So, is it just the 2 global matrices leftovers from the example you were based or is there any particular reason you did so?

Sorry for my wrong understanding of the parameter,

George

[Ying_H_Intel]

Hi George, 

You are right, the problem in your code is MPI_Bcast(&lda, 1, MPI_INT, 0, MPI_COMM_WORLD);  After remove, you will see it works. 

I just quick modify the code in example.  So keep the original code, the 2 Global matrix is not needed, just for test.  An d  Yes, it is ok for you to  letting the user set a value for all kind of parameters.  

Best Regards,
Ying