<div dir="ltr">Deiego, what is the "p" (n^2 -> np) in the Log message?</div><div class="gmail_extra"><br><br><div class="gmail_quote">On Fri, Feb 7, 2014 at 6:32 PM,  <span dir="ltr"><<a href="mailto:noreply@r-forge.r-project.org" target="_blank">noreply@r-forge.r-project.org</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Author: dfabregat<br>
Date: 2014-02-07 18:32:24 +0100 (Fri, 07 Feb 2014)<br>
New Revision: 1600<br>
<br>
Modified:<br>
   pkg/OmicABEL/src/REML.c<br>
Log:<br>
Performance improvement for REML estimation.<br>
Reusing precomputed data to replace an expensive<br>
n^2 gemv for a cheaper n*p gemv.<br>
<br>
<br>
Modified: pkg/OmicABEL/src/REML.c<br>
===================================================================<br>
--- pkg/OmicABEL/src/REML.c     2014-02-06 21:29:08 UTC (rev 1599)<br>
+++ pkg/OmicABEL/src/REML.c     2014-02-07 17:32:24 UTC (rev 1600)<br>
@@ -194,10 +194,15 @@<br>
        // loglik = a + b<br>
        //  a -> log(det(M))<br>
        //  b -> YmXB' inv(M) YmXB<br>
-    dgemv_(TRANS,<br>
+    /*dgemv_(TRANS,<br>
             &n, &n,<br>
             &ONE, Z, &n, YmXB, &iONE,<br>
-            &ZERO, ZtY_upd, &iONE);<br>
+            &ZERO, ZtY_upd, &iONE);*/<br>
+    memcpy( ZtY_upd, ZtY, n * sizeof(double) );<br>
+    dgemv_( NO_TRANS,<br>
+            &n, &wXL,<br>
+            &MINUS_ONE, ZtX, &n, beta, &iONE,<br>
+            &ONE, ZtY_upd, &iONE );<br>
        // YmXB' * inv( M ) * YmXB<br>
        *loglik = 0.0;<br>
        for (i = 0; i < n; i++ )<br>
<br>
_______________________________________________<br>
Genabel-commits mailing list<br>
<a href="mailto:Genabel-commits@lists.r-forge.r-project.org">Genabel-commits@lists.r-forge.r-project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/genabel-commits" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/genabel-commits</a><br>
</blockquote></div><br><br clear="all"><div><br></div>-- <br>-----------------------------------------------------<br>Yurii S. Aulchenko<br><div><br></div><div>[ <a href="http://nl.linkedin.com/in/yuriiaulchenko" target="_blank">LinkedIn</a> ] [ <a href="http://twitter.com/YuriiAulchenko" target="_blank">Twitter</a> ] [ <a href="http://yurii-aulchenko.blogspot.nl/" target="_blank">Blog</a> ]</div>

</div>