Oh right, I just realized in the man that dppsv very likely decomposes
its A argument - instead of requiring a decomposed mat as I first
thought... So I was actually performing two successive Cholesky
decompositions ^^
On Sat, Dec 20, 2014 at 10:57 PM, peter dalgaard <pda...@gmail.com> wrote:
> This isn't the help list for LAPACK, but as far as I can tell, dppsv expects
> a symmetric matrix input compacted as triangular, not a Choleski decomposed
> one. So try assigning lmat before the call to dpotrf.
>
> -pd
>
>
>> On 20 Dec 2014, at 22:06 , Pierrick Bruneau <pbrun...@gmail.com> wrote:
>>
>> Dear R contributors,
>>
>> Considering the following sample C code, that illustrates two possible
>> uses of a Cholesky decomp for inverting a matrix, equally valid at
>> least in theory:
>>
>> SEXP test() {
>>
>> int d = 2;
>> int info = 0;
>> double mat[4] = {2.5, 0.4, 0.4, 1.7};
>> double id[4] = {1.0, 0.0, 0.0, 1.0};
>> double lmat[3];
>> F77_CALL(dpotrf)("L", &d, mat, &d, &info);
>> lmat[0] = mat[0];
>> lmat[1] = mat[1];
>> lmat[2] = mat[3];
>> F77_CALL(dppsv)("L", &d, &d, lmat, id, &d, &info);
>> // id now contains L^(-T)
>> F77_CALL(dpotri)("L", &d, mat, &d, &info);
>> // mat contains mat^(-1)
>>
>> Rprintf("%f\n", id[0] * id[0]);
>> // owing to that id is lower triangular
>> Rprintf("%f\n", mat[0]);
>>
>> return(R_NilValue);
>>
>> }
>>
>> I expected both printed values to be identical, or almost so. But
>> issuing .Call("test") prints:
>> 0.426571
>> 0.415648
>>
>> Difference is thus many degrees of magnitude above numerical
>> precision. What am I missing that explains it?
>>
>> Thanks by advance for the kind answers,
>> Pierrick
>>
>> ______________________________________________
>> R-devel@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>
> --
> Peter Dalgaard, Professor,
> Center for Statistics, Copenhagen Business School
> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> Phone: (+45)38153501
> Email: pd....@cbs.dk Priv: pda...@gmail.com
>
>
>
>
>
>
>
>
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
