[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
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [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 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
Re: [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 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 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
[Rd] loadNamespace and versionChecking and the otherpackage::otherfun syntax
This is an enquiry not so much about what the code for loadNamespace does,
but
rather about the intent and design of loadNamespace, and how it interacts
with
the `::` function, which seems to me to follow a slightly different
philosophy.
It is not an urgent question - the issue that started me wondering has been
resolved another way - but I would like to complete my understanding of this
aspect of how R's packaging mechanism is meant to operate.
It's also rather a long query - if it's too long please don't waste your
time
- just ignore it. To try and make it slightly more digestible
it is divided into sections, as follows
SCENARIO
THE QUERY
MY OWN ATTEMPTS TO UNDERSTAND
AN ASIDE ABOUT loadNamespace AND Depends:
VERSION CHECKS ON otherpackage::otherfun AT LOAD TIME?
VERSION CHECKS ON otherpackage::otherfun AT EXECUTION TIME?
POSSIBLE ANSWER 1 - THIS IS TOO COMPLEX AND HYPOTHETICAL
POSSIBLE ANSWER 2 - THIS IS AN ISSUE FOR THE `::` FUNCTION
POSSIBLE ANSWER 3 - loadNamespace SHOULD BE VERSION AWARE EVEN WHEN
`::` IS USED
Many thanks in advance for any insights that are able to be offered.
Geoff
SCENARIO
The scenario is that `mypackage` uses an `otherpackage` via the the
`otherpackage::otherfun` syntax, and that the version of `otherpackage` must
be
(say) (>= 2.0). The constraint on the version of `anotherpackage` is
specified
in the DESCRIPTION file of mypackage using either Imports or Depends.
THE QUERY
=
The query is :
a) Is it intended that loadNamespace should check the version of
otherpackage
(if one is specified) when it is `loadNamespace`ing mypackage?
and if so
b) Is it intended that the process of loading mypackage should ensure that
the
correct version of otherfun cannot be accidentally masked (for example by
using
.libPaths() to change the search path between the time when mypackage is
loaded,
and the time when the component of mypackage that calls
otherpackage::otherfun
is executed) ?
MY OWN ATTEMPTS TO UNDERSTAND
=
What I've done so far
I've read the documentation I could find, stepped through loadNamespace
using
debugonce (several times) with toy packages to gains some insight into what
loadNamespace does at the moment (using R3.1.2 on a Windows 64 bit machine),
and read the loadNamespace code a few times (though I can't yet claim to
follow
all of its neat tricks and complexities).
My understanding ths far is that
a) loadNamespace learns about version constraints on `otherpackage`
dependencies
when it is processing the DESCRIPTION file, viz
vI <- pkgInfo$Imports
and
b) defers the checking of these dependencies till later, when it is
processing
the NAMESPACE file (to create the imports:mypackage
namespace/environment/frame which
encloses the namespace:mypackage namespace/environment/frame). The checking
occurs
inside 3 loops, which all use an appropriate entry from vI as the
versionCheck
argument in a recursive call to loadNamespace, viz.
for (i in nsInfo$imports) {
...etc...
}
for (imp in nsInfo$importClasses) ...etc...
for (imp in nsInfo$importMethods) ...etc...
AN ASIDE ABOUT loadNamespace AND Depends:
=
As an aside it appears that any version dependencies specified in the
Depends
field are overlooked - I suspect that loadNamespace might be more complete
if
there were something like
vD <- pkgInfo$Depends
vID <- c(vD, vI) #possibly with an unlist thrown in somewhere?
and the versionCheck used vID instead of just vI
VERSION CHECKS ON otherpackage::otherfun AT LOAD TIME?
===
The problem with this elegant recursive approach to checking the version of
depended upon packages is that (as far as I know thus far) using the
implicit
loading syntax otherpackage::otherfun does not involve any entry in the
mypackage NAMESPACE file, and hence loadNamespace never checks the version
specification for otherpackage.
Hence the first part of my query - is it intended that loadNamespace perform
such a check?
My initial thought was that the answer should be yes, loadNamespace should
do
such a check, and so I wrote a little function which checked all the
versions
specified in the DESCRIPTION file, at the time the file was initially
encountered. (After a bit of debugging) it seemed to do what I wanted, in
that
if the right version of otherpackage was not available, my amended
loadNamespace
threw an error.
But then I started to think about how I could test it thoroughly - by which
I
mean not does my code do what I think it should do, but does it achieve the
outcome that motivated that coding in the first place.
VERSION CHECKS ON otherpackage::otherfun AT EXECUTION TIME?
That led to the second part of my query - should / how could loadNamespace
ensure that I actually get the otherfun from the ver
