Pierrick Bruneau
2014-Dec-20 21:06 UTC
[Rd] Unexplained difference between results of dppsv and dpotri LAPACK routines
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
peter dalgaard
2014-Dec-20 21:57 UTC
[Rd] Unexplained difference between results of dppsv and dpotri LAPACK routines
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
Pierrick Bruneau
2014-Dec-20 22:09 UTC
[Rd] Unexplained difference between results of dppsv and dpotri LAPACK routines
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 <pdalgd at 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 <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 > > > > > > > >