[Genabel-commits] r1638 - branches/ProbABEL-0.50/src

[noreply at r-forge.r-project.org]

Author: maartenk
Date: 2014-03-11 23:12:35 +0100 (Tue, 11 Mar 2014)
New Revision: 1638

Modified:
   branches/ProbABEL-0.50/src/reg1.cpp
Log:
MMSCORE: same method for multiplying two predictors as introduced in r1636 for 1 predictor. The methods for 1 predictor is more than 2 to 2.5 faster than method used before. This method is faster because the line 
tXW= X.transpose()*W caused X to have not enough rows to let vectorizing work efficient. When using W*X the left-hand side has enough rows for vectorization. 

Modified: branches/ProbABEL-0.50/src/reg1.cpp
===================================================================
--- branches/ProbABEL-0.50/src/reg1.cpp	2014-03-09 21:47:13 UTC (rev 1637)
+++ branches/ProbABEL-0.50/src/reg1.cpp	2014-03-11 22:12:35 UTC (rev 1638)
@@ -317,11 +317,11 @@
     VectorXd Y = reg_data.Y.data.col(0);
     if (X.data.cols() == 3)
     {
-        Matrix<double, 3, Dynamic> tXW = X.data.transpose()* W_masked.masked_data->data;
-        Matrix3d xWx = tXW * X.data;
+        Matrix<double, Dynamic, 3> tXW = W_masked.masked_data->data * X.data;
+        Matrix2d xWx = tXW.transpose() * X.data;
         Ch = LDLT<MatrixXd>(xWx);
-        Vector3d beta_3f = Ch.solve(tXW * Y);
-        sigma2 = (Y - tXW.transpose() * beta_3f).squaredNorm();
+        Vector3d beta_3f = Ch.solve(tXW.transpose() * Y);
+        sigma2 = (Y - tXW * beta_3f).squaredNorm();
         beta.data = beta_3f;
     }
     else if (X.data.cols() == 2)