diff --git a/cpp/Matrix.cpp b/cpp/Matrix.cpp
index 765b68586..cabcb8210 100644
--- a/cpp/Matrix.cpp
+++ b/cpp/Matrix.cpp
@@ -1021,9 +1021,9 @@ Matrix inverse_square_root(const Matrix& A) {
 // inv(B) * inv(B)' == A
 // inv(B' * B) == A
 Matrix inverse_square_root(const Matrix& A) {
-	Matrix L = LLt(A);
+	Matrix R = RtR(A);
         Matrix inv(boost::numeric::ublas::identity_matrix<double>(A.size1()));
-        inplace_solve (L, inv, BNU::lower_tag ());
+        inplace_solve (R, inv, BNU::upper_tag ());
 	return inv;
 }
 
diff --git a/cpp/testMatrix.cpp b/cpp/testMatrix.cpp
index 1eaf6995a..3f6a5801e 100644
--- a/cpp/testMatrix.cpp
+++ b/cpp/testMatrix.cpp
@@ -817,18 +817,19 @@ TEST( matrix, inverse_square_root )
 	// use the same inverse routing inside inverse_square_root()
 	// as we use here to check it.
 
-	Matrix M = Matrix_(5, 5, 0.0785892, 0.0137923, -0.0142219, -0.0171880,
-			0.0028726, 0.0137923, 0.0908911, 0.0020775, -0.0101952, 0.0175868,
-			-0.0142219, 0.0020775, 0.0973051, 0.0054906, 0.0047064, -0.0171880,
-			-0.0101952, 0.0054906, 0.0892453, -0.0059468, 0.0028726, 0.0175868,
-			0.0047064, -0.0059468, 0.0816517);
+	Matrix M = Matrix_(5, 5, 
+			   0.0785892, 0.0137923, -0.0142219, -0.0171880, 0.0028726, 
+                           0.0137923, 0.0908911, 0.0020775, -0.0101952, 0.0175868,
+		  	  -0.0142219, 0.0020775, 0.0973051, 0.0054906, 0.0047064, 
+                          -0.0171880,-0.0101952, 0.0054906, 0.0892453, -0.0059468, 
+			   0.0028726, 0.0175868, 0.0047064, -0.0059468, 0.0816517);
 
-	expected = Matrix_(5, 5,
-			3.567126953241796, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000,
-			-0.590030436566913, 3.362022286742925, 0.000000000000000, 0.000000000000000, 0.000000000000000,
-			0.618207860252376, -0.168166020746503, 3.253086082942785, 0.000000000000000, 0.000000000000000,
-			0.683045380655496, 0.283773848115276, -0.099969232183396, 3.433537147891568, 0.000000000000000,
-			-0.006740136923185, -0.669325697387650, -0.169716689114923, 0.171493059476284, 3.583921085468937);
+	expected = trans(Matrix_(5, 5,
+				 3.567126953241796, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000,
+				 -0.590030436566913, 3.362022286742925, 0.000000000000000, 0.000000000000000, 0.000000000000000,
+				 0.618207860252376, -0.168166020746503, 3.253086082942785, 0.000000000000000, 0.000000000000000,
+				 0.683045380655496, 0.283773848115276, -0.099969232183396, 3.433537147891568, 0.000000000000000,
+				 -0.006740136923185, -0.669325697387650, -0.169716689114923, 0.171493059476284, 3.583921085468937));
 	EQUALITY(expected, inverse_square_root(M));
 
 }
@@ -845,11 +846,11 @@ TEST( matrix, LLt )
 			-0.0021113, 0.0036415, 0.0909464);
 
 	Matrix expected = Matrix_(5, 5,
-			0.295668226226627, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000,
-		 -0.010437374483502, 0.295235094820875, 0.000000000000000, 0.000000000000000, 0.000000000000000,
-			0.039560896175007, 0.063407813693827, 0.301721866387571, 0.000000000000000, 0.000000000000000,
-			0.027552165831157, 0.043423266737274, 0.021695600982708, 0.267613525371710, 0.000000000000000,
-			0.016485031422565, -0.012072546984405, -0.006621889326331, 0.014405837566082, 0.300462176944247);
+					0.295668226226627, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000,
+					-0.010437374483502, 0.295235094820875, 0.000000000000000, 0.000000000000000, 0.000000000000000,
+					0.039560896175007, 0.063407813693827, 0.301721866387571, 0.000000000000000, 0.000000000000000,
+					0.027552165831157, 0.043423266737274, 0.021695600982708, 0.267613525371710, 0.000000000000000,
+					0.016485031422565, -0.012072546984405, -0.006621889326331, 0.014405837566082, 0.300462176944247);
 
 	EQUALITY(expected, LLt(M));
 }