Author: Andrius Merkys <merkys@debian.org>
Description: Migrating to libejml-java v0.38.
--- a/src/dr/evomodel/treedatalikelihood/continuous/cdi/ContinuousDiffusionIntegrator.java
+++ b/src/dr/evomodel/treedatalikelihood/continuous/cdi/ContinuousDiffusionIntegrator.java
@@ -413,7 +413,7 @@
             final double vi = branchLengths[imo];
             final double vj = branchLengths[jmo];
 //
-//            final DenseMatrix64F Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait);
+//            final DMatrixRMaj Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait);
 //
             if (DEBUG) {
                 System.err.println("updatePreOrderPartial for node " + iBuffer);
@@ -427,33 +427,33 @@
                 final double pk = prePartials[kbo + dimTrait];
                 final double pj = partials[jbo + dimTrait];
 
-//                final DenseMatrix64F Pk = wrap(prePartials, kbo + dimTrait, dimTrait, dimTrait);
-//    //                final DenseMatrix64F Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait);
+//                final DMatrixRMaj Pk = wrap(prePartials, kbo + dimTrait, dimTrait, dimTrait);
+//    //                final DMatrixRMaj Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait);
 //
-//    //                final DenseMatrix64F Vk = wrap(prePartials, kbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
-//                final DenseMatrix64F Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+//    //                final DMatrixRMaj Vk = wrap(prePartials, kbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+//                final DMatrixRMaj Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
 //
 
                 // B. Inflate variance along sibling branch using matrix inversion
                 final double pjp = Double.isInfinite(pj) ?
                         1.0 / vj : pj / (1.0 + pj * vj);
 
-//    //                final DenseMatrix64F Vjp = new DenseMatrix64F(dimTrait, dimTrait);
-//                final DenseMatrix64F Vjp = matrix0;
-//                CommonOps.add(Vj, vj, Vd, Vjp);
+//    //                final DMatrixRMaj Vjp = new DMatrixRMaj(dimTrait, dimTrait);
+//                final DMatrixRMaj Vjp = matrix0;
+//                CommonOps_DDRM.add(Vj, vj, Vd, Vjp);
 //
-//    //                final DenseMatrix64F Pjp = new DenseMatrix64F(dimTrait, dimTrait);
-//                final DenseMatrix64F Pjp = matrix1;
+//    //                final DMatrixRMaj Pjp = new DMatrixRMaj(dimTrait, dimTrait);
+//                final DMatrixRMaj Pjp = matrix1;
 //                InversionResult cj = safeInvert(Vjp, Pjp, false);
 //
-//    //                final DenseMatrix64F Pip = new DenseMatrix64F(dimTrait, dimTrait);
-//                final DenseMatrix64F Pip = matrix2;
-//                CommonOps.add(Pk, Pjp, Pip);
+//    //                final DMatrixRMaj Pip = new DMatrixRMaj(dimTrait, dimTrait);
+//                final DMatrixRMaj Pip = matrix2;
+//                CommonOps_DDRM.add(Pk, Pjp, Pip);
 
                 final double pip = pjp + pk;
 //
-//    //                final DenseMatrix64F Vip = new DenseMatrix64F(dimTrait, dimTrait);
-//                final DenseMatrix64F Vip = matrix3;
+//    //                final DMatrixRMaj Vip = new DMatrixRMaj(dimTrait, dimTrait);
+//                final DMatrixRMaj Vip = matrix3;
 //                InversionResult cip = safeInvert(Pip, Vip, false);
 //
 //                // C. Compute prePartial mean
@@ -484,11 +484,11 @@
                 final double pi  = Double.isInfinite(pip) ?
                         1.0 / vi : pip / (1.0 + pip * vi);
 
-//                final DenseMatrix64F Vi = Vip;
-//                CommonOps.add(vi, Vd, Vip, Vi);
+//                final DMatrixRMaj Vi = Vip;
+//                CommonOps_DDRM.add(vi, Vd, Vip, Vi);
 //
-//    //                final DenseMatrix64F Pi = new DenseMatrix64F(dimTrait, dimTrait);
-//                final DenseMatrix64F Pi = matrix4;
+//    //                final DMatrixRMaj Pi = new DMatrixRMaj(dimTrait, dimTrait);
+//                final DMatrixRMaj Pi = matrix4;
 //                InversionResult ci = safeInvert(Vi, Pi, false);
 //
 //                // X. Store precision results for node
--- a/src/dr/evomodel/treedatalikelihood/continuous/cdi/MultivariateIntegrator.java
+++ b/src/dr/evomodel/treedatalikelihood/continuous/cdi/MultivariateIntegrator.java
@@ -2,8 +2,8 @@
 
 import dr.math.matrixAlgebra.WrappedVector;
 import dr.math.matrixAlgebra.missingData.InversionResult;
-import org.ejml.data.DenseMatrix64F;
-import org.ejml.ops.CommonOps;
+import org.ejml.data.DMatrixRMaj;
+import org.ejml.dense.row.CommonOps_DDRM;
 
 import java.util.HashMap;
 import java.util.Map;
@@ -56,13 +56,13 @@
 
     private final Map<String, Long> times;
 
-    DenseMatrix64F matrix0;
-    DenseMatrix64F matrix1;
-    DenseMatrix64F matrix2;
-    DenseMatrix64F matrix3;
-    DenseMatrix64F matrix4;
-    DenseMatrix64F matrix5;
-    DenseMatrix64F matrix6;
+    DMatrixRMaj matrix0;
+    DMatrixRMaj matrix1;
+    DMatrixRMaj matrix2;
+    DMatrixRMaj matrix3;
+    DMatrixRMaj matrix4;
+    DMatrixRMaj matrix5;
+    DMatrixRMaj matrix6;
 
     double[] vector0;
 
@@ -70,13 +70,13 @@
         inverseDiffusions = new double[dimTrait * dimTrait * diffusionCount];
 
         vector0 = new double[dimTrait];
-        matrix0 = new DenseMatrix64F(dimTrait, dimTrait);
-        matrix1 = new DenseMatrix64F(dimTrait, dimTrait);
-        matrix2 = new DenseMatrix64F(dimTrait, dimTrait);
-        matrix3 = new DenseMatrix64F(dimTrait, dimTrait);
-        matrix4 = new DenseMatrix64F(dimTrait, dimTrait);
-        matrix5 = new DenseMatrix64F(dimTrait, dimTrait);
-        matrix6 = new DenseMatrix64F(dimTrait, dimTrait);
+        matrix0 = new DMatrixRMaj(dimTrait, dimTrait);
+        matrix1 = new DMatrixRMaj(dimTrait, dimTrait);
+        matrix2 = new DMatrixRMaj(dimTrait, dimTrait);
+        matrix3 = new DMatrixRMaj(dimTrait, dimTrait);
+        matrix4 = new DMatrixRMaj(dimTrait, dimTrait);
+        matrix5 = new DMatrixRMaj(dimTrait, dimTrait);
+        matrix6 = new DMatrixRMaj(dimTrait, dimTrait);
     }
 
     @Override
@@ -86,9 +86,9 @@
         assert (inverseDiffusions != null);
 
         final int offset = dimTrait * dimTrait * precisionIndex;
-        DenseMatrix64F precision = wrap(diffusions, offset, dimTrait, dimTrait);
-        DenseMatrix64F variance = new DenseMatrix64F(dimTrait, dimTrait);
-        CommonOps.invert(precision, variance);
+        DMatrixRMaj precision = wrap(diffusions, offset, dimTrait, dimTrait);
+        DMatrixRMaj variance = new DMatrixRMaj(dimTrait, dimTrait);
+        CommonOps_DDRM.invert(precision, variance);
         unwrap(variance, inverseDiffusions, offset);
 
         if (DEBUG) {
@@ -124,7 +124,7 @@
         final double vi = branchLengths[imo];
         final double vj = branchLengths[jmo];
 
-        final DenseMatrix64F Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait);
+        final DMatrixRMaj Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait);
 
         if (DEBUG) {
             System.err.println("updatePreOrderPartial for node " + iBuffer);
@@ -137,27 +137,27 @@
         for (int trait = 0; trait < numTraits; ++trait) {
 
             // A. Get current precision of k and j
-            final DenseMatrix64F Pk = wrap(prePartials, kbo + dimTrait, dimTrait, dimTrait);
-//                final DenseMatrix64F Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait);
+            final DMatrixRMaj Pk = wrap(prePartials, kbo + dimTrait, dimTrait, dimTrait);
+//                final DMatrixRMaj Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait);
 
-//                final DenseMatrix64F Vk = wrap(prePartials, kbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
-            final DenseMatrix64F Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+//                final DMatrixRMaj Vk = wrap(prePartials, kbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+            final DMatrixRMaj Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
 
             // B. Inflate variance along sibling branch using matrix inversion
-//                final DenseMatrix64F Vjp = new DenseMatrix64F(dimTrait, dimTrait);
-            final DenseMatrix64F Vjp = matrix0;
-            CommonOps.add(Vj, vj, Vd, Vjp);
+//                final DMatrixRMaj Vjp = new DMatrixRMaj(dimTrait, dimTrait);
+            final DMatrixRMaj Vjp = matrix0;
+            CommonOps_DDRM.add(Vj, vj, Vd, Vjp);
 
-//                final DenseMatrix64F Pjp = new DenseMatrix64F(dimTrait, dimTrait);
-            final DenseMatrix64F Pjp = matrix1;
+//                final DMatrixRMaj Pjp = new DMatrixRMaj(dimTrait, dimTrait);
+            final DMatrixRMaj Pjp = matrix1;
             InversionResult cj = safeInvert(Vjp, Pjp, false);
 
-//                final DenseMatrix64F Pip = new DenseMatrix64F(dimTrait, dimTrait);
-            final DenseMatrix64F Pip = matrix2;
-            CommonOps.add(Pk, Pjp, Pip);
+//                final DMatrixRMaj Pip = new DMatrixRMaj(dimTrait, dimTrait);
+            final DMatrixRMaj Pip = matrix2;
+            CommonOps_DDRM.add(Pk, Pjp, Pip);
 
-//                final DenseMatrix64F Vip = new DenseMatrix64F(dimTrait, dimTrait);
-            final DenseMatrix64F Vip = matrix3;
+//                final DMatrixRMaj Vip = new DMatrixRMaj(dimTrait, dimTrait);
+            final DMatrixRMaj Vip = matrix3;
             InversionResult cip = safeInvert(Pip, Vip, false);
 
             // C. Compute prePartial mean
@@ -180,11 +180,11 @@
             }
 
             // C. Inflate variance along node branch
-            final DenseMatrix64F Vi = Vip;
-            CommonOps.add(vi, Vd, Vip, Vi);
+            final DMatrixRMaj Vi = Vip;
+            CommonOps_DDRM.add(vi, Vd, Vip, Vi);
 
-//                final DenseMatrix64F Pi = new DenseMatrix64F(dimTrait, dimTrait);
-            final DenseMatrix64F Pi = matrix4;
+//                final DMatrixRMaj Pi = new DMatrixRMaj(dimTrait, dimTrait);
+            final DMatrixRMaj Pi = matrix4;
             InversionResult ci = safeInvert(Vi, Pi, false);
 
             // X. Store precision results for node
@@ -242,7 +242,7 @@
         final double vi = branchLengths[imo];
         final double vj = branchLengths[jmo];
 
-        final DenseMatrix64F Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait);
+        final DMatrixRMaj Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait);
 
         if (DEBUG) {
             System.err.println("variance diffusion: " + Vd);
@@ -269,11 +269,11 @@
             final double lpi = partials[ibo + dimTrait + 2 * dimTrait * dimTrait];
             final double lpj = partials[jbo + dimTrait + 2 * dimTrait * dimTrait];
 
-            final DenseMatrix64F Pi = wrap(partials, ibo + dimTrait, dimTrait, dimTrait);
-            final DenseMatrix64F Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait);
+            final DMatrixRMaj Pi = wrap(partials, ibo + dimTrait, dimTrait, dimTrait);
+            final DMatrixRMaj Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait);
 
