LAPACK/ScaLAPACK Development

Hi,

I have implemented my own piece of code for computing a cholesky decomposition for symmetric packed matrices. To test it, I compared it with the LAPACK-routine dpptrf, but unfortunately the results were not similar.
However, compared to the dpotrf-routine, the results are exactly the same.
Besides the use of packed format in dpptrf is there any other algorithmic difference with dpotrf?
It seems that the second diagonal element (L(2,2)) is already quite different. In dpptrf it is just the square root of the original A(2,2)-element while it should be: sqrt(A(2,2)-A(2,1)²/A(1,1)).

Is there some sort of block algorithm used in dpptrf?

Best regards,
Arne

edit: I must say that the matrix was a singular matrix, does that make any difference?

Hi Arne,

Your formula for A(2,2) is correct. Are you sure you are storing the matrix in Lower Packed format as LAPACK expected it to be? See: http://www.netlib.org/lapack/lug/node123.html
I can have a deeper loop at the problem if any, but can you please check first that you use the good data layout.

Also, this might be of interest to you: LAPACK has another packed storage Cholesky factorization: DPFTRF (the name of the data storage is rectangular packed).
See: Fred G. Gustavson, Jerzy Waśniewski, Jack J. Dongarra, and Julien Langou. Rectangular full packed format for Cholesky's algorithm: Factorization, solution and inversion.
ACM Trans. Math. Software, 37(2):1-21, 2010.

Also, there are some routines to help convert matrices in various data storage, for example to go from full to packed this is DTRTTP, to go from full to rectangular packed, this is DTRTTF, etc.
(Note: Since LAPACK store only the triangular part of symmetric matrices, the triangular data layout conversion subroutines can be used there.)

Cheers,
Julien.

Hi Julien,

I use the C-interface of Lapack (LAPACKE) and it seems that DTRTTP is not working correctly (or I make a mistake in using it).
I have stored my symmetric matrix in full row-wise format (since it is symmetric row- or column-wise doesn't really matter).
When I use DTRTTP I do not get the correct matrix stored in packed format. The formula I use is:

info=LAPACKE_dtrttp(LAPACK_ROW_MAJOR,'L',n,fullA,n,symmA);

where fullA is the n x n matrix stored in full format, and symmA should be the lower triangular part of fullA stored in packed row-wise format.
So the third element of symmA should be the second diagonal element of fullA (fullA(2,2)). In stead it returns the element on position (1,4). This is even not the correct value if it was stored column wise.
However when using LAPACK_COL_MAJOR, symmA is entirely correct.
Perhaps this is also the cause of the incorrect functioning of dpptrf.

I tracked it down a bit further and according to me there's an error in the function dtp_trans.

Code: Select all

if( ( colmaj || upper ) && !( colmaj && upper ) ) {
        for( j = st; j < n; j++ ) {
            for( i = 0; i < j+1-st; i++ ) {
                out[ j-i + (i*(2*n-i+1))/2 ] = in[ ((j+1)*j)/2 + i ];
            }
        }
    } else {
        for( j = 0; j < n-st; j++ ) {
            for( i = j+st; i < n; i++ ) {
                out[ j + ((i+1)*i)/2 ] = in[ (j*(2*n-j+1))/2 + i-j ];
            }
        }
    }


When converting a column-major lower packed matrix, the second formula should be used (and not the first one), the same accounts for row-major upper packed matrix. It thus seems the formulas in the if-statement should be switched.

best regards,
Arne