- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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.

Link Copied

9 Replies

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

*> 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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

*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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

[bash]! TAU = TAU_I^{-1}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.

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]

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Thanks,

D.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

*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.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page