[GenABEL-dev] [Genabel-commits] r1625 - branches/ProbABEL-0.50/src

[L.C. Karssen]

Hi Maarten,

On 27-02-14 22:56, noreply at r-forge.r-project.org wrote:
> Author: maartenk
> Date: 2014-02-27 22:56:38 +0100 (Thu, 27 Feb 2014)
> New Revision: 1625
> 
> Modified:
>    branches/ProbABEL-0.50/src/reg1.cpp
> Log:
> using predefined matrices in mmscore method. Initial testing showed a 35% speedup in the regression part with(in a separate toy program to test regression)

So you're saying that using predefined Matrix sizes speeds up the
computation by roughly a third? Wow! SO if I understand it correctly,
the allocation of an fixed/known sized matrix followed by matrix
multiplication is much faster than in the case that dimensions of the
matrix are not fixed... Right?

> 
> Modified: branches/ProbABEL-0.50/src/reg1.cpp
> ===================================================================
> --- branches/ProbABEL-0.50/src/reg1.cpp	2014-02-26 21:43:42 UTC (rev 1624)
> +++ branches/ProbABEL-0.50/src/reg1.cpp	2014-02-27 21:56:38 UTC (rev 1625)
> @@ -357,6 +357,7 @@
>      double sigma2_internal;
>  
>  #if EIGEN
> +
>      LDLT <MatrixXd> Ch;
>  #else
>      mematrix<double> tXX_i;
> @@ -373,12 +374,30 @@
>          //Oct 26, 2009
>  
>  #if EIGEN
> -        // next line is  5997000 flops
> -        MatrixXd tXW = X.data.transpose() * invvarmatrixin.masked_data->data;
> -        Ch = LDLT <MatrixXd>(tXW * X.data); // 17991 flops
> -        beta.data = Ch.solve(tXW * reg_data.Y.data);//5997 flops
> -        //next line is: 1000+5000+3000= 9000 flops
> -        sigma2 = (reg_data.Y.data - tXW.transpose() * beta.data).squaredNorm();
> +        cout << "BB"<<X.data.cols()<<"AAAAAAa"<<endl;
> +        if (X.data.cols()== 3){
> +            Matrix<double,3,Dynamic> tXW = X.data.transpose()*invvarmatrixin.masked_data->data;
> +            Matrix3d xWx = tXW * X.data;
> +            Ch = LDLT <MatrixXd>  (xWx );
> +            Vector3d beta_3f = Ch.solve(tXW * reg_data.Y.data);
> +            sigma2 = (reg_data.Y.data - tXW.transpose() * beta_3f).squaredNorm();
> +            beta.data = beta_3f;
	
	The two lines above are repeated almost identically in the cases below.
Wouldn't it make more sense to simple call the Vector{3d,2d,beta.data}
in each of the cases 'betaf'. Then you could take the sigma2 calculation
and assignment to beta.data out of the ifs. Like this:
         Vector3d betaf = Ch.solve(tXW * reg_data.Y.data);
      }
	... the else if and else here...
      }
         sigma2 = (reg_data.Y.data - tXW.transpose() *
                    betaf).squaredNorm();
         beta.data = betaf;
  #else


Thanks!

Lennart.


> +        }
> +       else if(X.data.cols()== 2){
> +            Matrix<double,2,Dynamic> tXW = X.data.transpose()*invvarmatrixin.masked_data->data;
> +            Matrix2d xWx = tXW * X.data;
> +            Ch = LDLT <MatrixXd>  (xWx );
> +            Vector2d beta_2f = Ch.solve(tXW * reg_data.Y.data);
> +            sigma2 = (reg_data.Y.data - tXW.transpose() * beta_2f).squaredNorm();
> +            beta.data = beta_2f;
> +        }else{
> +            // next line is  5997000 flops
> +            MatrixXd tXW = X.data.transpose() * invvarmatrixin.masked_data->data;
> +            Ch = LDLT <MatrixXd>(tXW * X.data); // 17991 flops
> +            beta.data = Ch.solve(tXW * reg_data.Y.data);//5997 flops
> +            //next line is: 1000+5000+3000= 9000 flops
> +            sigma2 = (reg_data.Y.data - tXW.transpose() * beta.data).squaredNorm();
> +        }
>  #else
>          // next line is  5997000 flops
>          mematrix<double> tXW = transpose(X) * invvarmatrixin.masked_data;
> 
> _______________________________________________
> Genabel-commits mailing list
> Genabel-commits at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/genabel-commits
> 

-- 
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
L.C. Karssen
Utrecht
The Netherlands

lennart at karssen.org
http://blog.karssen.org
GPG key ID: A88F554A
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 230 bytes
Desc: OpenPGP digital signature
URL: <http://lists.r-forge.r-project.org/pipermail/genabel-devel/attachments/20140227/120e565b/attachment.sig>