-            final DenseMatrix64F Vi = wrap(partials, ibo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
-            final DenseMatrix64F Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+            final DMatrixRMaj Vi = wrap(partials, ibo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+            final DMatrixRMaj Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
 
             if (TIMING) {
                 endTime("peel1");
@@ -286,23 +286,23 @@
             final double lpjp = Double.isInfinite(lpj) ?
                     1.0 / vj : lpj / (1.0 + lpj * vj);
 
-//                final DenseMatrix64F Vip = new DenseMatrix64F(dimTrait, dimTrait);
-//                final DenseMatrix64F Vjp = new DenseMatrix64F(dimTrait, dimTrait);
-            final DenseMatrix64F Vip = matrix0;
-            final DenseMatrix64F Vjp = matrix1;
+//                final DMatrixRMaj Vip = new DMatrixRMaj(dimTrait, dimTrait);
+//                final DMatrixRMaj Vjp = new DMatrixRMaj(dimTrait, dimTrait);
+            final DMatrixRMaj Vip = matrix0;
+            final DMatrixRMaj Vjp = matrix1;
 
-            CommonOps.add(Vi, vi, Vd, Vip);
-            CommonOps.add(Vj, vj, Vd, Vjp);
+            CommonOps_DDRM.add(Vi, vi, Vd, Vip);
+            CommonOps_DDRM.add(Vj, vj, Vd, Vjp);
 
             if (TIMING) {
                 endTime("peel2");
                 startTime("peel2a");
             }
 
-//                final DenseMatrix64F Pip = new DenseMatrix64F(dimTrait, dimTrait);
-//                final DenseMatrix64F Pjp = new DenseMatrix64F(dimTrait, dimTrait);
-            final DenseMatrix64F Pip = matrix2;
-            final DenseMatrix64F Pjp = matrix3;
+//                final DMatrixRMaj Pip = new DMatrixRMaj(dimTrait, dimTrait);
+//                final DMatrixRMaj Pjp = new DMatrixRMaj(dimTrait, dimTrait);
+            final DMatrixRMaj Pip = matrix2;
+            final DMatrixRMaj Pjp = matrix3;
 
             InversionResult ci = safeInvert(Vip, Pip, true);
             InversionResult cj = safeInvert(Vjp, Pjp, true);
@@ -317,13 +317,13 @@
             // A. Partial precision and variance (for later use) using one matrix inversion
             final double lpk = lpip + lpjp;
 
-//                final DenseMatrix64F Pk = new DenseMatrix64F(dimTrait, dimTrait);
-            final DenseMatrix64F Pk = matrix4;
+//                final DMatrixRMaj Pk = new DMatrixRMaj(dimTrait, dimTrait);
+            final DMatrixRMaj Pk = matrix4;
 
-            CommonOps.add(Pip, Pjp, Pk);
+            CommonOps_DDRM.add(Pip, Pjp, Pk);
 
-//                final DenseMatrix64F Vk = new DenseMatrix64F(dimTrait, dimTrait);
-            final DenseMatrix64F Vk = matrix5;
+//                final DMatrixRMaj Vk = new DMatrixRMaj(dimTrait, dimTrait);
+            final DMatrixRMaj Vk = matrix5;
             InversionResult ck = safeInvert(Pk, Vk, true);
 
             // B. Partial mean
@@ -433,9 +433,9 @@
                     }
                 }
 
-//                    final DenseMatrix64F Vt = new DenseMatrix64F(dimTrait, dimTrait);
-                final DenseMatrix64F Vt = matrix6;
-                CommonOps.add(Vip, Vjp, Vt);
+//                    final DMatrixRMaj Vt = new DMatrixRMaj(dimTrait, dimTrait);
+                final DMatrixRMaj Vt = matrix6;
+                CommonOps_DDRM.add(Vip, Vjp, Vt);
 
                 if (DEBUG) {
                     System.err.println("Vt: " + Vt);
@@ -445,7 +445,7 @@
                         - ck.getEffectiveDimension();
                 
                 remainder += -dimensionChange * LOG_SQRT_2_PI - 0.5 *
-//                            (Math.log(CommonOps.det(Vip)) + Math.log(CommonOps.det(Vjp)) - Math.log(CommonOps.det(Vk)))
+//                            (Math.log(CommonOps_DDRM.det(Vip)) + Math.log(CommonOps_DDRM.det(Vjp)) - Math.log(CommonOps_DDRM.det(Vk)))
                         (Math.log(ci.getDeterminant()) + Math.log(cj.getDeterminant()) + Math.log(ck.getDeterminant()))
                         - 0.5 * (SSi + SSj - SSk);
 
@@ -469,7 +469,7 @@
 
                     assert (false);
 
-                    final DenseMatrix64F Pt = new DenseMatrix64F(dimTrait, dimTrait);
+                    final DMatrixRMaj Pt = new DMatrixRMaj(dimTrait, dimTrait);
                     InversionResult ct = safeInvert(Vt, Pt, false);
 
                     int opo = dimTrait * dimTrait * trait;
@@ -551,29 +551,29 @@
         int rootOffset = dimPartial * rootBufferIndex;
         int priorOffset = dimPartial * priorBufferIndex;
 
-        final DenseMatrix64F Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait);
+        final DMatrixRMaj Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait);
 
         // TODO For each trait in parallel
         for (int trait = 0; trait < numTraits; ++trait) {
 
-            final DenseMatrix64F Proot = wrap(partials, rootOffset + dimTrait, dimTrait, dimTrait);
-            final DenseMatrix64F Pprior = wrap(partials, priorOffset + dimTrait, dimTrait, dimTrait);
+            final DMatrixRMaj Proot = wrap(partials, rootOffset + dimTrait, dimTrait, dimTrait);
+            final DMatrixRMaj Pprior = wrap(partials, priorOffset + dimTrait, dimTrait, dimTrait);
 
-            final DenseMatrix64F Vroot = wrap(partials, rootOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
-            final DenseMatrix64F Vprior = wrap(partials, priorOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+            final DMatrixRMaj Vroot = wrap(partials, rootOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+            final DMatrixRMaj Vprior = wrap(partials, priorOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
 
             // TODO Block below is for the conjugate prior ONLY
             {
-                final DenseMatrix64F Vtmp = new DenseMatrix64F(dimTrait, dimTrait);
-                CommonOps.mult(Vd, Vprior, Vtmp);
+                final DMatrixRMaj Vtmp = new DMatrixRMaj(dimTrait, dimTrait);
+                CommonOps_DDRM.mult(Vd, Vprior, Vtmp);
                 Vprior.set(Vtmp);
             }
 
-            final DenseMatrix64F Vtotal = new DenseMatrix64F(dimTrait, dimTrait);
-            CommonOps.add(Vroot, Vprior, Vtotal);
+            final DMatrixRMaj Vtotal = new DMatrixRMaj(dimTrait, dimTrait);
+            CommonOps_DDRM.add(Vroot, Vprior, Vtotal);
 
-            final DenseMatrix64F Ptotal = new DenseMatrix64F(dimTrait, dimTrait);
-            CommonOps.invert(Vtotal, Ptotal);  // TODO Can return determinant at same time to avoid extra QR decomp
+            final DMatrixRMaj Ptotal = new DMatrixRMaj(dimTrait, dimTrait);
+            CommonOps_DDRM.invert(Vtotal, Ptotal);  // TODO Can return determinant at same time to avoid extra QR decomp
 
             double SS = 0;
             for (int g = 0; g < dimTrait; ++g) {
@@ -586,7 +586,7 @@
                 }
             }
 
-            final double logLike = -dimTrait * LOG_SQRT_2_PI - 0.5 * Math.log(CommonOps.det(Vtotal)) - 0.5 * SS;
+            final double logLike = -dimTrait * LOG_SQRT_2_PI - 0.5 * Math.log(CommonOps_DDRM.det(Vtotal)) - 0.5 * SS;
 
             final double remainder = remainders[rootBufferIndex * numTraits + trait];
             logLikelihoods[trait] = logLike + remainder;
--- a/src/dr/evomodel/treedatalikelihood/continuous/cdi/SafeMultivariateIntegrator.java
+++ b/src/dr/evomodel/treedatalikelihood/continuous/cdi/SafeMultivariateIntegrator.java
@@ -2,8 +2,8 @@
 
 import dr.math.matrixAlgebra.WrappedVector;
 import dr.math.matrixAlgebra.missingData.InversionResult;
-import org.ejml.data.DenseMatrix64F;
-import org.ejml.ops.CommonOps;
+import org.ejml.data.DMatrixRMaj;
+import org.ejml.dense.row.CommonOps_DDRM;
 
 import static dr.math.matrixAlgebra.missingData.InversionResult.Code.NOT_OBSERVED;
 import static dr.math.matrixAlgebra.missingData.MissingOps.*;
@@ -53,13 +53,13 @@
 
 //    private final Map<String, Long> times;
 //
-//    private DenseMatrix64F matrix0;
-//    private DenseMatrix64F matrix1;
-//    private DenseMatrix64F matrix2;
-//    private DenseMatrix64F matrix3;
-//    private DenseMatrix64F matrix4;
-//    private DenseMatrix64F matrix5;
-//    private DenseMatrix64F matrix6;
+//    private DMatrixRMaj matrix0;
+//    private DMatrixRMaj matrix1;
+//    private DMatrixRMaj matrix2;
+//    private DMatrixRMaj matrix3;
+//    private DMatrixRMaj matrix4;
+//    private DMatrixRMaj matrix5;
+//    private DMatrixRMaj matrix6;
 //
 //    private double[] vector0;
 
@@ -67,13 +67,13 @@
 //        inverseDiffusions = new double[dimTrait * dimTrait * diffusionCount];
 //
 //        vector0 = new double[dimTrait];
-//        matrix0 = new DenseMatrix64F(dimTrait, dimTrait);
-//        matrix1 = new DenseMatrix64F(dimTrait, dimTrait);
-//        matrix2 = new DenseMatrix64F(dimTrait, dimTrait);
-//        matrix3 = new DenseMatrix64F(dimTrait, dimTrait);
-//        matrix4 = new DenseMatrix64F(dimTrait, dimTrait);
-//        matrix5 = new DenseMatrix64F(dimTrait, dimTrait);
-//        matrix6 = new DenseMatrix64F(dimTrait, dimTrait);
+//        matrix0 = new DMatrixRMaj(dimTrait, dimTrait);
+//        matrix1 = new DMatrixRMaj(dimTrait, dimTrait);
+//        matrix2 = new DMatrixRMaj(dimTrait, dimTrait);
+//        matrix3 = new DMatrixRMaj(dimTrait, dimTrait);
+//        matrix4 = new DMatrixRMaj(dimTrait, dimTrait);
+//        matrix5 = new DMatrixRMaj(dimTrait, dimTrait);
+//        matrix6 = new DMatrixRMaj(dimTrait, dimTrait);
 //    }
 
 //    @Override
@@ -83,9 +83,9 @@
 //        assert (inverseDiffusions != null);
 //
 //        final int offset = dimTrait * dimTrait * precisionIndex;
-//        DenseMatrix64F precision = wrap(diffusions, offset, dimTrait, dimTrait);
-//        DenseMatrix64F variance = new DenseMatrix64F(dimTrait, dimTrait);
-//        CommonOps.invert(precision, variance);
+//        DMatrixRMaj precision = wrap(diffusions, offset, dimTrait, dimTrait);
+//        DMatrixRMaj variance = new DMatrixRMaj(dimTrait, dimTrait);
+//        CommonOps_DDRM.invert(precision, variance);
 //        unwrap(variance, inverseDiffusions, offset);
 //
 //        if (DEBUG) {
@@ -121,7 +121,7 @@
 //        final double vi = branchLengths[imo];
 //        final double vj = branchLengths[jmo];
 //
-//        final DenseMatrix64F Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait);
+//        final DMatrixRMaj Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait);
 //
 //        if (DEBUG) {
 //            System.err.println("updatePreOrderPartial for node " + iBuffer);
@@ -134,27 +134,27 @@
 //        for (int trait = 0; trait < numTraits; ++trait) {
 //
 //            // A. Get current precision of k and j
-//            final DenseMatrix64F Pk = wrap(prePartials, kbo + dimTrait, dimTrait, dimTrait);
-////                final DenseMatrix64F Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait);
+//            final DMatrixRMaj Pk = wrap(prePartials, kbo + dimTrait, dimTrait, dimTrait);
+////                final DMatrixRMaj Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait);
 //
-////                final DenseMatrix64F Vk = wrap(prePartials, kbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
-//            final DenseMatrix64F Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+////                final DMatrixRMaj Vk = wrap(prePartials, kbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+//            final DMatrixRMaj Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
 //
 //            // B. Inflate variance along sibling branch using matrix inversion
-////                final DenseMatrix64F Vjp = new DenseMatrix64F(dimTrait, dimTrait);
-//            final DenseMatrix64F Vjp = matrix0;
-//            CommonOps.add(Vj, vj, Vd, Vjp);
+////                final DMatrixRMaj Vjp = new DMatrixRMaj(dimTrait, dimTrait);
+//            final DMatrixRMaj Vjp = matrix0;
+//            CommonOps_DDRM.add(Vj, vj, Vd, Vjp);
 //
-////                final DenseMatrix64F Pjp = new DenseMatrix64F(dimTrait, dimTrait);
-//            final DenseMatrix64F Pjp = matrix1;
+////                final DMatrixRMaj Pjp = new DMatrixRMaj(dimTrait, dimTrait);
+//            final DMatrixRMaj Pjp = matrix1;
 //            InversionResult cj = safeInvert(Vjp, Pjp, false);
 //
-////                final DenseMatrix64F Pip = new DenseMatrix64F(dimTrait, dimTrait);
-//            final DenseMatrix64F Pip = matrix2;
-//            CommonOps.add(Pk, Pjp, Pip);
+////                final DMatrixRMaj Pip = new DMatrixRMaj(dimTrait, dimTrait);
+//            final DMatrixRMaj Pip = matrix2;
+//            CommonOps_DDRM.add(Pk, Pjp, Pip);
 //
-////                final DenseMatrix64F Vip = new DenseMatrix64F(dimTrait, dimTrait);
-//            final DenseMatrix64F Vip = matrix3;
+////                final DMatrixRMaj Vip = new DMatrixRMaj(dimTrait, dimTrait);
+//            final DMatrixRMaj Vip = matrix3;
 //            InversionResult cip = safeInvert(Pip, Vip, false);
 //
 //            // C. Compute prePartial mean
@@ -177,11 +177,11 @@
 //            }
 //
 //            // C. Inflate variance along node branch
-//            final DenseMatrix64F Vi = Vip;
-//            CommonOps.add(vi, Vd, Vip, Vi);
+//            final DMatrixRMaj Vi = Vip;
+//            CommonOps_DDRM.add(vi, Vd, Vip, Vi);
 //
-////                final DenseMatrix64F Pi = new DenseMatrix64F(dimTrait, dimTrait);
-//            final DenseMatrix64F Pi = matrix4;
+////                final DMatrixRMaj Pi = new DMatrixRMaj(dimTrait, dimTrait);
+//            final DMatrixRMaj Pi = matrix4;
 //            InversionResult ci = safeInvert(Vi, Pi, false);
 //
 //            // X. Store precision results for node
@@ -239,8 +239,8 @@
         final double vi = branchLengths[imo];
         final double vj = branchLengths[jmo];
 
-        final DenseMatrix64F Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait);
-        final DenseMatrix64F Pd = wrap(diffusions, precisionOffset, dimTrait, dimTrait);
+        final DMatrixRMaj Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait);
+        final DMatrixRMaj Pd = wrap(diffusions, precisionOffset, dimTrait, dimTrait);
 
         if (DEBUG) {
             System.err.println("variance diffusion: " + Vd);
@@ -267,8 +267,8 @@
             final double lpi = partials[ibo + dimTrait + 2 * dimTrait * dimTrait];
             final double lpj = partials[jbo + dimTrait + 2 * dimTrait * dimTrait];
 
-            final DenseMatrix64F Pi = wrap(partials, ibo + dimTrait, dimTrait, dimTrait);
-            final DenseMatrix64F Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait);
+            final DMatrixRMaj Pi = wrap(partials, ibo + dimTrait, dimTrait, dimTrait);
+            final DMatrixRMaj Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait);
 
             if (TIMING) {
                 endTime("peel1");
@@ -284,8 +284,8 @@
             InversionResult ci;
             InversionResult cj;
 
-            final DenseMatrix64F Pip = matrix2;
-            final DenseMatrix64F Pjp = matrix3;
+            final DMatrixRMaj Pip = matrix2;
+            final DMatrixRMaj Pjp = matrix3;
 
 //            boolean useVariance = anyDiagonalInfinities(Pi) || anyDiagonalInfinities(Pj);
             final boolean useVariancei = anyDiagonalInfinities(Pi);
@@ -293,77 +293,77 @@
 
             if (useVariancei) {
 
-                final DenseMatrix64F Vip = matrix0;
-                final DenseMatrix64F Vi = wrap(partials, ibo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
-                CommonOps.add(Vi, vi, Vd, Vip);
+                final DMatrixRMaj Vip = matrix0;
+                final DMatrixRMaj Vi = wrap(partials, ibo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+                CommonOps_DDRM.add(Vi, vi, Vd, Vip);
                 ci = safeInvert(Vip, Pip, true);
 
             } else {
 
-                final DenseMatrix64F PiPlusPd = matrix0;
-                CommonOps.add(Pi, 1.0 / vi, Pd, PiPlusPd);
-                final DenseMatrix64F PiPlusPdInv = new DenseMatrix64F(dimTrait, dimTrait);
+                final DMatrixRMaj PiPlusPd = matrix0;
+                CommonOps_DDRM.add(Pi, 1.0 / vi, Pd, PiPlusPd);
+                final DMatrixRMaj PiPlusPdInv = new DMatrixRMaj(dimTrait, dimTrait);
                 safeInvert(PiPlusPd, PiPlusPdInv, false);
-                CommonOps.mult(PiPlusPdInv, Pi, Pip);
-                CommonOps.mult(Pi, Pip, PiPlusPdInv);
-                CommonOps.add(Pi, -1, PiPlusPdInv, Pip);
+                CommonOps_DDRM.mult(PiPlusPdInv, Pi, Pip);
+                CommonOps_DDRM.mult(Pi, Pip, PiPlusPdInv);
+                CommonOps_DDRM.add(Pi, -1, PiPlusPdInv, Pip);
                 ci = safeDeterminant(Pip, false);
             }
 
             if (useVariancej) {
 
-                final DenseMatrix64F Vjp = matrix1;
-                final DenseMatrix64F Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
-                CommonOps.add(Vj, vj, Vd, Vjp);
+                final DMatrixRMaj Vjp = matrix1;
+                final DMatrixRMaj Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+                CommonOps_DDRM.add(Vj, vj, Vd, Vjp);
                 cj = safeInvert(Vjp, Pjp, true);
 
             } else {
 
-                final DenseMatrix64F PjPlusPd = matrix1;
-                CommonOps.add(Pj, 1.0 / vj, Pd, PjPlusPd);
-                final DenseMatrix64F PjPlusPdInv = new DenseMatrix64F(dimTrait, dimTrait);
+                final DMatrixRMaj PjPlusPd = matrix1;
+                CommonOps_DDRM.add(Pj, 1.0 / vj, Pd, PjPlusPd);
+                final DMatrixRMaj PjPlusPdInv = new DMatrixRMaj(dimTrait, dimTrait);
                 safeInvert(PjPlusPd, PjPlusPdInv, false);
-                CommonOps.mult(PjPlusPdInv, Pj, Pjp);
-                CommonOps.mult(Pj, Pjp, PjPlusPdInv);
-                CommonOps.add(Pj, -1, PjPlusPdInv, Pjp);
+                CommonOps_DDRM.mult(PjPlusPdInv, Pj, Pjp);
+                CommonOps_DDRM.mult(Pj, Pjp, PjPlusPdInv);
+                CommonOps_DDRM.add(Pj, -1, PjPlusPdInv, Pjp);
                 cj = safeDeterminant(Pjp, false);
             }
 
 //            if (useVariance) {
 //
-//                final DenseMatrix64F Vip = matrix0;
-//                final DenseMatrix64F Vjp = matrix1;
+//                final DMatrixRMaj Vip = matrix0;
+//                final DMatrixRMaj Vjp = matrix1;
 //
-//                final DenseMatrix64F Vi = wrap(partials, ibo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
-//                final DenseMatrix64F Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+//                final DMatrixRMaj Vi = wrap(partials, ibo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+//                final DMatrixRMaj Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
 //
-//                CommonOps.add(Vi, vi, Vd, Vip);
-//                CommonOps.add(Vj, vj, Vd, Vjp);
+//                CommonOps_DDRM.add(Vi, vi, Vd, Vip);
+//                CommonOps_DDRM.add(Vj, vj, Vd, Vjp);
 //
 //                ci = safeInvert(Vip, Pip, true);
 //                cj = safeInvert(Vjp, Pjp, true);
 //            } else {
 //
-//                final DenseMatrix64F PiPlusPd = matrix0;
-//                final DenseMatrix64F PjPlusPd = matrix1;
+//                final DMatrixRMaj PiPlusPd = matrix0;
+//                final DMatrixRMaj PjPlusPd = matrix1;
 //
-//                CommonOps.add(Pi, 1.0 / vi, Pd, PiPlusPd);
-//                CommonOps.add(Pj, 1.0 / vj, Pd, PjPlusPd);
+//                CommonOps_DDRM.add(Pi, 1.0 / vi, Pd, PiPlusPd);
+//                CommonOps_DDRM.add(Pj, 1.0 / vj, Pd, PjPlusPd);
 //
-//                final DenseMatrix64F PiPlusPdInv = new DenseMatrix64F(dimTrait, dimTrait);
-//                final DenseMatrix64F PjPlusPdInv = new DenseMatrix64F(dimTrait, dimTrait);
+//                final DMatrixRMaj PiPlusPdInv = new DMatrixRMaj(dimTrait, dimTrait);
+//                final DMatrixRMaj PjPlusPdInv = new DMatrixRMaj(dimTrait, dimTrait);
 //
 //                safeInvert(PiPlusPd, PiPlusPdInv, false);
 //                safeInvert(PjPlusPd, PjPlusPdInv, false);
 //
-//                CommonOps.mult(PiPlusPdInv, Pi, Pip);
-//                CommonOps.mult(PjPlusPdInv, Pj, Pjp);
+//                CommonOps_DDRM.mult(PiPlusPdInv, Pi, Pip);
+//                CommonOps_DDRM.mult(PjPlusPdInv, Pj, Pjp);
 //
-//                CommonOps.mult(Pi, Pip, PiPlusPdInv);
-//                CommonOps.mult(Pj, Pjp, PjPlusPdInv);
+//                CommonOps_DDRM.mult(Pi, Pip, PiPlusPdInv);
+//                CommonOps_DDRM.mult(Pj, Pjp, PjPlusPdInv);
 //
-//                CommonOps.add(Pi, -1, PiPlusPdInv, Pip);
-//                CommonOps.add(Pj, -1, PjPlusPdInv, Pjp);
+//                CommonOps_DDRM.add(Pi, -1, PiPlusPdInv, Pip);
+//                CommonOps_DDRM.add(Pj, -1, PjPlusPdInv, Pjp);
 //
 //                ci = safeDeterminant(Pip, false);
 //                cj = safeDeterminant(Pjp, false);
@@ -385,12 +385,12 @@
             // A. Partial precision and variance (for later use) using one matrix inversion
             final double lpk = lpip + lpjp;
 
-//                final DenseMatrix64F Pk = new DenseMatrix64F(dimTrait, dimTrait);
-            final DenseMatrix64F Pk = matrix4;
-            CommonOps.add(Pip, Pjp, Pk);
+//                final DMatrixRMaj Pk = new DMatrixRMaj(dimTrait, dimTrait);
+            final DMatrixRMaj Pk = matrix4;
+            CommonOps_DDRM.add(Pip, Pjp, Pk);
 
-//                final DenseMatrix64F Vk = new DenseMatrix64F(dimTrait, dimTrait);
-//            final DenseMatrix64F Vk = matrix5;
+//                final DMatrixRMaj Vk = new DMatrixRMaj(dimTrait, dimTrait);
+//            final DMatrixRMaj Vk = matrix5;
 
 //            if (useVariance) {
 ////            InversionResult ck =
@@ -524,9 +524,9 @@
                     }
                 }
 
-//                    final DenseMatrix64F Vt = new DenseMatrix64F(dimTrait, dimTrait);
-//                final DenseMatrix64F Vt = matrix6;
-//                CommonOps.add(Vip, Vjp, Vt);
+//                    final DMatrixRMaj Vt = new DMatrixRMaj(dimTrait, dimTrait);
+//                final DMatrixRMaj Vt = matrix6;
+//                CommonOps_DDRM.add(Vip, Vjp, Vt);
 
 //                if (DEBUG) {
 //                    System.err.println("Vt: " + Vt);
@@ -536,16 +536,16 @@
                         - ck.getEffectiveDimension();
 
 //                    System.err.println(ci.getDeterminant());
-//                    System.err.println(CommonOps.det(Vip));
+//                    System.err.println(CommonOps_DDRM.det(Vip));
 //
 //                    System.err.println(cj.getDeterminant());
-//                    System.err.println(CommonOps.det(Vjp));
+//                    System.err.println(CommonOps_DDRM.det(Vjp));
 //
 //                    System.err.println(1.0 / ck.getDeterminant());
-//                    System.err.println(CommonOps.det(Vk));
+//                    System.err.println(CommonOps_DDRM.det(Vk));
 
                 remainder += -dimensionChange * LOG_SQRT_2_PI - 0.5 *
-//                            (Math.log(CommonOps.det(Vip)) + Math.log(CommonOps.det(Vjp)) - Math.log(CommonOps.det(Vk)))
+//                            (Math.log(CommonOps_DDRM.det(Vip)) + Math.log(CommonOps_DDRM.det(Vjp)) - Math.log(CommonOps_DDRM.det(Vk)))
                         (Math.log(ci.getDeterminant()) + Math.log(cj.getDeterminant()) + Math.log(ck.getDeterminant()))
                         - 0.5 * (SSi + SSj - SSk);
 
@@ -632,42 +632,42 @@
         int rootOffset = dimPartial * rootBufferIndex;
         int priorOffset = dimPartial * priorBufferIndex;
 
-        final DenseMatrix64F Pd = wrap(diffusions, precisionOffset, dimTrait, dimTrait);
-//        final DenseMatrix64F Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait);
+        final DMatrixRMaj Pd = wrap(diffusions, precisionOffset, dimTrait, dimTrait);
+//        final DMatrixRMaj Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait);
 
         // TODO For each trait in parallel
         for (int trait = 0; trait < numTraits; ++trait) {
 
-            final DenseMatrix64F Proot = wrap(partials, rootOffset + dimTrait, dimTrait, dimTrait);
-            final DenseMatrix64F Pprior = wrap(partials, priorOffset + dimTrait, dimTrait, dimTrait);
+            final DMatrixRMaj Proot = wrap(partials, rootOffset + dimTrait, dimTrait, dimTrait);
+            final DMatrixRMaj Pprior = wrap(partials, priorOffset + dimTrait, dimTrait, dimTrait);
 
-//            final DenseMatrix64F Vroot = wrap(partials, rootOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
-//            final DenseMatrix64F Vprior = wrap(partials, priorOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+//            final DMatrixRMaj Vroot = wrap(partials, rootOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+//            final DMatrixRMaj Vprior = wrap(partials, priorOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
 
             // TODO Block below is for the conjugate prior ONLY
             {
-//                final DenseMatrix64F Vtmp = new DenseMatrix64F(dimTrait, dimTrait);
-//                CommonOps.mult(Vd, Vprior, Vtmp);
+//                final DMatrixRMaj Vtmp = new DMatrixRMaj(dimTrait, dimTrait);
+//                CommonOps_DDRM.mult(Vd, Vprior, Vtmp);
 //                Vprior.set(Vtmp);
 
-                final DenseMatrix64F Ptmp = new DenseMatrix64F(dimTrait, dimTrait);
-                CommonOps.mult(Pd, Pprior, Ptmp);
+                final DMatrixRMaj Ptmp = new DMatrixRMaj(dimTrait, dimTrait);
+                CommonOps_DDRM.mult(Pd, Pprior, Ptmp);
                 Pprior.set(Ptmp); // TODO What does this do?
             }
 
-            final DenseMatrix64F Vtotal = new DenseMatrix64F(dimTrait, dimTrait);
-//            CommonOps.add(Vroot, Vprior, Vtotal);
+            final DMatrixRMaj Vtotal = new DMatrixRMaj(dimTrait, dimTrait);
+//            CommonOps_DDRM.add(Vroot, Vprior, Vtotal);
 
-            final DenseMatrix64F Ptotal = new DenseMatrix64F(dimTrait, dimTrait);
-            CommonOps.invert(Vtotal, Ptotal);  // TODO Can return determinant at same time to avoid extra QR decomp
+            final DMatrixRMaj Ptotal = new DMatrixRMaj(dimTrait, dimTrait);
+            CommonOps_DDRM.invert(Vtotal, Ptotal);  // TODO Can return determinant at same time to avoid extra QR decomp
 
-            final DenseMatrix64F tmp1 = new DenseMatrix64F(dimTrait, dimTrait);
-            final DenseMatrix64F tmp2 = new DenseMatrix64F(dimTrait, dimTrait);
-            CommonOps.add(Proot, Pprior, Ptotal);
-            CommonOps.invert(Ptotal, Vtotal);
-            CommonOps.mult(Vtotal, Proot, tmp1);
-            CommonOps.mult(Proot, tmp1, tmp2);
-            CommonOps.add(Proot, -1.0, tmp2, Ptotal);
+            final DMatrixRMaj tmp1 = new DMatrixRMaj(dimTrait, dimTrait);
+            final DMatrixRMaj tmp2 = new DMatrixRMaj(dimTrait, dimTrait);
+            CommonOps_DDRM.add(Proot, Pprior, Ptotal);
+            CommonOps_DDRM.invert(Ptotal, Vtotal);
+            CommonOps_DDRM.mult(Vtotal, Proot, tmp1);
+            CommonOps_DDRM.mult(Proot, tmp1, tmp2);
+            CommonOps_DDRM.add(Proot, -1.0, tmp2, Ptotal);
 
             double SS = 0;
             for (int g = 0; g < dimTrait; ++g) {
@@ -681,8 +681,8 @@
             }
 
             final double logLike = -dimTrait * LOG_SQRT_2_PI
-//                    - 0.5 * Math.log(CommonOps.det(Vtotal))
-                    + 0.5 * Math.log(CommonOps.det(Ptotal))
+//                    - 0.5 * Math.log(CommonOps_DDRM.det(Vtotal))
+                    + 0.5 * Math.log(CommonOps_DDRM.det(Ptotal))
                     - 0.5 * SS;
 
             final double remainder = remainders[rootBufferIndex * numTraits + trait];
--- a/src/dr/evomodel/treedatalikelihood/continuous/cdi/SafeMultivariateWithDriftIntegrator.java
+++ b/src/dr/evomodel/treedatalikelihood/continuous/cdi/SafeMultivariateWithDriftIntegrator.java
@@ -2,8 +2,8 @@
 
 import dr.math.matrixAlgebra.WrappedVector;
 import dr.math.matrixAlgebra.missingData.InversionResult;
-import org.ejml.data.DenseMatrix64F;
-import org.ejml.ops.CommonOps;
+import org.ejml.data.DMatrixRMaj;
+import org.ejml.dense.row.CommonOps_DDRM;
 
 import static dr.math.matrixAlgebra.missingData.InversionResult.Code.NOT_OBSERVED;
 import static dr.math.matrixAlgebra.missingData.MissingOps.*;
@@ -65,9 +65,9 @@
         assert (inverseDiffusions != null);
 
         final int offset = dimTrait * dimTrait * precisionIndex;
-        DenseMatrix64F precision = wrap(diffusions, offset, dimTrait, dimTrait);
-        DenseMatrix64F variance = new DenseMatrix64F(dimTrait, dimTrait);
-        CommonOps.invert(precision, variance);
+        DMatrixRMaj precision = wrap(diffusions, offset, dimTrait, dimTrait);
+        DMatrixRMaj variance = new DMatrixRMaj(dimTrait, dimTrait);
+        CommonOps_DDRM.invert(precision, variance);
         unwrap(variance, inverseDiffusions, offset);
 
         if (DEBUG) {
@@ -199,14 +199,14 @@
 //        final double vi = variances[imo];
 //        final double vj = variances[jmo];
 
-        final DenseMatrix64F Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait);
-//        final DenseMatrix64F Pd = wrap(diffusions, precisionOffset, dimTrait, dimTrait);
+        final DMatrixRMaj Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait);
+//        final DMatrixRMaj Pd = wrap(diffusions, precisionOffset, dimTrait, dimTrait);
 
-        final DenseMatrix64F Vdi = wrap(variances, imo, dimTrait, dimTrait);
-        final DenseMatrix64F Vdj = wrap(variances, jmo, dimTrait, dimTrait);
+        final DMatrixRMaj Vdi = wrap(variances, imo, dimTrait, dimTrait);
+        final DMatrixRMaj Vdj = wrap(variances, jmo, dimTrait, dimTrait);
 
-        final DenseMatrix64F Pdi = wrap(precisions, imo, dimTrait, dimTrait); // TODO Only if needed
-        final DenseMatrix64F Pdj = wrap(precisions, jmo, dimTrait, dimTrait); // TODO Only if needed
+        final DMatrixRMaj Pdi = wrap(precisions, imo, dimTrait, dimTrait); // TODO Only if needed
+        final DMatrixRMaj Pdj = wrap(precisions, jmo, dimTrait, dimTrait); // TODO Only if needed
 
         // TODO End fix
 
@@ -237,8 +237,8 @@
 //            final double lpi = partials[ibo + dimTrait + 2 * dimTrait * dimTrait];
 //            final double lpj = partials[jbo + dimTrait + 2 * dimTrait * dimTrait];
 
-            final DenseMatrix64F Pi = wrap(partials, ibo + dimTrait, dimTrait, dimTrait);
-            final DenseMatrix64F Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait);
+            final DMatrixRMaj Pi = wrap(partials, ibo + dimTrait, dimTrait, dimTrait);
+            final DMatrixRMaj Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait);
 
             if (TIMING) {
                 endTime("peel1");
@@ -254,8 +254,8 @@
             InversionResult ci;
             InversionResult cj;
 
-            final DenseMatrix64F Pip = matrix2;
-            final DenseMatrix64F Pjp = matrix3;
+            final DMatrixRMaj Pip = matrix2;
+            final DMatrixRMaj Pjp = matrix3;
 
 //            boolean useVariance = anyDiagonalInfinities(Pi) || anyDiagonalInfinities(Pj);
             final boolean useVariancei = anyDiagonalInfinities(Pi);
@@ -263,43 +263,43 @@
 
             if (useVariancei) {
 
-                final DenseMatrix64F Vip = matrix0;
-                final DenseMatrix64F Vi = wrap(partials, ibo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
-//                CommonOps.add(Vi, vi, Vd, Vip);  // TODO Fix
-                CommonOps.add(Vi, Vdi, Vip);
+                final DMatrixRMaj Vip = matrix0;
+                final DMatrixRMaj Vi = wrap(partials, ibo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+//                CommonOps_DDRM.add(Vi, vi, Vd, Vip);  // TODO Fix
+                CommonOps_DDRM.add(Vi, Vdi, Vip);
                 ci = safeInvert(Vip, Pip, true);
 
             } else {
 
-                final DenseMatrix64F PiPlusPd = matrix0;
-//                CommonOps.add(Pi, 1.0 / vi, Pd, PiPlusPd); // TODO Fix
-                CommonOps.add(Pi, Pdi, PiPlusPd);
-                final DenseMatrix64F PiPlusPdInv = new DenseMatrix64F(dimTrait, dimTrait);
+                final DMatrixRMaj PiPlusPd = matrix0;
+//                CommonOps_DDRM.add(Pi, 1.0 / vi, Pd, PiPlusPd); // TODO Fix
+                CommonOps_DDRM.add(Pi, Pdi, PiPlusPd);
+                final DMatrixRMaj PiPlusPdInv = new DMatrixRMaj(dimTrait, dimTrait);
                 safeInvert(PiPlusPd, PiPlusPdInv, false);
-                CommonOps.mult(PiPlusPdInv, Pi, Pip);
-                CommonOps.mult(Pi, Pip, PiPlusPdInv);
-                CommonOps.add(Pi, -1, PiPlusPdInv, Pip);
+                CommonOps_DDRM.mult(PiPlusPdInv, Pi, Pip);
+                CommonOps_DDRM.mult(Pi, Pip, PiPlusPdInv);
+                CommonOps_DDRM.add(Pi, -1, PiPlusPdInv, Pip);
                 ci = safeDeterminant(Pip, false);
             }
 
             if (useVariancej) {
 
-                final DenseMatrix64F Vjp = matrix1;
-                final DenseMatrix64F Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
-//                CommonOps.add(Vj, vj, Vd, Vjp); // TODO Fix
-                CommonOps.add(Vj, Vdj, Vjp);
+                final DMatrixRMaj Vjp = matrix1;
+                final DMatrixRMaj Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+//                CommonOps_DDRM.add(Vj, vj, Vd, Vjp); // TODO Fix
+                CommonOps_DDRM.add(Vj, Vdj, Vjp);
                 cj = safeInvert(Vjp, Pjp, true);
 
             } else {
 
-                final DenseMatrix64F PjPlusPd = matrix1;
-//                CommonOps.add(Pj, 1.0 / vj, Pd, PjPlusPd); // TODO Fix
-                CommonOps.add(Pj, Pdj, PjPlusPd);
-                final DenseMatrix64F PjPlusPdInv = new DenseMatrix64F(dimTrait, dimTrait);
+                final DMatrixRMaj PjPlusPd = matrix1;
+//                CommonOps_DDRM.add(Pj, 1.0 / vj, Pd, PjPlusPd); // TODO Fix
+                CommonOps_DDRM.add(Pj, Pdj, PjPlusPd);
+                final DMatrixRMaj PjPlusPdInv = new DMatrixRMaj(dimTrait, dimTrait);
                 safeInvert(PjPlusPd, PjPlusPdInv, false);
-                CommonOps.mult(PjPlusPdInv, Pj, Pjp);
-                CommonOps.mult(Pj, Pjp, PjPlusPdInv);
-                CommonOps.add(Pj, -1, PjPlusPdInv, Pjp);
+                CommonOps_DDRM.mult(PjPlusPdInv, Pj, Pjp);
+                CommonOps_DDRM.mult(Pj, Pjp, PjPlusPdInv);
+                CommonOps_DDRM.add(Pj, -1, PjPlusPdInv, Pjp);
                 cj = safeDeterminant(Pjp, false);
             }
 
@@ -313,12 +313,12 @@
             // A. Partial precision and variance (for later use) using one matrix inversion
 //            final double lpk = lpip + lpjp;
 
-//                final DenseMatrix64F Pk = new DenseMatrix64F(dimTrait, dimTrait);
-            final DenseMatrix64F Pk = matrix4;
-            CommonOps.add(Pip, Pjp, Pk);
+//                final DMatrixRMaj Pk = new DMatrixRMaj(dimTrait, dimTrait);
+            final DMatrixRMaj Pk = matrix4;
+            CommonOps_DDRM.add(Pip, Pjp, Pk);
 
-//                final DenseMatrix64F Vk = new DenseMatrix64F(dimTrait, dimTrait);
-//            final DenseMatrix64F Vk = matrix5;
+//                final DMatrixRMaj Vk = new DMatrixRMaj(dimTrait, dimTrait);
+//            final DMatrixRMaj Vk = matrix5;
 
 //            if (useVariance) {
 ////            InversionResult ck =
@@ -475,9 +475,9 @@
                     }
                 }
 
-//                    final DenseMatrix64F Vt = new DenseMatrix64F(dimTrait, dimTrait);
-//                final DenseMatrix64F Vt = matrix6;
-//                CommonOps.add(Vip, Vjp, Vt);
+//                    final DMatrixRMaj Vt = new DMatrixRMaj(dimTrait, dimTrait);
+//                final DMatrixRMaj Vt = matrix6;
+//                CommonOps_DDRM.add(Vip, Vjp, Vt);
 
 //                if (DEBUG) {
 //                    System.err.println("Vt: " + Vt);
@@ -487,16 +487,16 @@
                         - ck.getEffectiveDimension();
 
 //                    System.err.println(ci.getDeterminant());
-//                    System.err.println(CommonOps.det(Vip));
+//                    System.err.println(CommonOps_DDRM.det(Vip));
 //
 //                    System.err.println(cj.getDeterminant());
-//                    System.err.println(CommonOps.det(Vjp));
+//                    System.err.println(CommonOps_DDRM.det(Vjp));
 //
 //                    System.err.println(1.0 / ck.getDeterminant());
-//                    System.err.println(CommonOps.det(Vk));
+//                    System.err.println(CommonOps_DDRM.det(Vk));
 
                 remainder += -dimensionChange * LOG_SQRT_2_PI - 0.5 *
-//                            (Math.log(CommonOps.det(Vip)) + Math.log(CommonOps.det(Vjp)) - Math.log(CommonOps.det(Vk)))
+//                            (Math.log(CommonOps_DDRM.det(Vip)) + Math.log(CommonOps_DDRM.det(Vjp)) - Math.log(CommonOps_DDRM.det(Vk)))
                         (Math.log(ci.getDeterminant()) + Math.log(cj.getDeterminant()) + Math.log(ck.getDeterminant()))
                         - 0.5 * (SSi + SSj - SSk);
 
--- a/src/dr/evomodel/treedatalikelihood/continuous/IntegratedFactorAnalysisLikelihood.java
+++ b/src/dr/evomodel/treedatalikelihood/continuous/IntegratedFactorAnalysisLikelihood.java
@@ -38,7 +38,7 @@
 import dr.math.matrixAlgebra.WrappedVector;
 import dr.math.matrixAlgebra.missingData.InversionResult;
 import dr.xml.*;
-import org.ejml.data.DenseMatrix64F;
+import org.ejml.data.DMatrixRMaj;
 
 import java.util.ArrayList;
 import java.util.Arrays;
@@ -229,7 +229,7 @@
 
     private final double nuggetPrecision;
 
-    private void computePrecisionForTaxon(final DenseMatrix64F precision, final int taxon,
+    private void computePrecisionForTaxon(final DMatrixRMaj precision, final int taxon,
                                            final int numFactors) {
 
         final double[] observed = observedIndicators[taxon];
@@ -249,7 +249,7 @@
         }
     }
 
-    private InversionResult fillInMeanForTaxon(final WrappedVector output, final DenseMatrix64F precision, final int taxon) {
+    private InversionResult fillInMeanForTaxon(final WrappedVector output, final DMatrixRMaj precision, final int taxon) {
 
         final double[] observed = observedIndicators[taxon];
         final Parameter Y = traitParameter.getParameter(taxon);
@@ -269,8 +269,8 @@
             tmp[row] = sum;
         }
 
-        DenseMatrix64F B = DenseMatrix64F.wrap(numFactors, 1, tmp);
-        DenseMatrix64F X = DenseMatrix64F.wrap(numFactors, 1, tmp2);
+        DMatrixRMaj B = DMatrixRMaj.wrap(numFactors, 1, tmp);
+        DMatrixRMaj X = DMatrixRMaj.wrap(numFactors, 1, tmp2);
 
         InversionResult ci = safeSolve(precision, B, X, true);
 
@@ -294,7 +294,7 @@
         return sum;
     }
 
-    private double computeFactorInnerProduct(final WrappedVector mean, final DenseMatrix64F precision) {
+    private double computeFactorInnerProduct(final WrappedVector mean, final DMatrixRMaj precision) {
         // Compute \mu_i^t P_i \mu^t
         double sum = 0;
         for (int row = 0; row < numFactors; ++row) {
@@ -334,7 +334,7 @@
         return logDet;
     }
 
-    private void makeCompletedUnobserved(final DenseMatrix64F matrix, double diagonal) {
+    private void makeCompletedUnobserved(final DMatrixRMaj matrix, double diagonal) {
         for (int row = 0; row < numFactors; ++row) {
             for (int col = 0; col < numFactors; ++col) {
                 double x = (row == col) ? diagonal : 0.0;
@@ -345,8 +345,8 @@
 
     private void computePartialsAndRemainders() {
         
-        final DenseMatrix64F precision = new DenseMatrix64F(numFactors, numFactors);
-        final DenseMatrix64F variance = new DenseMatrix64F(numFactors, numFactors);
+        final DMatrixRMaj precision = new DMatrixRMaj(numFactors, numFactors);
+        final DMatrixRMaj variance = new DMatrixRMaj(numFactors, numFactors);
 
         int partialsOffset = 0;
         for (int taxon = 0; taxon < numTaxa; ++taxon) {
--- a/src/dr/evomodel/treedatalikelihood/continuous/RepeatedMeasuresTraitLikelihood.java
+++ b/src/dr/evomodel/treedatalikelihood/continuous/RepeatedMeasuresTraitLikelihood.java
@@ -39,7 +39,7 @@
 import dr.math.matrixAlgebra.WrappedVector;
 import dr.math.matrixAlgebra.missingData.InversionResult;
 import dr.xml.*;
-import org.ejml.data.DenseMatrix64F;
+import org.ejml.data.DMatrixRMaj;
 
 import java.util.ArrayList;
 import java.util.Arrays;
@@ -230,7 +230,7 @@
 
     private final double nuggetPrecision;
 
-//    private void computePrecisionForTaxon(final DenseMatrix64F precision, final int taxon,
+//    private void computePrecisionForTaxon(final DMatrixRMaj precision, final int taxon,
 //                                           final int numFactors) {
 //
 //        final double[] observed = observedIndicators[taxon];
@@ -250,7 +250,7 @@
 //        }
 //    }
 
-//    private InversionResult fillInMeanForTaxon(final WrappedVector output, final DenseMatrix64F precision, final int taxon) {
+//    private InversionResult fillInMeanForTaxon(final WrappedVector output, final DMatrixRMaj precision, final int taxon) {
 //
 //        final double[] observed = observedIndicators[taxon];
 //        final Parameter Y = traitParameter.getParameter(taxon);
@@ -270,8 +270,8 @@
 //            tmp[row] = sum;
 //        }
 //
-//        DenseMatrix64F B = DenseMatrix64F.wrap(numFactors, 1, tmp);
-//        DenseMatrix64F X = DenseMatrix64F.wrap(numFactors, 1, tmp2);
+//        DMatrixRMaj B = DMatrixRMaj.wrap(numFactors, 1, tmp);
+//        DMatrixRMaj X = DMatrixRMaj.wrap(numFactors, 1, tmp2);
 //
 //        InversionResult ci = safeSolve(precision, B, X, true);
 //
@@ -295,7 +295,7 @@
 //        return sum;
 //    }
 
-//    private double computeFactorInnerProduct(final WrappedVector mean, final DenseMatrix64F precision) {
+//    private double computeFactorInnerProduct(final WrappedVector mean, final DMatrixRMaj precision) {
 //        // Compute \mu_i^t P_i \mu^t
 //        double sum = 0;
 //        for (int row = 0; row < numFactors; ++row) {
@@ -320,7 +320,7 @@
 //        return det;
 //    }
 
-//    private void makeCompletedUnobserved(final DenseMatrix64F matrix, double diagonal) {
+//    private void makeCompletedUnobserved(final DMatrixRMaj matrix, double diagonal) {
 //        for (int row = 0; row < numFactors; ++row) {
 //            for (int col = 0; col < numFactors; ++col) {
 //                double x = (row == col) ? diagonal : 0.0;
@@ -331,8 +331,8 @@
 
     private void computePartialsAndRemainders() {
         
-//        final DenseMatrix64F precision = new DenseMatrix64F(numFactors, numFactors);
-//        final DenseMatrix64F variance = new DenseMatrix64F(numFactors, numFactors);
+//        final DMatrixRMaj precision = new DMatrixRMaj(numFactors, numFactors);
+//        final DMatrixRMaj variance = new DMatrixRMaj(numFactors, numFactors);
 //
 //        int partialsOffset = 0;
 //        for (int taxon = 0; taxon < numTaxa; ++taxon) {
--- a/src/dr/evomodel/treedatalikelihood/preorder/AbstractValuesViaFullConditionalDelegate.java
+++ b/src/dr/evomodel/treedatalikelihood/preorder/AbstractValuesViaFullConditionalDelegate.java
@@ -5,7 +5,7 @@
 import dr.evomodel.continuous.MultivariateDiffusionModel;
 import dr.evomodel.treedatalikelihood.continuous.*;
 import dr.math.matrixAlgebra.WrappedVector;
-import org.ejml.data.DenseMatrix64F;
+import org.ejml.data.DMatrixRMaj;
 
 import static dr.math.matrixAlgebra.missingData.MissingOps.wrap;
 
@@ -63,7 +63,7 @@
                 System.err.println("Missing tip = " + node.getNumber() + " (" + nodeBuffer + "), trait = " + trait);
 
                 final WrappedVector preMean = new WrappedVector.Raw(conditionalNodeBuffer, partialOffset, dimTrait);
-                final DenseMatrix64F preVar = wrap(conditionalNodeBuffer, partialOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
+                final DMatrixRMaj preVar = wrap(conditionalNodeBuffer, partialOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
 
                 final WrappedVector postObs = new WrappedVector.Raw(partialNodeBuffer, partialOffset, dimTrait);
 
--- a/src/dr/evomodel/treedatalikelihood/preorder/ConditionalVarianceAndTransform2.java
+++ b/src/dr/evomodel/treedatalikelihood/preorder/ConditionalVarianceAndTransform2.java
@@ -1,8 +1,8 @@
 package dr.evomodel.treedatalikelihood.preorder;
 
 import dr.math.matrixAlgebra.WrappedVector;
-import org.ejml.data.DenseMatrix64F;
-import org.ejml.ops.CommonOps;
+import org.ejml.data.DMatrixRMaj;
+import org.ejml.dense.row.CommonOps_DDRM;
 
 import static dr.math.matrixAlgebra.missingData.MissingOps.gatherRowsAndColumns;
 
@@ -23,8 +23,8 @@
      * \bar{\Sigma} = \Sigma_{11} - \Sigma_{12}\Sigma_{22}^1\Sigma{21}
      */
 
-    final private DenseMatrix64F sBar;
-    final private DenseMatrix64F affineTransform;
+    final private DMatrixRMaj sBar;
+    final private DMatrixRMaj affineTransform;
 
     private final int[] missingIndices;
     private final int[] notMissingIndices;
@@ -36,9 +36,9 @@
     private static final boolean DEBUG = false;
 
     private double[][] cholesky = null;
-    private DenseMatrix64F sBarInv = null;
+    private DMatrixRMaj sBarInv = null;
 
-    ConditionalVarianceAndTransform2(final DenseMatrix64F variance,
+    ConditionalVarianceAndTransform2(final DMatrixRMaj variance,
                                             final int[] missingIndices, final int[] notMissingIndices) {
 
         assert (missingIndices.length + notMissingIndices.length == variance.getNumRows());
@@ -51,44 +51,44 @@
             System.err.println("variance:\n" + variance);
         }
 
-        DenseMatrix64F S22 = new DenseMatrix64F(notMissingIndices.length, notMissingIndices.length);
+        DMatrixRMaj S22 = new DMatrixRMaj(notMissingIndices.length, notMissingIndices.length);
         gatherRowsAndColumns(variance, S22, notMissingIndices, notMissingIndices);
 
         if (DEBUG) {
             System.err.println("S22:\n" + S22);
         }
 
-        DenseMatrix64F S22Inv = new DenseMatrix64F(notMissingIndices.length, notMissingIndices.length);
-        CommonOps.invert(S22, S22Inv);
+        DMatrixRMaj S22Inv = new DMatrixRMaj(notMissingIndices.length, notMissingIndices.length);
+        CommonOps_DDRM.invert(S22, S22Inv);
 
         if (DEBUG) {
             System.err.println("S22Inv:\n" + S22Inv);
         }
 
-        DenseMatrix64F S12 = new DenseMatrix64F(missingIndices.length, notMissingIndices.length);
+        DMatrixRMaj S12 = new DMatrixRMaj(missingIndices.length, notMissingIndices.length);
         gatherRowsAndColumns(variance, S12, missingIndices, notMissingIndices);
 
         if (DEBUG) {
             System.err.println("S12:\n" + S12);
         }
 
-        DenseMatrix64F S12S22Inv = new DenseMatrix64F(missingIndices.length, notMissingIndices.length);
-        CommonOps.mult(S12, S22Inv, S12S22Inv);
+        DMatrixRMaj S12S22Inv = new DMatrixRMaj(missingIndices.length, notMissingIndices.length);
+        CommonOps_DDRM.mult(S12, S22Inv, S12S22Inv);
 
         if (DEBUG) {
             System.err.println("S12S22Inv:\n" + S12S22Inv);
         }
 
-        DenseMatrix64F S12S22InvS21 = new DenseMatrix64F(missingIndices.length, missingIndices.length);
-        CommonOps.multTransB(S12S22Inv, S12, S12S22InvS21);
+        DMatrixRMaj S12S22InvS21 = new DMatrixRMaj(missingIndices.length, missingIndices.length);
+        CommonOps_DDRM.multTransB(S12S22Inv, S12, S12S22InvS21);
 
         if (DEBUG) {
             System.err.println("S12S22InvS21:\n" + S12S22InvS21);
         }
 
-        sBar = new DenseMatrix64F(missingIndices.length, missingIndices.length);
+        sBar = new DMatrixRMaj(missingIndices.length, missingIndices.length);
         gatherRowsAndColumns(variance, sBar, missingIndices, missingIndices);
-        CommonOps.subtract(sBar, S12S22InvS21, sBar);
+        CommonOps_DDRM.subtract(sBar, S12S22InvS21, sBar);
 
 
         if (DEBUG) {
@@ -141,18 +141,18 @@
         return cholesky;
     }
 
-//    public final DenseMatrix64F getAffineTransform() {
+//    public final DMatrixRMaj getAffineTransform() {
 //        return affineTransform;
 //    }
 
-    final DenseMatrix64F getConditionalVariance() {
+    final DMatrixRMaj getConditionalVariance() {
         return sBar;
     }
 
-    final DenseMatrix64F getConditionalPrecision() {
+    final DMatrixRMaj getConditionalPrecision() {
         if (sBarInv == null) {
-            sBarInv = new DenseMatrix64F(numMissing, numMissing);
-            CommonOps.invert(sBar, sBarInv);
+            sBarInv = new DMatrixRMaj(numMissing, numMissing);
+            CommonOps_DDRM.invert(sBar, sBarInv);
         }
         return sBarInv;
     }
--- a/src/dr/evomodel/treedatalikelihood/preorder/MultivariateConditionalOnTipsRealizedDelegate.java
+++ b/src/dr/evomodel/treedatalikelihood/preorder/MultivariateConditionalOnTipsRealizedDelegate.java
@@ -6,8 +6,8 @@
 import dr.math.distributions.MultivariateNormalDistribution;
 import dr.math.matrixAlgebra.Matrix;
 import dr.math.matrixAlgebra.WrappedVector;
-import org.ejml.data.DenseMatrix64F;
-import org.ejml.ops.CommonOps;
+import org.ejml.data.DMatrixRMaj;
+import org.ejml.dense.row.CommonOps_DDRM;
 
 import static dr.math.matrixAlgebra.missingData.MissingOps.*;
 
@@ -41,14 +41,14 @@
         // scalar, dT + 2 * dT * dT, 1
 
         // Integrate out against prior
-        final DenseMatrix64F rootPrec = wrap(partialNodeBuffer, offsetPartial + dimTrait, dimTrait, dimTrait);
-        final DenseMatrix64F priorPrec = new DenseMatrix64F(dimTrait, dimTrait);
-        CommonOps.mult(Pd, wrap(partialPriorBuffer, offsetPartial + dimTrait, dimTrait, dimTrait), priorPrec);
+        final DMatrixRMaj rootPrec = wrap(partialNodeBuffer, offsetPartial + dimTrait, dimTrait, dimTrait);
+        final DMatrixRMaj priorPrec = new DMatrixRMaj(dimTrait, dimTrait);
+        CommonOps_DDRM.mult(Pd, wrap(partialPriorBuffer, offsetPartial + dimTrait, dimTrait, dimTrait), priorPrec);
 
-        final DenseMatrix64F totalPrec = new DenseMatrix64F(dimTrait, dimTrait);
-        CommonOps.add(rootPrec, priorPrec, totalPrec);
+        final DMatrixRMaj totalPrec = new DMatrixRMaj(dimTrait, dimTrait);
+        CommonOps_DDRM.add(rootPrec, priorPrec, totalPrec);
 
-        final DenseMatrix64F totalVar = new DenseMatrix64F(dimTrait, dimTrait);
+        final DMatrixRMaj totalVar = new DMatrixRMaj(dimTrait, dimTrait);
         safeInvert(totalPrec, totalVar, false);
 
         final double[] tmp = new double[dimTrait];
@@ -93,7 +93,7 @@
         }
     }
 
-//    boolean extremeValue(final DenseMatrix64F x) {
+//    boolean extremeValue(final DMatrixRMaj x) {
 //        return extremeValue(new WrappedVector.Raw(x.getData(), 0, x.getNumElements()));
 //    }
 //
@@ -130,7 +130,7 @@
                                               final int offsetPartial,
                                               final double branchPrecision) {
 
-        final DenseMatrix64F P0 = wrap(partialNodeBuffer, offsetPartial + dimTrait, dimTrait, dimTrait);
+        final DMatrixRMaj P0 = wrap(partialNodeBuffer, offsetPartial + dimTrait, dimTrait, dimTrait);
         final int missingCount = countFiniteDiagonals(P0);
 
         if (missingCount == 0) { // Completely observed
@@ -167,26 +167,26 @@
                     final int[] observed = indices.getComplement();
                     final int[] missing = indices.getArray();
 
-                    final DenseMatrix64F V1 = new DenseMatrix64F(dimTrait, dimTrait);
-                    CommonOps.scale(1.0 / branchPrecision, Vd, V1);
+                    final DMatrixRMaj V1 = new DMatrixRMaj(dimTrait, dimTrait);
+                    CommonOps_DDRM.scale(1.0 / branchPrecision, Vd, V1);
 
                     ConditionalVarianceAndTransform2 transform =
                             new ConditionalVarianceAndTransform2(
                                     V1, missing, observed
                             ); // TODO Cache (via delegated function)
 
-                    final DenseMatrix64F cP0 = new DenseMatrix64F(missing.length, missing.length);
+                    final DMatrixRMaj cP0 = new DMatrixRMaj(missing.length, missing.length);
                     gatherRowsAndColumns(P0, cP0, missing, missing);
 
                     final WrappedVector cM2 = transform.getConditionalMean(
                             partialNodeBuffer, offsetPartial, // Tip value
                             sample, offsetParent); // Parent value
 
-                    final DenseMatrix64F cP1 = transform.getConditionalPrecision();
+                    final DMatrixRMaj cP1 = transform.getConditionalPrecision();
 
-                    final DenseMatrix64F cP2 = new DenseMatrix64F(missing.length, missing.length);
-                    final DenseMatrix64F cV2 = new DenseMatrix64F(missing.length, missing.length);
-                    CommonOps.add(cP0, cP1, cP2);
+                    final DMatrixRMaj cP2 = new DMatrixRMaj(missing.length, missing.length);
+                    final DMatrixRMaj cV2 = new DMatrixRMaj(missing.length, missing.length);
+                    CommonOps_DDRM.add(cP0, cP1, cP2);
 
                     safeInvert(cP2, cV2, false);
                     double[][] cC2 = getCholeskyOfVariance(cV2.getData(), missing.length);
@@ -205,8 +205,8 @@
                         final WrappedVector M0 = new WrappedVector.Raw(partialNodeBuffer, offsetPartial, dimTrait);
 
                         final WrappedVector M1 = new WrappedVector.Raw(sample, offsetParent, dimTrait);
-                        final DenseMatrix64F P1 = new DenseMatrix64F(dimTrait, dimTrait);
-                        CommonOps.scale(branchPrecision, Pd, P1);
+                        final DMatrixRMaj P1 = new DMatrixRMaj(dimTrait, dimTrait);
+                        CommonOps_DDRM.scale(branchPrecision, Pd, P1);
 
                         final WrappedVector newSample = new WrappedVector.Raw(sample, offsetSample, dimTrait);
 
@@ -246,25 +246,25 @@
         if (!Double.isInfinite(branchPrecision)) {
 
             final WrappedVector M0 = new WrappedVector.Raw(partialNodeBuffer, offsetPartial, dimTrait);
-            final DenseMatrix64F P0 = wrap(partialNodeBuffer, offsetPartial + dimTrait, dimTrait, dimTrait);
+            final DMatrixRMaj P0 = wrap(partialNodeBuffer, offsetPartial + dimTrait, dimTrait, dimTrait);
 
             final WrappedVector M1;
-            final DenseMatrix64F P1;
+            final DMatrixRMaj P1;
 
             if (hasNoDrift) {
                 M1 = new WrappedVector.Raw(sample, offsetParent, dimTrait);
-                P1 = new DenseMatrix64F(dimTrait, dimTrait);
-                CommonOps.scale(branchPrecision, Pd, P1);
+                P1 = new DMatrixRMaj(dimTrait, dimTrait);
+                CommonOps_DDRM.scale(branchPrecision, Pd, P1);
             } else {
                 M1 = getMeanWithDrift(sample, offsetParent, displacementBuffer, dimTrait);
-                P1 = DenseMatrix64F.wrap(dimTrait, dimTrait, precisionBuffer);
+                P1 = DMatrixRMaj.wrap(dimTrait, dimTrait, precisionBuffer);
             }
 
 //            boolean DEBUG_PRECISION = false;
 //
 //            if (DEBUG_PRECISION) {
-//                DenseMatrix64F tP1 = new DenseMatrix64F(dimTrait, dimTrait);
-//                CommonOps.scale(branchPrecision, Pd, tP1);
+//                DMatrixRMaj tP1 = new DMatrixRMaj(dimTrait, dimTrait);
+//                CommonOps_DDRM.scale(branchPrecision, Pd, tP1);
 //
 //                for (int i = 0; i < dimTrait; ++i) {
 //                    for (int j = 0; j < dimTrait; ++j) {
@@ -277,10 +277,10 @@
 //            }
 
             final WrappedVector M2 = new WrappedVector.Raw(tmpMean, 0, dimTrait);
-            final DenseMatrix64F P2 = new DenseMatrix64F(dimTrait, dimTrait);
-            final DenseMatrix64F V2 = new DenseMatrix64F(dimTrait, dimTrait);
+            final DMatrixRMaj P2 = new DMatrixRMaj(dimTrait, dimTrait);
+            final DMatrixRMaj V2 = new DMatrixRMaj(dimTrait, dimTrait);
 
-            CommonOps.add(P0, P1, P2);
+            CommonOps_DDRM.add(P0, P1, P2);
             safeInvert(P2, V2, false);
             weightedAverage(M0, P0, M1, P1, M2, V2, dimTrait);
 
--- a/src/dr/evomodel/treedatalikelihood/preorder/ProcessSimulationDelegate.java
+++ b/src/dr/evomodel/treedatalikelihood/preorder/ProcessSimulationDelegate.java
@@ -35,7 +35,7 @@
 import dr.inference.model.Model;
 import dr.inference.model.ModelListener;
 import dr.math.matrixAlgebra.*;
-import org.ejml.data.DenseMatrix64F;
+import org.ejml.data.DMatrixRMaj;
 
 import java.util.List;
 import java.util.Map;
@@ -161,8 +161,8 @@
         final ContinuousDataLikelihoodDelegate likelihoodDelegate;
 
         double[] diffusionVariance;
-        DenseMatrix64F Vd;
-        DenseMatrix64F Pd;
+        DMatrixRMaj Vd;
+        DMatrixRMaj Pd;
 
         double[][] cholesky;
         Map<PartiallyMissingInformation.HashedIntArray,
@@ -221,7 +221,7 @@
                 double[][] diffusionPrecision = diffusionModel.getPrecisionmatrix();
                 diffusionVariance = getVectorizedVarianceFromPrecision(diffusionPrecision);
                 Vd = wrap(diffusionVariance, 0, dimTrait, dimTrait);
-                Pd = new DenseMatrix64F(diffusionPrecision);
+                Pd = new DMatrixRMaj(diffusionPrecision);
             }
             if (cholesky == null) {
                 cholesky = getCholeskyOfVariance(diffusionVariance, dimTrait);
--- a/src/dr/inferencexml/operators/EllipticalSliceOperatorParser.java
+++ b/src/dr/inferencexml/operators/EllipticalSliceOperatorParser.java
@@ -35,7 +35,7 @@
 import dr.math.distributions.GaussianProcessRandomGenerator;
 import dr.math.distributions.MultivariateNormalDistribution;
 import dr.xml.*;
-import org.ejml.ops.CommonOps;
+import org.ejml.dense.row.CommonOps_DDRM;
 
 /**
  */
--- a/src/dr/math/matrixAlgebra/missingData/MissingOps.java
+++ b/src/dr/math/matrixAlgebra/missingData/MissingOps.java
@@ -3,17 +3,17 @@
 import dr.math.matrixAlgebra.Vector;
 import dr.math.matrixAlgebra.WrappedVector;
 import dr.util.Transform;
-import org.ejml.alg.dense.decomposition.lu.LUDecompositionAlt_D64;
-import org.ejml.alg.dense.decomposition.qr.QRColPivDecompositionHouseholderColumn_D64;
-import org.ejml.alg.dense.linsol.lu.LinearSolverLu_D64;
-import org.ejml.alg.dense.misc.UnrolledDeterminantFromMinor;
-import org.ejml.alg.dense.misc.UnrolledInverseFromMinor;
-import org.ejml.data.DenseMatrix64F;
-import org.ejml.factory.LinearSolverFactory;
-import org.ejml.interfaces.decomposition.QRPDecomposition;
-import org.ejml.interfaces.decomposition.SingularValueDecomposition;
+import org.ejml.dense.row.decomposition.lu.LUDecompositionAlt_DDRM;
+import org.ejml.dense.row.decomposition.qr.QRColPivDecompositionHouseholderColumn_DDRM;
+import org.ejml.dense.row.linsol.lu.LinearSolverLu_DDRM;
+import org.ejml.dense.row.misc.UnrolledDeterminantFromMinor_DDRM;
+import org.ejml.dense.row.misc.UnrolledInverseFromMinor_DDRM;
+import org.ejml.data.DMatrixRMaj;
+import org.ejml.dense.row.factory.LinearSolverFactory_DDRM;
+import org.ejml.interfaces.decomposition.QRPDecomposition_F64;
+import org.ejml.interfaces.decomposition.SingularValueDecomposition_F64;
 import org.ejml.interfaces.linsol.LinearSolver;
-import org.ejml.ops.CommonOps;
+import org.ejml.dense.row.CommonOps_DDRM;
 
 import java.util.Arrays;
 
@@ -26,21 +26,21 @@
  */
 public class MissingOps {
 
-    public static DenseMatrix64F wrap(final double[] source, final int offset,
+    public static DMatrixRMaj wrap(final double[] source, final int offset,
                                       final int numRows, final int numCols) {
         double[] buffer = new double[numRows * numCols];
         return wrap(source, offset, numRows, numCols, buffer);
     }
 
-    public static DenseMatrix64F wrap(final double[] source, final int offset,
+    public static DMatrixRMaj wrap(final double[] source, final int offset,
                                               final int numRows, final int numCols,
                                               final double[] buffer) {
         System.arraycopy(source, offset, buffer, 0, numRows * numCols);
-        return DenseMatrix64F.wrap(numRows, numCols, buffer);
+        return DMatrixRMaj.wrap(numRows, numCols, buffer);
     }
 
 
-    public static void gatherRowsAndColumns(final DenseMatrix64F source, final DenseMatrix64F destination,
+    public static void gatherRowsAndColumns(final DMatrixRMaj source, final DMatrixRMaj destination,
                                                       final int[] rowIndices, final int[] colIndices) {
 
         final int rowLength = rowIndices.length;
@@ -57,7 +57,7 @@
         }
     }
 
-    public static void scatterRowsAndColumns(final DenseMatrix64F source, final DenseMatrix64F destination,
+    public static void scatterRowsAndColumns(final DMatrixRMaj source, final DMatrixRMaj destination,
                                              final int[] rowIdices, final int[] colIndices, final boolean clear) {
         if (clear) {
             Arrays.fill(destination.getData(), 0.0);
@@ -77,11 +77,11 @@
         }
     }
 
-    public static void unwrap(final DenseMatrix64F source, final double[] destination, final int offset) {
+    public static void unwrap(final DMatrixRMaj source, final double[] destination, final int offset) {
         System.arraycopy(source.getData(), 0, destination, offset, source.getNumElements());
     }
 
-    public static boolean anyDiagonalInfinities(DenseMatrix64F source) {
+    public static boolean anyDiagonalInfinities(DMatrixRMaj source) {
         boolean anyInfinities = false;
         for (int i = 0; i < source.getNumCols() && !anyInfinities; ++i) {
             if (Double.isInfinite(source.unsafe_get(i, i))) {
@@ -91,7 +91,7 @@
         return anyInfinities;
     }
 
-    public static boolean allFiniteDiagonals(DenseMatrix64F source) {
+    public static boolean allFiniteDiagonals(DMatrixRMaj source) {
         boolean allFinite = true;
 
         final int length = source.getNumCols();
@@ -101,7 +101,7 @@
         return allFinite;
     }
 
-    public static int countFiniteDiagonals(DenseMatrix64F source) {
+    public static int countFiniteDiagonals(DMatrixRMaj source) {
         final int length = source.getNumCols();
 
         int count = 0;
@@ -114,7 +114,7 @@
         return count;
     }
 
-    public static int countZeroDiagonals(DenseMatrix64F source) {
+    public static int countZeroDiagonals(DMatrixRMaj source) {
         final int length = source.getNumCols();
 
         int count = 0;
@@ -127,7 +127,7 @@
         return count;
     }
 
-    public static void getFiniteDiagonalIndices(final DenseMatrix64F source, final int[] indices) {
+    public static void getFiniteDiagonalIndices(final DMatrixRMaj source, final int[] indices) {
         final int length = source.getNumCols();
 
         int index = 0;
@@ -140,7 +140,7 @@
         }
     }
 
-    public static int countFiniteNonZeroDiagonals(DenseMatrix64F source) {
+    public static int countFiniteNonZeroDiagonals(DMatrixRMaj source) {
         final int length = source.getNumCols();
 
         int count = 0;
@@ -153,7 +153,7 @@
         return count;
     }
 
-    public static void getFiniteNonZeroDiagonalIndices(final DenseMatrix64F source, final int[] indices) {
+    public static void getFiniteNonZeroDiagonalIndices(final DMatrixRMaj source, final int[] indices) {
         final int length = source.getNumCols();
 
         int index = 0;
@@ -166,22 +166,22 @@
         }
     }
 
-    public static void addToDiagonal(DenseMatrix64F source, double increment) {
+    public static void addToDiagonal(DMatrixRMaj source, double increment) {
         final int width = source.getNumRows();
         for (int i = 0; i < width; ++i) {
             source.unsafe_set(i,i, source.unsafe_get(i, i) + increment);
         }
     }
 
-    public static double det(DenseMatrix64F mat) {
+    public static double det(DMatrixRMaj mat) {
         int numCol = mat.getNumCols();
         int numRow = mat.getNumRows();
         if(numCol != numRow) {
             throw new IllegalArgumentException("Must be a square matrix.");
         } else if(numCol <= 6) {
-            return numCol >= 2? UnrolledDeterminantFromMinor.det(mat):mat.get(0);
+            return numCol >= 2? UnrolledDeterminantFromMinor_DDRM.det(mat):mat.get(0);
         } else {
-            LUDecompositionAlt_D64 alg = new LUDecompositionAlt_D64();
+            LUDecompositionAlt_DDRM alg = new LUDecompositionAlt_DDRM();
             if(alg.inputModified()) {
                 mat = mat.copy();
             }
@@ -190,7 +190,7 @@
         }
     }
 
-    public static double invertAndGetDeterminant(DenseMatrix64F mat, DenseMatrix64F result) {
+    public static double invertAndGetDeterminant(DMatrixRMaj mat, DMatrixRMaj result) {
 
         final int numCol = mat.getNumCols();
         final int numRow = mat.getNumRows();
@@ -201,18 +201,18 @@
         if (numCol <= 5) {
 
             if (numCol >= 2) {
-                UnrolledInverseFromMinor.inv(mat, result);
+                UnrolledInverseFromMinor_DDRM.inv(mat, result);
             } else {
                 result.set(0, 1.0D / mat.get(0));
             }
 
             return numCol >= 2 ?
-                    UnrolledDeterminantFromMinor.det(mat) :
+                    UnrolledDeterminantFromMinor_DDRM.det(mat) :
                     mat.get(0);
         } else {
 
-            LUDecompositionAlt_D64 alg = new LUDecompositionAlt_D64();
-            LinearSolverLu_D64 solver = new LinearSolverLu_D64(alg);
+            LUDecompositionAlt_DDRM alg = new LUDecompositionAlt_DDRM();
+            LinearSolverLu_DDRM solver = new LinearSolverLu_DDRM(alg);
             if (solver.modifiesA()) {
                 mat = mat.copy();
             }
@@ -228,7 +228,7 @@
         }
     }
 
-    public static InversionResult safeDeterminant(DenseMatrix64F source, boolean invert) {
+    public static InversionResult safeDeterminant(DMatrixRMaj source, boolean invert) {
         final int finiteCount = countFiniteNonZeroDiagonals(source);
 
         InversionResult result;
@@ -236,10 +236,10 @@
         if (finiteCount == 0) {
             result = new InversionResult(NOT_OBSERVED, 0, 0);
         } else {
-            LinearSolver<DenseMatrix64F> solver = LinearSolverFactory.pseudoInverse(true);
+            LinearSolver<DMatrixRMaj, DMatrixRMaj> solver = LinearSolverFactory_DDRM.pseudoInverse(true);
             solver.setA(source);
 
-            SingularValueDecomposition<DenseMatrix64F> svd = solver.getDecomposition();
+            SingularValueDecomposition_F64<DMatrixRMaj> svd = solver.getDecomposition();
             double[] values = svd.getSingularValues();
 
 
@@ -267,7 +267,7 @@
         return result;
     }
 
-    public static InversionResult safeSolve(DenseMatrix64F A,
+    public static InversionResult safeSolve(DMatrixRMaj A,
                                             WrappedVector b,
                                             WrappedVector x,
                                             boolean getDeterminat) {
@@ -275,8 +275,8 @@
 
         assert(A.getNumRows() == dim && A.getNumCols() == dim);
 
-        final DenseMatrix64F B = wrap(b.getBuffer(), b.getOffset(), dim, 1);
-        final DenseMatrix64F X = new DenseMatrix64F(dim, 1);
+        final DMatrixRMaj B = wrap(b.getBuffer(), b.getOffset(), dim, 1);
+        final DMatrixRMaj X = new DMatrixRMaj(dim, 1);
 
         InversionResult ir = safeSolve(A, B, X, getDeterminat);
 
@@ -288,7 +288,7 @@
         return ir;
     }
 
-    public static InversionResult safeSolve(DenseMatrix64F A, DenseMatrix64F B, DenseMatrix64F X, boolean getDeterminant) {
+    public static InversionResult safeSolve(DMatrixRMaj A, DMatrixRMaj B, DMatrixRMaj X, boolean getDeterminant) {
 
         final int finiteCount = countFiniteNonZeroDiagonals(A);
 
@@ -298,7 +298,7 @@
             result = new InversionResult(NOT_OBSERVED, 0, 0);
         } else {
 
-            LinearSolver<DenseMatrix64F> solver = LinearSolverFactory.pseudoInverse(true);
+            LinearSolver<DMatrixRMaj, DMatrixRMaj> solver = LinearSolverFactory_DDRM.pseudoInverse(true);
             solver.setA(A);
             solver.solve(B, X);
 
@@ -306,7 +306,7 @@
             double det = 1;
 
             if (getDeterminant) {
-                SingularValueDecomposition<DenseMatrix64F> svd = solver.getDecomposition();
+                SingularValueDecomposition_F64<DMatrixRMaj> svd = solver.getDecomposition();
                 double[] values = svd.getSingularValues();
 
                 for (int i = 0; i < values.length; ++i) {
@@ -325,7 +325,7 @@
     }
 
 
-    public static InversionResult safeInvert(DenseMatrix64F source, DenseMatrix64F destination, boolean getDeterminant) {
+    public static InversionResult safeInvert(DMatrixRMaj source, DMatrixRMaj destination, boolean getDeterminant) {
 
         final int dim = source.getNumCols();
         final int finiteCount = countFiniteNonZeroDiagonals(source);
@@ -335,7 +335,7 @@
             if (getDeterminant) {
                 det = invertAndGetDeterminant(source, destination);
             } else {
-                CommonOps.invert(source, destination);
+                CommonOps_DDRM.invert(source, destination);
             }
             return new InversionResult(FULLY_OBSERVED, dim, det);
         } else {
@@ -346,14 +346,14 @@
                 final int[] finiteIndices = new int[finiteCount];
                 getFiniteNonZeroDiagonalIndices(source, finiteIndices);
 
-                final DenseMatrix64F subSource = new DenseMatrix64F(finiteCount, finiteCount);
+                final DMatrixRMaj subSource = new DMatrixRMaj(finiteCount, finiteCount);
                 gatherRowsAndColumns(source, subSource, finiteIndices, finiteIndices);
 
-                final DenseMatrix64F inverseSubSource = new DenseMatrix64F(finiteCount, finiteCount);
+                final DMatrixRMaj inverseSubSource = new DMatrixRMaj(finiteCount, finiteCount);
                 if (getDeterminant) {
                     det = invertAndGetDeterminant(subSource, inverseSubSource);
                 } else {
-                    CommonOps.invert(subSource, inverseSubSource);
+                    CommonOps_DDRM.invert(subSource, inverseSubSource);
                 }
 
                 scatterRowsAndColumns(inverseSubSource, destination, finiteIndices, finiteIndices, true);
@@ -364,11 +364,11 @@
     }
 
     public static void weightedAverage(final WrappedVector mi,
-                                       final DenseMatrix64F Pi,
+                                       final DMatrixRMaj Pi,
                                        final WrappedVector mj,
-                                       final DenseMatrix64F Pj,
+                                       final DMatrixRMaj Pj,
                                        final WrappedVector mk,
-                                       final DenseMatrix64F Vk,
+                                       final DMatrixRMaj Vk,
                                        final int dimTrait) {
         final double[] tmp = new double[dimTrait];
         for (int g = 0; g < dimTrait; ++g) {
@@ -390,13 +390,13 @@
 
     public static void weightedAverage(final double[] ipartial,
                                        final int ibo,
-                                       final DenseMatrix64F Pi,
+                                       final DMatrixRMaj Pi,
                                        final double[] jpartial,
                                        final int jbo,
-                                       final DenseMatrix64F Pj,
+                                       final DMatrixRMaj Pj,
                                        final double[] kpartial,
                                        final int kbo,
-                                       final DenseMatrix64F Vk,
+                                       final DMatrixRMaj Vk,
                                        final int dimTrait) {
         final double[] tmp = new double[dimTrait];
         for (int g = 0; g < dimTrait; ++g) {