matrix inversion using packed storage

[Danesh_Daroui]

Hi,

I have a real symmetric matrix that I want to calculate inverse of it. For testing, I first create the full matrix and then save it in packed storage to see if everything is correct. Then I call packed storage routines for choleskt factorization (I have symmetric matrix) and then matrix inversion. The code I use is:

[cpp]for (int i=0; i

where the matrix I am trying to invert is a sub-matrix "nx-by-nx" which is in the bigger matrix "nCells-by-nCells". The code is clear enough how I read data from the larger matrix and store it as a packed sub-matrix. The manual of MKL suggests the addressing of the data in one-dimension vector when lower triangle is used as "ap(i+(2*n-j)*(j-1)/2)" so I just modified it to change to to zero-based addressing as above. I hope it is correct!

When I run the code, the routine "dpptrf" always return the number "22" in "info" variable. The strange thing is that the code works OK (info is always zero) for smaller matrices but as the dimensions grow the error comes up. Does anybody know where my code can be faulty? I have read it several time but everything seem to be correct. The manual says that positive value of "info" can be because of error formatting the matrix in packed storage format but if my formatting is incorrect why it does work for smaller matrices? Besides, I have checked the addressing many times and it seems to be correct.

Thanks in advance,

D.

[mecej4]

I do not see anything wrong with the packing of the triangular matrix or the call to Lapack. You used 'always' twice in the second paragraph in a misleading way.

It is possible that the matrix as passed is not positive definite for certain values of nx. The returned value of 22 suggests just that possibility, i.e., it is not the arrangement of the matrix elements that is causing the problem but the values of certain elements. You could check this by dumping the smallest matrix that elicits info=22 and trying the Cholesky factorization of that using, e.g., Matlab or Octave.

[Danesh_Daroui]

The matrix is not triangular but symmetric and rectangular. The problem is that using "dgetrf + dgetri" the sub-matrix can be inverted successfully but when I pack it and call Cholesky factorization, the routine complains as I wrote. If this is not positive definite, the factorization should fail when I use "dgetrf" too. Isn't it true?

[mecej4]

> The matrix is not triangular but symmetric and rectangular

Well, in mathematics all matrices are rectangular but "triangular matrix" is standard in numerical analysis as representing the lower or upper triangle including the main diagonal.

> ...If this is not positive definite, the factorization should fail when I use "dgetrf" too..

Not true. Routine DGETRF is for general matrices. and computes the A = P . L . U factorization. This factorization exists for symmetric indefinite and unsymmetric matrices, whereas the Cholesky factorization exists only for positive definite matrices.

Here is a simple matrix that DGETRF will factor but not DPOTRF:

A = [1 2; 2 1]

and the vector x = [1; -2] will prove that A is indefinite ( x^T . A . x = -3 ).

Your results show that for small values of nx the submatrix that you hand over to DPOTRF is positive definite but for larger values of nx it is not. You will have to investigate whether/why that is indeed so, or whether the values of the matrix elements are incorrect.

[Danesh_Daroui]

I am not agree about your definition of "triangular matrices"! A triangular matrix is a special kind of square matrix where all its elements below/above the diagonal are zero, but it is not the issue now.

Anyway, I checked more and found out that as you mentioned my matrix is symmetric and indefinite so I used Bunch-Kaufman factorization for packed storage and it worked correctly.

Thanks,

D.

[mecej4]

> A triangular matrix is a special kind of square matrix where all its elements below/above the diagonal are zero

I don't think there is a disagreement here. The one extra feature in numerical analysis is that the elements in the zero part of a triangular matrix are not even referenced in an algorithm dealing with such matrices. Not only are needless operations with zeroes avoided, the storage of the matrix consumes a little over half that of the full matrix.

[junuylia]

For symmetric matrix, I used packed format and the routines for symmetric matrices, and it worked pretty well, you might have a try.

[bash]!  TAU  = TAU_I^{-1}
          CALL COPY ( TAU_I, TAU)

!  CHOLESKY FACTORIZATION, THE MATRIX IS UPDATED BY THE LOWER MATRIX
          CALL PPTRF( TAU, UPLO, INFO)

!  INVERSION AFTER THE CHOLESKY FACTORIZATION, THE MATRIX IS OVERWRITTEN
!    BY THE INVERSION
          CALL PPTRI( TAU, UPLO, INFO)
[/bash]

I'm trying to get the inversion of tau_i, since the original matrix would be overwritten, I used a substitute. Also, make sure the triangular matrix is stored by column.

[yuriisig]

Unfortunately, algorithms of operation with the packed matrixes in Intel MKL are ineffective.

[Danesh_Daroui]

I guess PPTRF and PPTRI use Cholesky factorization, so as I wrote my matrices are not always positive definite and therefore I could not use Cholesky factorization.

Thanks,

D.

[Danesh_Daroui]

Is that true that MKL is inefficient for packed storage matrices? Why would it be? Packed storage only differs in addressing and should not affect the performance. However, I realized that the factorization works in parallel but matrix inversion is totally serial (only one core is used). Would it be same for other inversions like ?getri or this is only for packed storage?

Regards,

D.