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
Unexplained difference between results of dppsv and dpotri LAPACK routines
3 messages · Peter Dalgaard, Pierrick Bruneau
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
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