[Rd] Unexplained difference between results of dppsv and dpotri LAPACK routines

[peter dalgaard]

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 <pbruneau at 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 at 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.mes at cbs.dk  Priv: PDalgd at gmail.com