Use in-place left division (#196)

Author: junyuan-chen

src/fit.jl

@@ -298,15 +298,15 @@ function reg(
 
         # Build
         newZ = hcat(Xexo, Z)
-        Pi = cholesky!(Symmetric(newZ' * newZ)) \ (newZ' * Xendo)
+        Pi = ldiv!(cholesky!(Symmetric(newZ'newZ)), newZ'Xendo)
         Xhat = hcat(Xexo, newZ * Pi)
         X = hcat(Xexo, Xendo)
 
         # prepare residuals used for first stage F statistic
         ## partial out Xendo in place wrt (Xexo, Z)
         Xendo_res = BLAS.gemm!('N', 'N', -1.0, newZ, Pi, 1.0, Xendo)
         ## partial out Z in place wrt Xexo
-        Pi2 = cholesky!(Symmetric(Xexo' * Xexo)) \ (Xexo' * Z)
+        Pi2 = ldiv!(cholesky!(Symmetric(Xexo'Xexo)), Xexo'Z)
         Z_res = BLAS.gemm!('N', 'N', -1.0, Xexo, Pi2, 1.0, Z)
     else
         # get linearly independent columns
@@ -324,8 +324,8 @@ function reg(
     ##
     ##############################################################################
 
-    crossx = cholesky!(Symmetric(Xhat' * Xhat))
-    coef = crossx \ (Xhat' * y)
+    crossx = cholesky!(Symmetric(Xhat'Xhat))
+    coef = ldiv!(crossx, Xhat'y)
 
     ##############################################################################
     ##

src/partial_out.jl

@@ -121,10 +121,9 @@ function partial_out(
     end
     # Compute residuals
     if size(X, 2) > 0
-        residuals = Y .- X * (X \ Y)
-    else
-        residuals = Y
+        mul!(Y, X, X\Y, -1.0, 1.0)
     end
+    residuals = Y
 
     # rescale residuals
     if has_weights