diff --git a/.cproject b/.cproject
index d5a6ca4d4..1783efa98 100644
--- a/.cproject
+++ b/.cproject
@@ -568,7 +568,6 @@
 			
 			
 				make
-				
 				tests/testBayesTree.run
 				true
 				false
@@ -576,7 +575,6 @@
 			
 			
 				make
-				
 				testBinaryBayesNet.run
 				true
 				false
@@ -624,7 +622,6 @@
 			
 			
 				make
-				
 				testSymbolicBayesNet.run
 				true
 				false
@@ -632,7 +629,6 @@
 			
 			
 				make
-				
 				tests/testSymbolicFactor.run
 				true
 				false
@@ -640,7 +636,6 @@
 			
 			
 				make
-				
 				testSymbolicFactorGraph.run
 				true
 				false
@@ -656,12 +651,19 @@
 			
 			
 				make
-				
 				tests/testBayesTree
 				true
 				false
 				true
 			
+			
+				make
+				-j2
+				testGaussianFactor.run
+				true
+				true
+				true
+			
 			
 				make
 				-j5
@@ -678,6 +680,142 @@
 				true
 				true
 			
+			
+				make
+				-j5
+				testCal3Bundler.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testCal3DS2.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testCalibratedCamera.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testEssentialMatrix.run
+				true
+				true
+				true
+			
+			
+				make
+				-j1 VERBOSE=1
+				testHomography2.run
+				true
+				false
+				true
+			
+			
+				make
+				-j5
+				testPinholeCamera.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testPoint2.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testPoint3.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testPose2.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testPose3.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testRot3M.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testSphere2.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testStereoCamera.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				timeCalibratedCamera.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				timePinholeCamera.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				timeStereoCamera.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testCal3Unified.run
+				true
+				true
+				true
+			
 			
 				make
 				-j5
@@ -718,6 +856,62 @@
 				true
 				true
 			
+			
+				make
+				-j2
+				all
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				check
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				clean
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				all
+				true
+				false
+				true
+			
+			
+				make
+				-j5
+				check
+				true
+				false
+				true
+			
+			
+				make
+				-j2
+				all
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				clean
+				true
+				true
+				true
+			
 			
 				make
 				-j2
@@ -734,6 +928,30 @@
 				true
 				true
 			
+			
+				make
+				-j2
+				all
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				clean
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				clean all
+				true
+				true
+				true
+			
 			
 				make
 				-j2
@@ -806,6 +1024,46 @@
 				true
 				true
 			
+			
+				make
+				-j5
+				testWrap.run
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				check
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				timeCalibratedCamera.run
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				timeRot3.run
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				clean
+				true
+				true
+				true
+			
 			
 				make
 				-j5
@@ -974,6 +1232,307 @@
 				true
 				true
 			
+			
+				make
+				-j5
+				testDiscreteFactor.run
+				true
+				true
+				true
+			
+			
+				make
+				-j1
+				testDiscreteBayesTree.run
+				true
+				false
+				true
+			
+			
+				make
+				-j5
+				testDiscreteFactorGraph.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testDiscreteConditional.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testDiscreteMarginals.run
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				vSFMexample.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testMatrix.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testVector.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testInvDepthCamera3.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testTriangulation.run
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				testVSLAMGraph
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				check.tests
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				timeGaussianFactorGraph.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testMarginals.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testGaussianISAM2.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testSymbolicFactorGraphB.run
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				timeSequentialOnDataset.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testGradientDescentOptimizer.run
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				testGaussianFactor.run
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				testNonlinearOptimizer.run
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				testGaussianBayesNet.run
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				testNonlinearISAM.run
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				testNonlinearEquality.run
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				testExtendedKalmanFilter.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				timing.tests
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testNonlinearFactor.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				clean
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testGaussianJunctionTreeB.run
+				true
+				true
+				true
+			
+			
+				make
+				testGraph.run
+				true
+				false
+				true
+			
+			
+				make
+				testJunctionTree.run
+				true
+				false
+				true
+			
+			
+				make
+				testSymbolicBayesNetB.run
+				true
+				false
+				true
+			
+			
+				make
+				-j5
+				testGaussianISAM.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testDoglegOptimizer.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testNonlinearFactorGraph.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testIterative.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testSubgraphSolver.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testGaussianFactorGraphB.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testSummarization.run
+				true
+				true
+				true
+			
 			
 				make
 				-j5
@@ -1080,7 +1639,6 @@
 			
 			
 				make
-				
 				testErrors.run
 				true
 				false
@@ -1126,6 +1684,14 @@
 				true
 				true
 			
+			
+				make
+				-j2
+				testGaussianFactor.run
+				true
+				true
+				true
+			
 			
 				make
 				-j5
@@ -1310,6 +1876,54 @@
 				true
 				true
 			
+			
+				make
+				-j2
+				install
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				clean
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				all
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				clean
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				all
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				clean
+				true
+				true
+				true
+			
 			
 				make
 				-j2
@@ -1318,6 +1932,70 @@
 				true
 				true
 			
+			
+				make
+				-j5
+				testAntiFactor.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testBetweenFactor.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testDataset.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testEssentialMatrixFactor.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testGeneralSFMFactor_Cal3Bundler.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testGeneralSFMFactor.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testProjectionFactor.run
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				testRotateFactor.run
+				true
+				true
+				true
+			
 			
 				make
 				-j5
@@ -1416,6 +2094,7 @@
 			
 			
 				make
+				
 				testSimulated2DOriented.run
 				true
 				false
@@ -1455,6 +2134,7 @@
 			
 			
 				make
+				
 				testSimulated2D.run
 				true
 				false
@@ -1462,6 +2142,7 @@
 			
 			
 				make
+				
 				testSimulated3D.run
 				true
 				false
@@ -1515,1140 +2196,6 @@
 				false
 				true
 			
-			
-				make
-				-j2
-				check
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				tests/testLieConfig.run
-				true
-				true
-				true
-			
-			
-				make
-				-j3
-				install
-				true
-				false
-				true
-			
-			
-				make
-				-j2
-				clean
-				true
-				true
-				true
-			
-			
-				make
-				-j1
-				check
-				true
-				false
-				true
-			
-			
-				make
-				-j5
-				all
-				true
-				true
-				true
-			
-			
-				cmake
-				..
-				true
-				false
-				true
-			
-			
-				make
-				-j5
-				gtsam-shared
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				gtsam-static
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				timing
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				examples
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				VERBOSE=1 all
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				VERBOSE=1 check
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				check.base
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				timing.base
-				true
-				true
-				true
-			
-			
-				make
-				-j2 VERBOSE=1
-				check.geometry
-				true
-				false
-				true
-			
-			
-				make
-				-j5
-				timing.geometry
-				true
-				true
-				true
-			
-			
-				make
-				-j2 VERBOSE=1
-				check.inference
-				true
-				false
-				true
-			
-			
-				make
-				-j5
-				timing.inference
-				true
-				true
-				true
-			
-			
-				make
-				-j2 VERBOSE=1
-				check.linear
-				true
-				false
-				true
-			
-			
-				make
-				-j5
-				timing.linear
-				true
-				true
-				true
-			
-			
-				make
-				-j2 VERBOSE=1
-				check.nonlinear
-				true
-				false
-				true
-			
-			
-				make
-				-j5
-				timing.nonlinear
-				true
-				true
-				true
-			
-			
-				make
-				-j2 VERBOSE=1
-				check.slam
-				true
-				false
-				true
-			
-			
-				make
-				-j5
-				timing.slam
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				wrap_gtsam
-				true
-				true
-				true
-			
-			
-				make
-				VERBOSE=1
-				wrap_gtsam
-				true
-				false
-				true
-			
-			
-				cpack
-				
-				-G DEB
-				true
-				false
-				true
-			
-			
-				cpack
-				
-				-G RPM
-				true
-				false
-				true
-			
-			
-				cpack
-				
-				-G TGZ
-				true
-				false
-				true
-			
-			
-				cpack
-				
-				--config CPackSourceConfig.cmake
-				true
-				false
-				true
-			
-			
-				make
-				-j5
-				check.discrete
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				check.discrete_unstable
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				check.base_unstable
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				check.dynamics_unstable
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				check.slam_unstable
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				check.unstable
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				wrap_gtsam_build
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				wrap_gtsam_unstable_build
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				wrap_gtsam_unstable
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				wrap
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				wrap_gtsam_distclean
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				wrap_gtsam_unstable_distclean
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				doc
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				doc_clean
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				check
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				check.geometry_unstable
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				check.linear_unstable
-				true
-				true
-				true
-			
-			
-				make
-				-j6 -j8
-				gtsam_unstable-shared
-				true
-				true
-				true
-			
-			
-				make
-				-j6 -j8
-				gtsam_unstable-static
-				true
-				true
-				true
-			
-			
-				make
-				-j6 -j8
-				check.nonlinear_unstable
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				check.tests
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				check
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				clean
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				install
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				all
-				true
-				true
-				true
-			
-			
-				cmake
-				..
-				true
-				false
-				true
-			
-			
-				make
-				-j2
-				testGaussianFactor.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testCal3Bundler.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testCal3DS2.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testCalibratedCamera.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testEssentialMatrix.run
-				true
-				true
-				true
-			
-			
-				make
-				-j1 VERBOSE=1
-				testHomography2.run
-				true
-				false
-				true
-			
-			
-				make
-				-j5
-				testPinholeCamera.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testPoint2.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testPoint3.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testPose2.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testPose3.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testRot3M.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testSphere2.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testStereoCamera.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				timeCalibratedCamera.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				timePinholeCamera.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				timeStereoCamera.run
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				all
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				check
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				clean
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				all
-				true
-				false
-				true
-			
-			
-				make
-				-j5
-				check
-				true
-				false
-				true
-			
-			
-				make
-				-j2
-				all
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				clean
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				all
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				clean
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				clean all
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				check
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				timeCalibratedCamera.run
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				timeRot3.run
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				clean
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testWrap.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testDiscreteFactor.run
-				true
-				true
-				true
-			
-			
-				make
-				-j1
-				testDiscreteBayesTree.run
-				true
-				false
-				true
-			
-			
-				make
-				-j5
-				testDiscreteFactorGraph.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testDiscreteConditional.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testDiscreteMarginals.run
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				vSFMexample.run
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				testVSLAMGraph
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testInvDepthCamera3.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testTriangulation.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testMatrix.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testVector.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				check.tests
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				timeGaussianFactorGraph.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testMarginals.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testGaussianISAM2.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testSymbolicFactorGraphB.run
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				timeSequentialOnDataset.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testGradientDescentOptimizer.run
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				testGaussianFactor.run
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				testNonlinearOptimizer.run
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				testGaussianBayesNet.run
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				testNonlinearISAM.run
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				testNonlinearEquality.run
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				testExtendedKalmanFilter.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				timing.tests
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testNonlinearFactor.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				clean
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testGaussianJunctionTreeB.run
-				true
-				true
-				true
-			
-			
-				make
-				
-				testGraph.run
-				true
-				false
-				true
-			
-			
-				make
-				
-				testJunctionTree.run
-				true
-				false
-				true
-			
-			
-				make
-				
-				testSymbolicBayesNetB.run
-				true
-				false
-				true
-			
-			
-				make
-				-j5
-				testGaussianISAM.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testDoglegOptimizer.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testNonlinearFactorGraph.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testIterative.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testSubgraphSolver.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testGaussianFactorGraphB.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testSummarization.run
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				testGaussianFactor.run
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				install
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				clean
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				all
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				clean
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				all
-				true
-				true
-				true
-			
-			
-				make
-				-j2
-				clean
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testAntiFactor.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testBetweenFactor.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testDataset.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testEssentialMatrixFactor.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testGeneralSFMFactor_Cal3Bundler.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testGeneralSFMFactor.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testProjectionFactor.run
-				true
-				true
-				true
-			
-			
-				make
-				-j5
-				testRotateFactor.run
-				true
-				true
-				true
-			
 			
 				make
 				-j2
@@ -2883,11 +2430,28 @@
 			
 			
 				make
+				
 				tests/testGaussianISAM2
 				true
 				false
 				true
 			
+			
+				make
+				-j2
+				check
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				tests/testLieConfig.run
+				true
+				true
+				true
+			
 			
 				make
 				-j2
@@ -2984,6 +2548,401 @@
 				true
 				true
 			
+			
+				make
+				-j3
+				install
+				true
+				false
+				true
+			
+			
+				make
+				-j2
+				clean
+				true
+				true
+				true
+			
+			
+				make
+				-j1
+				check
+				true
+				false
+				true
+			
+			
+				make
+				-j5
+				all
+				true
+				true
+				true
+			
+			
+				cmake
+				..
+				true
+				false
+				true
+			
+			
+				make
+				-j5
+				gtsam-shared
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				gtsam-static
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				timing
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				examples
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				VERBOSE=1 all
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				VERBOSE=1 check
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				check.base
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				timing.base
+				true
+				true
+				true
+			
+			
+				make
+				-j2 VERBOSE=1
+				check.geometry
+				true
+				false
+				true
+			
+			
+				make
+				-j5
+				timing.geometry
+				true
+				true
+				true
+			
+			
+				make
+				-j2 VERBOSE=1
+				check.inference
+				true
+				false
+				true
+			
+			
+				make
+				-j5
+				timing.inference
+				true
+				true
+				true
+			
+			
+				make
+				-j2 VERBOSE=1
+				check.linear
+				true
+				false
+				true
+			
+			
+				make
+				-j5
+				timing.linear
+				true
+				true
+				true
+			
+			
+				make
+				-j2 VERBOSE=1
+				check.nonlinear
+				true
+				false
+				true
+			
+			
+				make
+				-j5
+				timing.nonlinear
+				true
+				true
+				true
+			
+			
+				make
+				-j2 VERBOSE=1
+				check.slam
+				true
+				false
+				true
+			
+			
+				make
+				-j5
+				timing.slam
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				wrap_gtsam
+				true
+				true
+				true
+			
+			
+				make
+				VERBOSE=1
+				wrap_gtsam
+				true
+				false
+				true
+			
+			
+				cpack
+				-G DEB
+				true
+				false
+				true
+			
+			
+				cpack
+				-G RPM
+				true
+				false
+				true
+			
+			
+				cpack
+				-G TGZ
+				true
+				false
+				true
+			
+			
+				cpack
+				--config CPackSourceConfig.cmake
+				true
+				false
+				true
+			
+			
+				make
+				-j5
+				check.discrete
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				check.discrete_unstable
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				check.base_unstable
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				check.dynamics_unstable
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				check.slam_unstable
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				check.unstable
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				wrap_gtsam_build
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				wrap_gtsam_unstable_build
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				wrap_gtsam_unstable
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				wrap
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				wrap_gtsam_distclean
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				wrap_gtsam_unstable_distclean
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				doc
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				doc_clean
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				check
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				check.geometry_unstable
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				check.linear_unstable
+				true
+				true
+				true
+			
+			
+				make
+				-j6 -j8
+				gtsam_unstable-shared
+				true
+				true
+				true
+			
+			
+				make
+				-j6 -j8
+				gtsam_unstable-static
+				true
+				true
+				true
+			
+			
+				make
+				-j6 -j8
+				check.nonlinear_unstable
+				true
+				true
+				true
+			
+			
+				make
+				-j5
+				check.tests
+				true
+				true
+				true
+			
 			
 				make
 				-j5
@@ -3008,6 +2967,45 @@
 				true
 				true
 			
+			
+				make
+				-j2
+				check
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				clean
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				install
+				true
+				true
+				true
+			
+			
+				make
+				-j2
+				all
+				true
+				true
+				true
+			
+			
+				cmake
+				..
+				true
+				false
+				true
+			
 		
 	
 
diff --git a/CMakeLists.txt b/CMakeLists.txt
index d64814829..8ca25e132 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -34,7 +34,7 @@ if(${CMAKE_SOURCE_DIR} STREQUAL ${CMAKE_BINARY_DIR})
   message(FATAL_ERROR "In-source builds not allowed. Please make a new directory (called a build directory) and run CMake from there. You may need to remove CMakeCache.txt. ")
 endif()
 
-# See whether gtsam_unstable is available (it will be present only if we're using an SVN checkout)
+# See whether gtsam_unstable is available (it will be present only if we're using a git checkout)
 if(EXISTS "${PROJECT_SOURCE_DIR}/gtsam_unstable" AND IS_DIRECTORY "${PROJECT_SOURCE_DIR}/gtsam_unstable")
     set(GTSAM_UNSTABLE_AVAILABLE 1)
 else()
@@ -93,7 +93,7 @@ if(MSVC)
 	# Disable autolinking
 	if(NOT Boost_USE_STATIC_LIBS)
 		add_definitions(-DBOOST_ALL_NO_LIB)
-	    add_definitions(-DBOOST_ALL_DYN_LINK)
+	  add_definitions(-DBOOST_ALL_DYN_LINK)
 	endif()
 endif()
 
diff --git a/examples/Pose2SLAMExample_g2o.cpp b/examples/Pose2SLAMExample_g2o.cpp
index 3a4ee9cff..b46a53198 100644
--- a/examples/Pose2SLAMExample_g2o.cpp
+++ b/examples/Pose2SLAMExample_g2o.cpp
@@ -31,7 +31,7 @@ int main(const int argc, const char *argv[]) {
   // Read graph from file
   string g2oFile;
   if (argc < 2)
-    g2oFile = "../../examples/Data/noisyToyGraph.txt";
+    g2oFile = findExampleDataFile("noisyToyGraph.txt");
   else
     g2oFile = argv[1];
 
diff --git a/examples/Pose2SLAMExample_graph.cpp b/examples/Pose2SLAMExample_graph.cpp
index 449144fef..86be74c75 100644
--- a/examples/Pose2SLAMExample_graph.cpp
+++ b/examples/Pose2SLAMExample_graph.cpp
@@ -32,7 +32,8 @@ int main (int argc, char** argv) {
   NonlinearFactorGraph::shared_ptr graph;
   Values::shared_ptr initial;
   SharedDiagonal model = noiseModel::Diagonal::Sigmas((Vector(3) << 0.05, 0.05, 5.0 * M_PI / 180.0));
-  boost::tie(graph, initial) = load2D("../../examples/Data/w100.graph", model);
+  string graph_file = findExampleDataFile("w100.graph");
+  boost::tie(graph, initial) = load2D(graph_file, model);
   initial->print("Initial estimate:\n");
 
   // Add a Gaussian prior on first poses
diff --git a/examples/Pose2SLAMExample_lago.cpp b/examples/Pose2SLAMExample_lago.cpp
index 1f5d80d7b..bb434f3ce 100644
--- a/examples/Pose2SLAMExample_lago.cpp
+++ b/examples/Pose2SLAMExample_lago.cpp
@@ -32,7 +32,7 @@ int main(const int argc, const char *argv[]) {
   // Read graph from file
   string g2oFile;
   if (argc < 2)
-    g2oFile = "../../examples/Data/noisyToyGraph.txt";
+    g2oFile = findExampleDataFile("noisyToyGraph.txt");
   else
     g2oFile = argv[1];
 
diff --git a/examples/Pose2SLAMwSPCG.cpp b/examples/Pose2SLAMwSPCG.cpp
index 2b2f63064..8e8e28570 100644
--- a/examples/Pose2SLAMwSPCG.cpp
+++ b/examples/Pose2SLAMwSPCG.cpp
@@ -104,7 +104,7 @@ int main(int argc, char** argv) {
   LevenbergMarquardtParams parameters;
   parameters.verbosity = NonlinearOptimizerParams::ERROR;
   parameters.verbosityLM = LevenbergMarquardtParams::LAMBDA;
-  parameters.linearSolverType = NonlinearOptimizerParams::CONJUGATE_GRADIENT;
+  parameters.linearSolverType = NonlinearOptimizerParams::Iterative;
 
   {
     parameters.iterativeParams = boost::make_shared();
diff --git a/examples/RangeISAMExample_plaza2.cpp b/examples/RangeISAMExample_plaza2.cpp
index 29b02a271..252372f39 100644
--- a/examples/RangeISAMExample_plaza2.cpp
+++ b/examples/RangeISAMExample_plaza2.cpp
@@ -40,6 +40,7 @@
 #include 
 #include 
 #include 
+#include 
 
 // Standard headers, added last, so we know headers above work on their own
 #include 
@@ -59,9 +60,8 @@ namespace NM = gtsam::noiseModel;
 typedef pair TimedOdometry;
 list readOdometry() {
   list odometryList;
-  ifstream is("../../examples/Data/Plaza2_DR.txt");
-  if (!is)
-    throw runtime_error("../../examples/Data/Plaza2_DR.txt file not found");
+  string data_file = findExampleDataFile("Plaza2_DR.txt");
+  ifstream is(data_file.c_str());
 
   while (is) {
     double t, distance_traveled, delta_heading;
@@ -78,9 +78,8 @@ list readOdometry() {
 typedef boost::tuple RangeTriple;
 vector readTriples() {
   vector triples;
-  ifstream is("../../examples/Data/Plaza2_TD.txt");
-  if (!is)
-    throw runtime_error("../../examples/Data/Plaza2_TD.txt file not found");
+  string data_file = findExampleDataFile("Plaza2_TD.txt");
+  ifstream is(data_file.c_str());
 
   while (is) {
     double t, sender, receiver, range;
diff --git a/examples/StereoVOExample.cpp b/examples/StereoVOExample.cpp
index 8a04d13a3..19725798c 100644
--- a/examples/StereoVOExample.cpp
+++ b/examples/StereoVOExample.cpp
@@ -24,47 +24,48 @@
  */
  
 #include 
-#include 
+#include 
+#include 
+#include 
 #include 
 #include 
-#include 
-#include 
-#include 
 #include 
-#include 
-#include 
-
-
 
 using namespace std;
 using namespace gtsam;
 
 int main(int argc, char** argv){
+
   //create graph object, add first pose at origin with key '1'
   NonlinearFactorGraph graph;
-  Pose3 first_pose = Pose3();
-  graph.push_back(NonlinearEquality(1, first_pose));
+  Pose3 first_pose;
+  graph.push_back(NonlinearEquality(1, Pose3()));
   
   //create factor noise model with 3 sigmas of value 1
   const noiseModel::Isotropic::shared_ptr model = noiseModel::Isotropic::Sigma(3,1);
   //create stereo camera calibration object with .2m between cameras
   const Cal3_S2Stereo::shared_ptr K(new Cal3_S2Stereo(1000, 1000, 0, 320, 240, 0.2));
+
   //create and add stereo factors between first pose (key value 1) and the three landmarks
   graph.push_back(GenericStereoFactor(StereoPoint2(520, 480, 440), model, 1, 3, K));
   graph.push_back(GenericStereoFactor(StereoPoint2(120, 80, 440), model, 1, 4, K));
   graph.push_back(GenericStereoFactor(StereoPoint2(320, 280, 140), model, 1, 5, K));
+
   //create and add stereo factors between second pose and the three landmarks
   graph.push_back(GenericStereoFactor(StereoPoint2(570, 520, 490), model, 2, 3, K));
   graph.push_back(GenericStereoFactor(StereoPoint2(70, 20, 490), model, 2, 4, K));
   graph.push_back(GenericStereoFactor(StereoPoint2(320, 270, 115), model, 2, 5, K));
+
   //create Values object to contain initial estimates of camera poses and landmark locations
-  Values initial_estimate = Values();
+  Values initial_estimate;
+
   //create and add iniital estimates
   initial_estimate.insert(1, first_pose);
   initial_estimate.insert(2, Pose3(Rot3(), Point3(0.1, -0.1, 1.1)));
   initial_estimate.insert(3, Point3(1, 1, 5));
   initial_estimate.insert(4, Point3(-1, 1, 5));
   initial_estimate.insert(5, Point3(0, -0.5, 5));
+
   //create Levenberg-Marquardt optimizer for resulting factor graph, optimize
   LevenbergMarquardtOptimizer optimizer = LevenbergMarquardtOptimizer(graph, initial_estimate);
   Values result = optimizer.optimize();
@@ -72,4 +73,4 @@ int main(int argc, char** argv){
   result.print("Final result:\n");
 
   return 0;
-}
\ No newline at end of file
+}
diff --git a/gtsam.h b/gtsam.h
index b7178d534..678e2a3d6 100644
--- a/gtsam.h
+++ b/gtsam.h
@@ -1481,9 +1481,7 @@ class GaussianISAM {
 
 #include 
 virtual class IterativeOptimizationParameters {
-  string getKernel() const ;
   string getVerbosity() const;
-  void setKernel(string s) ;
   void setVerbosity(string s) ;
   void print() const;
 };
@@ -1516,7 +1514,7 @@ virtual class SubgraphSolverParameters : gtsam::ConjugateGradientParameters {
   void print() const;
 };
 
-class SubgraphSolver  {
+virtual class SubgraphSolver  {
   SubgraphSolver(const gtsam::GaussianFactorGraph &A, const gtsam::SubgraphSolverParameters ¶meters, const gtsam::Ordering& ordering);
   SubgraphSolver(const gtsam::GaussianFactorGraph &Ab1, const gtsam::GaussianFactorGraph &Ab2, const gtsam::SubgraphSolverParameters ¶meters, const gtsam::Ordering& ordering);
   gtsam::VectorValues optimize() const;
@@ -1855,12 +1853,12 @@ virtual class NonlinearOptimizerParams {
   
   void setLinearSolverType(string solver);
   void setOrdering(const gtsam::Ordering& ordering);
-  void setIterativeParams(const gtsam::SubgraphSolverParameters ¶ms);
+  void setIterativeParams(gtsam::IterativeOptimizationParameters* params);
 
   bool isMultifrontal() const;
   bool isSequential() const;
   bool isCholmod() const;
-  bool isCG() const;
+  bool isIterative() const;
 };
 
 bool checkConvergence(double relativeErrorTreshold,
@@ -2199,6 +2197,25 @@ virtual class GeneralSFMFactor2 : gtsam::NoiseModelFactor {
   void serialize() const;
 };
 
+#include 
+template
+virtual class SmartProjectionPoseFactor : gtsam::NonlinearFactor {
+
+  SmartProjectionPoseFactor(double rankTol, double linThreshold,
+      bool manageDegeneracy, bool enableEPI, const POSE& body_P_sensor);
+
+  SmartProjectionPoseFactor(double rankTol);
+  SmartProjectionPoseFactor();
+
+  void add(const gtsam::Point2& measured_i, size_t poseKey_i, const gtsam::noiseModel::Base* noise_i,
+      const CALIBRATION* K_i);
+
+  // enabling serialization functionality
+  //void serialize() const;
+};
+
+typedef gtsam::SmartProjectionPoseFactor SmartProjectionPose3Factor;
+
 
 #include 
 template
@@ -2301,7 +2318,8 @@ virtual class ConstantBias : gtsam::Value {
 #include 
 class ImuFactorPreintegratedMeasurements {
   // Standard Constructor
-  ImuFactorPreintegratedMeasurements(const gtsam::imuBias::ConstantBias& bias, Matrix measuredAccCovariance, Matrix measuredOmegaCovariance, Matrix integrationErrorCovariance);
+  ImuFactorPreintegratedMeasurements(const gtsam::imuBias::ConstantBias& bias, Matrix measuredAccCovariance,Matrix measuredOmegaCovariance, Matrix integrationErrorCovariance, bool use2ndOrderIntegration);
+  ImuFactorPreintegratedMeasurements(const gtsam::imuBias::ConstantBias& bias, Matrix measuredAccCovariance,Matrix measuredOmegaCovariance, Matrix integrationErrorCovariance);
   ImuFactorPreintegratedMeasurements(const gtsam::ImuFactorPreintegratedMeasurements& rhs);
 
   // Testable
@@ -2339,6 +2357,15 @@ class CombinedImuFactorPreintegratedMeasurements {
       Matrix biasAccCovariance,
       Matrix biasOmegaCovariance,
       Matrix biasAccOmegaInit);
+  CombinedImuFactorPreintegratedMeasurements(
+      const gtsam::imuBias::ConstantBias& bias,
+      Matrix measuredAccCovariance,
+      Matrix measuredOmegaCovariance,
+      Matrix integrationErrorCovariance,
+      Matrix biasAccCovariance,
+      Matrix biasOmegaCovariance,
+      Matrix biasAccOmegaInit,
+      bool use2ndOrderIntegration);
   CombinedImuFactorPreintegratedMeasurements(const gtsam::CombinedImuFactorPreintegratedMeasurements& rhs);
 
   // Testable
@@ -2352,8 +2379,7 @@ class CombinedImuFactorPreintegratedMeasurements {
 
 virtual class CombinedImuFactor : gtsam::NonlinearFactor {
   CombinedImuFactor(size_t pose_i, size_t vel_i, size_t pose_j, size_t vel_j, size_t bias_i, size_t bias_j,
-      const gtsam::CombinedImuFactorPreintegratedMeasurements& CombinedPreintegratedMeasurements, Vector gravity, Vector omegaCoriolis,
-      const gtsam::noiseModel::Base* model);
+      const gtsam::CombinedImuFactorPreintegratedMeasurements& CombinedPreintegratedMeasurements, Vector gravity, Vector omegaCoriolis);
 
   // Standard Interface
   gtsam::CombinedImuFactorPreintegratedMeasurements preintegratedMeasurements() const;
diff --git a/gtsam/geometry/Cal3DS2.cpp b/gtsam/geometry/Cal3DS2.cpp
index e6c97cb1a..c75eff022 100644
--- a/gtsam/geometry/Cal3DS2.cpp
+++ b/gtsam/geometry/Cal3DS2.cpp
@@ -25,7 +25,7 @@ namespace gtsam {
 
 /* ************************************************************************* */
 Cal3DS2::Cal3DS2(const Vector &v):
-    fx_(v[0]), fy_(v[1]), s_(v[2]), u0_(v[3]), v0_(v[4]), k1_(v[5]), k2_(v[6]), k3_(v[7]), k4_(v[8]){}
+    fx_(v[0]), fy_(v[1]), s_(v[2]), u0_(v[3]), v0_(v[4]), k1_(v[5]), k2_(v[6]), p1_(v[7]), p2_(v[8]){}
 
 /* ************************************************************************* */
 Matrix Cal3DS2::K() const {
@@ -34,32 +34,64 @@ Matrix Cal3DS2::K() const {
 
 /* ************************************************************************* */
 Vector Cal3DS2::vector() const {
-  return (Vector(9) << fx_, fy_, s_, u0_, v0_, k1_, k2_, k3_, k4_);
+  return (Vector(9) << fx_, fy_, s_, u0_, v0_, k1_, k2_, p1_, p2_);
 }
 
 /* ************************************************************************* */
-void Cal3DS2::print(const std::string& s) const {
-  gtsam::print(K(), s + ".K");
-  gtsam::print(Vector(k()), s + ".k");
+void Cal3DS2::print(const std::string& s_) const {
+  gtsam::print(K(), s_ + ".K");
+  gtsam::print(Vector(k()), s_ + ".k");
 }
 
 /* ************************************************************************* */
 bool Cal3DS2::equals(const Cal3DS2& K, double tol) const {
   if (fabs(fx_ - K.fx_) > tol || fabs(fy_ - K.fy_) > tol || fabs(s_ - K.s_) > tol ||
       fabs(u0_ - K.u0_) > tol || fabs(v0_ - K.v0_) > tol || fabs(k1_ - K.k1_) > tol ||
-      fabs(k2_ - K.k2_) > tol || fabs(k3_ - K.k3_) > tol || fabs(k4_ - K.k4_) > tol)
+      fabs(k2_ - K.k2_) > tol || fabs(p1_ - K.p1_) > tol || fabs(p2_ - K.p2_) > tol)
     return false;
   return true;
 }
 
 /* ************************************************************************* */
-Point2 Cal3DS2::uncalibrate(const Point2& p,
-       boost::optional H1,
-       boost::optional H2) const {
+static Eigen::Matrix D2dcalibration(double x, double y, double xx,
+    double yy, double xy, double rr, double r4, double pnx, double pny,
+    const Eigen::Matrix& DK) {
+  Eigen::Matrix DR1;
+  DR1 << pnx, 0.0, pny, 1.0, 0.0, 0.0, pny, 0.0, 0.0, 1.0;
+  Eigen::Matrix DR2;
+  DR2 << x * rr, x * r4, 2 * xy, rr + 2 * xx, //
+         y * rr, y * r4, rr + 2 * yy, 2 * xy;
+  Eigen::Matrix D;
+  D << DR1, DK * DR2;
+  return D;
+}
 
-  // parameters
-  const double fx = fx_, fy = fy_, s = s_;
-  const double k1 = k1_, k2 = k2_, k3 = k3_, k4 = k4_;
+/* ************************************************************************* */
+static Eigen::Matrix D2dintrinsic(double x, double y, double rr,
+    double g, double k1, double k2, double p1, double p2,
+    const Eigen::Matrix& DK) {
+  const double drdx = 2. * x;
+  const double drdy = 2. * y;
+  const double dgdx = k1 * drdx + k2 * 2. * rr * drdx;
+  const double dgdy = k1 * drdy + k2 * 2. * rr * drdy;
+
+  // Dx = 2*p1*xy + p2*(rr+2*xx);
+  // Dy = 2*p2*xy + p1*(rr+2*yy);
+  const double dDxdx = 2. * p1 * y + p2 * (drdx + 4. * x);
+  const double dDxdy = 2. * p1 * x + p2 * drdy;
+  const double dDydx = 2. * p2 * y + p1 * drdx;
+  const double dDydy = 2. * p2 * x + p1 * (drdy + 4. * y);
+
+  Eigen::Matrix DR;
+  DR << g + x * dgdx + dDxdx, x * dgdy + dDxdy, //
+        y * dgdx + dDydx, g + y * dgdy + dDydy;
+
+  return DK * DR;
+}
+
+/* ************************************************************************* */
+Point2 Cal3DS2::uncalibrate(const Point2& p, boost::optional H1,
+    boost::optional H2) const {
 
   //  rr = x^2 + y^2;
   //  g = (1 + k(1)*rr + k(2)*rr^2);
@@ -68,40 +100,29 @@ Point2 Cal3DS2::uncalibrate(const Point2& p,
   const double x = p.x(), y = p.y(), xy = x * y, xx = x * x, yy = y * y;
   const double rr = xx + yy;
   const double r4 = rr * rr;
-  const double g = 1. + k1 * rr + k2 * r4;
-  const double dx = 2. * k3 * xy + k4 * (rr + 2. * xx);
-  const double dy = 2. * k4 * xy + k3 * (rr + 2. * yy);
+  const double g = 1. + k1_ * rr + k2_ * r4; // scaling factor
 
-  const double pnx = g*x + dx;
-  const double pny = g*y + dy;
+  // tangential component
+  const double dx = 2. * p1_ * xy + p2_ * (rr + 2. * xx);
+  const double dy = 2. * p2_ * xy + p1_ * (rr + 2. * yy);
 
-  // Inlined derivative for calibration
-  if (H1) {
-    *H1 = (Matrix(2, 9) <<  pnx, 0.0, pny, 1.0, 0.0, fx * x * rr + s * y * rr,
-        fx * x * r4 + s * y * r4, fx * 2. * xy + s * (rr + 2. * yy),
-        fx * (rr + 2. * xx) + s * (2. * xy), 0.0, pny, 0.0, 0.0, 1.0,
-        fy * y * rr, fy * y * r4, fy * (rr + 2. * yy), fy * (2. * xy));
-  }
-  // Inlined derivative for points
-  if (H2) {
-    const double dr_dx = 2. * x;
-    const double dr_dy = 2. * y;
-    const double dg_dx = k1 * dr_dx + k2 * 2. * rr * dr_dx;
-    const double dg_dy = k1 * dr_dy + k2 * 2. * rr * dr_dy;
+  // Radial and tangential distortion applied
+  const double pnx = g * x + dx;
+  const double pny = g * y + dy;
 
-    const double dDx_dx = 2. * k3 * y + k4 * (dr_dx + 4. * x);
-    const double dDx_dy = 2. * k3 * x + k4 * dr_dy;
-    const double dDy_dx = 2. * k4 * y + k3 * dr_dx;
-    const double dDy_dy = 2. * k4 * x + k3 * (dr_dy + 4. * y);
+  Eigen::Matrix DK;
+  if (H1 || H2) DK << fx_, s_, 0.0, fy_;
 
-    Matrix DK = (Matrix(2, 2) << fx, s_, 0.0, fy);
-    Matrix DR = (Matrix(2, 2) << g + x * dg_dx + dDx_dx, x * dg_dy + dDx_dy,
-        y * dg_dx + dDy_dx, g + y * dg_dy + dDy_dy);
+  // Derivative for calibration
+  if (H1)
+    *H1 = D2dcalibration(x, y, xx, yy, xy, rr, r4, pnx, pny, DK);
 
-    *H2 = DK * DR;
-  }
+  // Derivative for points
+  if (H2)
+    *H2 = D2dintrinsic(x, y, rr, g, k1_, k2_, p1_, p2_, DK);
 
-  return Point2(fx * pnx + s * pny + u0_, fy * pny + v0_);
+  // Regular uncalibrate after distortion
+  return Point2(fx_ * pnx + s_ * pny + u0_, fy_ * pny + v0_);
 }
 
 /* ************************************************************************* */
@@ -118,14 +139,14 @@ Point2 Cal3DS2::calibrate(const Point2& pi, const double tol) const {
   // iterate until the uncalibrate is close to the actual pixel coordinate
   const int maxIterations = 10;
   int iteration;
-  for ( iteration = 0; iteration < maxIterations; ++iteration ) {
-    if ( uncalibrate(pn).distance(pi) <= tol ) break;
-    const double x = pn.x(), y = pn.y(), xy = x*y, xx = x*x, yy = y*y;
+  for (iteration = 0; iteration < maxIterations; ++iteration) {
+    if (uncalibrate(pn).distance(pi) <= tol) break;
+    const double x = pn.x(), y = pn.y(), xy = x * y, xx = x * x, yy = y * y;
     const double rr = xx + yy;
-    const double g = (1+k1_*rr+k2_*rr*rr);
-    const double dx = 2*k3_*xy + k4_*(rr+2*xx);
-    const double dy = 2*k4_*xy + k3_*(rr+2*yy);
-    pn = (invKPi - Point2(dx,dy))/g;
+    const double g = (1 + k1_ * rr + k2_ * rr * rr);
+    const double dx = 2 * p1_ * xy + p2_ * (rr + 2 * xx);
+    const double dy = 2 * p2_ * xy + p1_ * (rr + 2 * yy);
+    pn = (invKPi - Point2(dx, dy)) / g;
   }
 
   if ( iteration >= maxIterations )
@@ -136,47 +157,28 @@ Point2 Cal3DS2::calibrate(const Point2& pi, const double tol) const {
 
 /* ************************************************************************* */
 Matrix Cal3DS2::D2d_intrinsic(const Point2& p) const {
-  //const double fx = fx_, fy = fy_, s = s_;
-  const double k1 = k1_, k2 = k2_, k3 = k3_, k4 = k4_;
-  //const double x = p.x(), y = p.y(), xx = x*x, yy = y*y, xy = x*y;
-  const double x = p.x(), y = p.y(), xx = x*x, yy = y*y;
+  const double x = p.x(), y = p.y(), xx = x * x, yy = y * y;
   const double rr = xx + yy;
-  const double dr_dx = 2*x;
-  const double dr_dy = 2*y;
-  const double r4 = rr*rr;
-  const double g = 1 + k1*rr + k2*r4;
-  const double dg_dx = k1*dr_dx + k2*2*rr*dr_dx;
-  const double dg_dy = k1*dr_dy + k2*2*rr*dr_dy;
-
-  // Dx = 2*k3*xy + k4*(rr+2*xx);
-  // Dy = 2*k4*xy + k3*(rr+2*yy);
-  const double dDx_dx = 2*k3*y + k4*(dr_dx + 4*x);
-  const double dDx_dy = 2*k3*x + k4*dr_dy;
-  const double dDy_dx = 2*k4*y + k3*dr_dx;
-  const double dDy_dy = 2*k4*x + k3*(dr_dy + 4*y);
-
-  Matrix DK = (Matrix(2, 2) << fx_, s_, 0.0, fy_);
-  Matrix DR = (Matrix(2, 2) << g + x*dg_dx + dDx_dx,     x*dg_dy + dDx_dy,
-             y*dg_dx + dDy_dx, g + y*dg_dy + dDy_dy);
-
-  return DK * DR;
+  const double r4 = rr * rr;
+  const double g = (1 + k1_ * rr + k2_ * r4);
+  Eigen::Matrix DK;
+  DK << fx_, s_, 0.0, fy_;
+  return D2dintrinsic(x, y, rr, g, k1_, k2_, p1_, p2_, DK);
 }
 
 /* ************************************************************************* */
 Matrix Cal3DS2::D2d_calibration(const Point2& p) const {
-  const double x = p.x(), y = p.y(), xx = x*x, yy = y*y, xy = x*y;
+  const double x = p.x(), y = p.y(), xx = x * x, yy = y * y, xy = x * y;
   const double rr = xx + yy;
-  const double r4 = rr*rr;
-  const double fx = fx_, fy = fy_, s = s_;
-  const double g = (1+k1_*rr+k2_*r4);
-  const double dx = 2*k3_*xy + k4_*(rr+2*xx);
-  const double dy = 2*k4_*xy + k3_*(rr+2*yy);
-  const double pnx = g*x + dx;
-  const double pny = g*y + dy;
-
-  return (Matrix(2, 9) <<
-  pnx, 0.0, pny, 1.0, 0.0, fx*x*rr + s*y*rr, fx*x*r4 + s*y*r4, fx*2*xy + s*(rr+2*yy), fx*(rr+2*xx) + s*(2*xy),
-  0.0, pny, 0.0, 0.0, 1.0, fy*y*rr      , fy*y*r4         , fy*(rr+2*yy)         , fy*(2*xy)  );
+  const double r4 = rr * rr;
+  const double g = (1 + k1_ * rr + k2_ * r4);
+  const double dx = 2 * p1_ * xy + p2_ * (rr + 2 * xx);
+  const double dy = 2 * p2_ * xy + p1_ * (rr + 2 * yy);
+  const double pnx = g * x + dx;
+  const double pny = g * y + dy;
+  Eigen::Matrix DK;
+  DK << fx_, s_, 0.0, fy_;
+  return D2dcalibration(x, y, xx, yy, xy, rr, r4, pnx, pny, DK);
 }
 
 /* ************************************************************************* */
diff --git a/gtsam/geometry/Cal3DS2.h b/gtsam/geometry/Cal3DS2.h
index ced05086b..51fe958d6 100644
--- a/gtsam/geometry/Cal3DS2.h
+++ b/gtsam/geometry/Cal3DS2.h
@@ -27,6 +27,15 @@ namespace gtsam {
  * @brief Calibration of a camera with radial distortion
  * @addtogroup geometry
  * \nosubgrouping
+ *
+ * Uses same distortionmodel as OpenCV, with
+ * http://docs.opencv.org/modules/calib3d/doc/camera_calibration_and_3d_reconstruction.html
+ * but using only k1,k2,p1, and p2 coefficients.
+ * K = [ fx s u0 ; 0 fy v0 ; 0 0 1 ]
+ * rr = Pn.x^2 + Pn.y^2
+ * \hat{pn} = (1 + k1*rr + k2*rr^2 ) pn + [ 2*k3 pn.x pn.y + k4 (rr + 2 Pn.x^2) ;
+ *                      k3 (rr + 2 Pn.y^2) + 2*k4 pn.x pn.y  ]
+ * pi = K*pn
  */
 class GTSAM_EXPORT Cal3DS2 : public DerivedValue {
 
@@ -34,28 +43,22 @@ protected:
 
   double fx_, fy_, s_, u0_, v0_ ; // focal length, skew and principal point
   double k1_, k2_ ; // radial 2nd-order and 4th-order
-  double k3_, k4_ ; // tangential distortion
-
-  // K = [ fx s u0 ; 0 fy v0 ; 0 0 1 ]
-  // rr = Pn.x^2 + Pn.y^2
-  // \hat{pn} = (1 + k1*rr + k2*rr^2 ) pn + [ 2*k3 pn.x pn.y + k4 (rr + 2 Pn.x^2) ;
-  //                      k3 (rr + 2 Pn.y^2) + 2*k4 pn.x pn.y  ]
-  // pi = K*pn
+  double p1_, p2_ ; // tangential distortion
 
 public:
   Matrix K() const ;
-  Eigen::Vector4d k() const { return Eigen::Vector4d(k1_, k2_, k3_, k4_); }
+  Eigen::Vector4d k() const { return Eigen::Vector4d(k1_, k2_, p1_, p2_); }
   Vector vector() const ;
 
   /// @name Standard Constructors
   /// @{
 
   /// Default Constructor with only unit focal length
-  Cal3DS2() : fx_(1), fy_(1), s_(0), u0_(0), v0_(0), k1_(0), k2_(0), k3_(0), k4_(0) {}
+  Cal3DS2() : fx_(1), fy_(1), s_(0), u0_(0), v0_(0), k1_(0), k2_(0), p1_(0), p2_(0) {}
 
   Cal3DS2(double fx, double fy, double s, double u0, double v0,
-      double k1, double k2, double k3 = 0.0, double k4 = 0.0) :
-  fx_(fx), fy_(fy), s_(s), u0_(u0), v0_(v0), k1_(k1), k2_(k2), k3_(k3), k4_(k4) {}
+      double k1, double k2, double p1 = 0.0, double p2 = 0.0) :
+  fx_(fx), fy_(fy), s_(s), u0_(u0), v0_(v0), k1_(k1), k2_(k2), p1_(p1), p2_(p2) {}
 
   /// @}
   /// @name Advanced Constructors
@@ -92,18 +95,30 @@ public:
   /// image center in y
   inline double py() const { return v0_;}
 
+  /// First distortion coefficient
+  inline double k1() const { return k1_;}
+
+  /// Second distortion coefficient
+  inline double k2() const { return k2_;}
+
+  /// First tangential distortion coefficient
+  inline double p1() const { return p1_;}
+
+  /// Second tangential distortion coefficient
+  inline double p2() const { return p2_;}
+
   /**
-   * convert intrinsic coordinates xy to image coordinates uv
+   * convert intrinsic coordinates xy to (distorted) image coordinates uv
    * @param p point in intrinsic coordinates
    * @param Dcal optional 2*9 Jacobian wrpt Cal3DS2 parameters
    * @param Dp optional 2*2 Jacobian wrpt intrinsic coordinates
-   * @return point in image coordinates
+   * @return point in (distorted) image coordinates
    */
   Point2 uncalibrate(const Point2& p,
       boost::optional Dcal = boost::none,
       boost::optional Dp = boost::none) const ;
 
-  /// Conver a pixel coordinate to ideal coordinate
+  /// Convert (distorted) image coordinates uv to intrinsic coordinates xy
   Point2 calibrate(const Point2& p, const double tol=1e-5) const;
 
   /// Derivative of uncalibrate wrpt intrinsic coordinates
@@ -148,8 +163,8 @@ private:
     ar & BOOST_SERIALIZATION_NVP(v0_);
     ar & BOOST_SERIALIZATION_NVP(k1_);
     ar & BOOST_SERIALIZATION_NVP(k2_);
-    ar & BOOST_SERIALIZATION_NVP(k3_);
-    ar & BOOST_SERIALIZATION_NVP(k4_);
+    ar & BOOST_SERIALIZATION_NVP(p1_);
+    ar & BOOST_SERIALIZATION_NVP(p2_);
   }
 
 
diff --git a/gtsam/geometry/Cal3Unified.cpp b/gtsam/geometry/Cal3Unified.cpp
index 808cb84a4..e7b408982 100644
--- a/gtsam/geometry/Cal3Unified.cpp
+++ b/gtsam/geometry/Cal3Unified.cpp
@@ -43,7 +43,7 @@ void Cal3Unified::print(const std::string& s) const {
 bool Cal3Unified::equals(const Cal3Unified& K, double tol) const {
   if (fabs(fx_ - K.fx_) > tol || fabs(fy_ - K.fy_) > tol || fabs(s_ - K.s_) > tol ||
       fabs(u0_ - K.u0_) > tol || fabs(v0_ - K.v0_) > tol || fabs(k1_ - K.k1_) > tol ||
-      fabs(k2_ - K.k2_) > tol || fabs(k3_ - K.k3_) > tol || fabs(k4_ - K.k4_) > tol ||
+      fabs(k2_ - K.k2_) > tol || fabs(p1_ - K.p1_) > tol || fabs(p2_ - K.p2_) > tol ||
       fabs(xi_ - K.xi_) > tol)
     return false;
   return true;
diff --git a/gtsam/geometry/Cal3Unified.h b/gtsam/geometry/Cal3Unified.h
index a1d47332a..58e024c27 100644
--- a/gtsam/geometry/Cal3Unified.h
+++ b/gtsam/geometry/Cal3Unified.h
@@ -31,6 +31,14 @@ namespace gtsam {
  * @brief Calibration of a omni-directional camera with mirror + lens radial distortion
  * @addtogroup geometry
  * \nosubgrouping
+ *
+ * Similar to Cal3DS2, does distortion but has additional mirror parameter xi
+ * K = [ fx s u0 ; 0 fy v0 ; 0 0 1 ]
+ * Pn = [ P.x / (1 + xi * \sqrt{P.x^2 + P.y^2 + 1}), P.y / (1 + xi * \sqrt{P.x^2 + P.y^2 + 1})]
+ * rr = Pn.x^2 + Pn.y^2
+ * \hat{pn} = (1 + k1*rr + k2*rr^2 ) pn + [ 2*k3 pn.x pn.y + k4 (rr + 2 Pn.x^2) ;
+ *                      k3 (rr + 2 Pn.y^2) + 2*k4 pn.x pn.y  ]
+ * pi = K*pn
  */
 class GTSAM_EXPORT Cal3Unified : public Cal3DS2 {
 
@@ -41,13 +49,6 @@ private:
 
   double xi_;  // mirror parameter
 
-  // K = [ fx s u0 ; 0 fy v0 ; 0 0 1 ]
-  // Pn = [ P.x / (1 + xi * \sqrt{P.x^2 + P.y^2 + 1}), P.y / (1 + xi * \sqrt{P.x^2 + P.y^2 + 1})]
-  // rr = Pn.x^2 + Pn.y^2
-  // \hat{pn} = (1 + k1*rr + k2*rr^2 ) pn + [ 2*k3 pn.x pn.y + k4 (rr + 2 Pn.x^2) ;
-  //                      k3 (rr + 2 Pn.y^2) + 2*k4 pn.x pn.y  ]
-  // pi = K*pn
-
 public:
   //Matrix K() const ;
   //Eigen::Vector4d k() const { return Base::k(); }
@@ -60,8 +61,8 @@ public:
   Cal3Unified() : Base(), xi_(0) {}
 
   Cal3Unified(double fx, double fy, double s, double u0, double v0,
-      double k1, double k2, double k3 = 0.0, double k4 = 0.0, double xi = 0.0) :
-        Base(fx, fy, s, u0, v0, k1, k2, k3, k4), xi_(xi) {}
+      double k1, double k2, double p1 = 0.0, double p2 = 0.0, double xi = 0.0) :
+        Base(fx, fy, s, u0, v0, k1, k2, p1, p2), xi_(xi) {}
 
   /// @}
   /// @name Advanced Constructors
diff --git a/gtsam/geometry/tests/testCal3DS2.cpp b/gtsam/geometry/tests/testCal3DS2.cpp
index fc457767e..eddd3a674 100644
--- a/gtsam/geometry/tests/testCal3DS2.cpp
+++ b/gtsam/geometry/tests/testCal3DS2.cpp
@@ -60,6 +60,8 @@ TEST( Cal3DS2, Duncalibrate1)
   K.uncalibrate(p, computed, boost::none);
   Matrix numerical = numericalDerivative21(uncalibrate_, K, p, 1e-7);
   CHECK(assert_equal(numerical,computed,1e-5));
+  Matrix separate = K.D2d_calibration(p);
+  CHECK(assert_equal(numerical,separate,1e-5));
 }
 
 /* ************************************************************************* */
@@ -68,6 +70,8 @@ TEST( Cal3DS2, Duncalibrate2)
   Matrix computed; K.uncalibrate(p, boost::none, computed);
   Matrix numerical = numericalDerivative22(uncalibrate_, K, p, 1e-7);
   CHECK(assert_equal(numerical,computed,1e-5));
+  Matrix separate = K.D2d_intrinsic(p);
+  CHECK(assert_equal(numerical,separate,1e-5));
 }
 
 /* ************************************************************************* */
diff --git a/gtsam/inference/Ordering.h b/gtsam/inference/Ordering.h
index ea2ad7eda..7b1a2bb2e 100644
--- a/gtsam/inference/Ordering.h
+++ b/gtsam/inference/Ordering.h
@@ -17,9 +17,11 @@
 
 #pragma once
 
+#include 
 #include 
 #include 
 
+#include 
 #include 
 #include 
 #include 
@@ -135,6 +137,15 @@ namespace gtsam {
     static GTSAM_EXPORT Ordering COLAMDConstrained(const VariableIndex& variableIndex,
       const FastMap& groups);
 
+    /// Return a natural Ordering. Typically used by iterative solvers
+    template 
+    static Ordering Natural(const FactorGraph &fg) {
+      FastSet src = fg.keys();
+      std::vector keys(src.begin(), src.end());
+      std::stable_sort(keys.begin(), keys.end());
+      return Ordering(keys);
+    }
+
     /// @}
 
     /// @name Testable @{
diff --git a/gtsam/linear/ConjugateGradientSolver.cpp b/gtsam/linear/ConjugateGradientSolver.cpp
new file mode 100644
index 000000000..0505f6c01
--- /dev/null
+++ b/gtsam/linear/ConjugateGradientSolver.cpp
@@ -0,0 +1,49 @@
+/*
+ * ConjugateGradientSolver.cpp
+ *
+ *  Created on: Jun 4, 2014
+ *      Author: ydjian
+ */
+
+#include 
+#include 
+#include 
+
+using namespace std;
+
+namespace gtsam {
+
+/*****************************************************************************/
+void ConjugateGradientParameters::print(ostream &os) const {
+   Base::print(os);
+   cout << "ConjugateGradientParameters" << endl
+        << "minIter:       " << minIterations_ << endl
+        << "maxIter:       " << maxIterations_ << endl
+        << "resetIter:     " << reset_ << endl
+        << "eps_rel:       " << epsilon_rel_ << endl
+        << "eps_abs:       " << epsilon_abs_ << endl;
+}
+
+/*****************************************************************************/
+std::string ConjugateGradientParameters::blasTranslator(const BLASKernel value) {
+  std::string s;
+  switch (value) {
+  case ConjugateGradientParameters::GTSAM:      s = "GTSAM" ;      break;
+  default:                                      s = "UNDEFINED" ;  break;
+  }
+  return s;
+}
+
+/*****************************************************************************/
+ConjugateGradientParameters::BLASKernel ConjugateGradientParameters::blasTranslator(const std::string &src) {
+  std::string s = src;  boost::algorithm::to_upper(s);
+  if (s == "GTSAM")  return ConjugateGradientParameters::GTSAM;
+
+  /* default is SBM */
+  return ConjugateGradientParameters::GTSAM;
+}
+
+
+}
+
+
diff --git a/gtsam/linear/ConjugateGradientSolver.h b/gtsam/linear/ConjugateGradientSolver.h
index d1b3b2c7e..6e8509309 100644
--- a/gtsam/linear/ConjugateGradientSolver.h
+++ b/gtsam/linear/ConjugateGradientSolver.h
@@ -12,14 +12,15 @@
 #pragma once
 
 #include 
+#include 
 
 namespace gtsam {
 
 /**
  * parameters for the conjugate gradient method
  */
+class GTSAM_EXPORT ConjugateGradientParameters : public IterativeOptimizationParameters {
 
-class ConjugateGradientParameters : public IterativeOptimizationParameters {
 public:
   typedef IterativeOptimizationParameters Base;
   typedef boost::shared_ptr shared_ptr;
@@ -30,14 +31,23 @@ public:
   double epsilon_rel_;    ///< threshold for relative error decrease
   double epsilon_abs_;    ///< threshold for absolute error decrease
 
-  ConjugateGradientParameters()
-  : minIterations_(1), maxIterations_(500), reset_(501), epsilon_rel_(1e-3), epsilon_abs_(1e-3){}
+  /* Matrix Operation Kernel */
+  enum BLASKernel {
+    GTSAM = 0,        ///< Jacobian Factor Graph of GTSAM
+  } blas_kernel_ ;
 
-  ConjugateGradientParameters(size_t minIterations, size_t maxIterations, size_t reset, double epsilon_rel, double epsilon_abs)
-    : minIterations_(minIterations), maxIterations_(maxIterations), reset_(reset), epsilon_rel_(epsilon_rel), epsilon_abs_(epsilon_abs){}
+  ConjugateGradientParameters()
+    : minIterations_(1), maxIterations_(500), reset_(501), epsilon_rel_(1e-3),
+      epsilon_abs_(1e-3), blas_kernel_(GTSAM) {}
+
+  ConjugateGradientParameters(size_t minIterations, size_t maxIterations, size_t reset,
+    double epsilon_rel, double epsilon_abs, BLASKernel blas)
+    : minIterations_(minIterations), maxIterations_(maxIterations), reset_(reset),
+      epsilon_rel_(epsilon_rel), epsilon_abs_(epsilon_abs), blas_kernel_(blas) {}
 
   ConjugateGradientParameters(const ConjugateGradientParameters &p)
-    : Base(p), minIterations_(p.minIterations_), maxIterations_(p.maxIterations_), reset_(p.reset_), epsilon_rel_(p.epsilon_rel_), epsilon_abs_(p.epsilon_abs_) {}
+    : Base(p), minIterations_(p.minIterations_), maxIterations_(p.maxIterations_), reset_(p.reset_),
+               epsilon_rel_(p.epsilon_rel_), epsilon_abs_(p.epsilon_abs_), blas_kernel_(GTSAM) {}
 
   /* general interface */
   inline size_t minIterations() const { return minIterations_; }
@@ -61,15 +71,81 @@ public:
   inline void setEpsilon_rel(double value) { epsilon_rel_ = value; }
   inline void setEpsilon_abs(double value) { epsilon_abs_ = value; }
 
-  virtual void print() const {
-    Base::print();
-    std::cout << "ConjugateGradientParameters" << std::endl
-              << "minIter:       " << minIterations_ << std::endl
-              << "maxIter:       " << maxIterations_ << std::endl
-              << "resetIter:     " << reset_ << std::endl
-              << "eps_rel:       " << epsilon_rel_ << std::endl
-              << "eps_abs:       " << epsilon_abs_ << std::endl;
-  }
+
+  void print() const { Base::print(); }
+  virtual void print(std::ostream &os) const;
+
+  static std::string blasTranslator(const BLASKernel k) ;
+  static BLASKernel blasTranslator(const std::string &s) ;
 };
 
+/*********************************************************************************************/
+/*
+ * A template of linear preconditioned conjugate gradient method.
+ * System class should support residual(v, g), multiply(v,Av), leftPrecondition(v, S^{-t}v,
+ * rightPrecondition(v, S^{-1}v), scal(alpha,v), dot(v,v), axpy(alpha,x,y)
+ * Note that the residual is in the preconditioned domain. Refer to Section 9.2 of Saad's book.
+ */
+template 
+V preconditionedConjugateGradient(const S &system, const V &initial, const ConjugateGradientParameters ¶meters) {
+
+  V estimate, residual, direction, q1, q2;
+  estimate = residual = direction = q1 = q2 = initial;
+
+  system.residual(estimate, q1);                /* q1 = b-Ax */
+  system.leftPrecondition(q1, residual);        /* r = S^{-T} (b-Ax) */
+  system.rightPrecondition(residual, direction);/* d = S^{-1} r */
+
+  double currentGamma = system.dot(residual, residual), prevGamma, alpha, beta;
+
+  const size_t iMaxIterations = parameters.maxIterations(),
+               iMinIterations = parameters.minIterations(),
+               iReset = parameters.reset() ;
+  const double threshold = std::max(parameters.epsilon_abs(),
+                                    parameters.epsilon() * parameters.epsilon() * currentGamma);
+
+  if (parameters.verbosity() >= ConjugateGradientParameters::COMPLEXITY )
+    std::cout << "[PCG] epsilon = " << parameters.epsilon()
+             << ", max = " << parameters.maxIterations()
+             << ", reset = " << parameters.reset()
+             << ", ||r0||^2 = " << currentGamma
+             << ", threshold = " << threshold << std::endl;
+
+  size_t k;
+  for ( k = 1 ; k <= iMaxIterations && (currentGamma > threshold || k <= iMinIterations) ; k++ ) {
+
+    if ( k % iReset == 0 ) {
+      system.residual(estimate, q1);                      /* q1 = b-Ax */
+      system.leftPrecondition(q1, residual);              /* r = S^{-T} (b-Ax) */
+      system.rightPrecondition(residual, direction);      /* d = S^{-1} r */
+      currentGamma = system.dot(residual, residual);
+    }
+    system.multiply(direction, q1);                       /* q1 = A d */
+    alpha = currentGamma / system.dot(direction, q1);     /* alpha = gamma / (d' A d) */
+    system.axpy(alpha, direction, estimate);              /* estimate += alpha * direction */
+    system.leftPrecondition(q1, q2);                      /* q2 = S^{-T} * q1 */
+    system.axpy(-alpha, q2, residual);                    /* residual -= alpha * q2 */
+    prevGamma = currentGamma;
+    currentGamma = system.dot(residual, residual);        /* gamma = |residual|^2 */
+    beta = currentGamma / prevGamma;
+    system.rightPrecondition(residual, q1);               /* q1 = S^{-1} residual */
+    system.scal(beta, direction);
+    system.axpy(1.0, q1, direction);                      /* direction = q1 + beta * direction */
+
+    if (parameters.verbosity() >= ConjugateGradientParameters::ERROR )
+       std::cout << "[PCG] k = " << k
+                 << ", alpha = " << alpha
+                 << ", beta = " << beta
+                 << ", ||r||^2 = " << currentGamma
+                 << std::endl;
+  }
+  if (parameters.verbosity() >= ConjugateGradientParameters::COMPLEXITY )
+     std::cout << "[PCG] iterations = " << k
+               << ", ||r||^2 = " << currentGamma
+               << std::endl;
+
+  return estimate;
+}
+
+
 }
diff --git a/gtsam/linear/GaussianFactorGraph.cpp b/gtsam/linear/GaussianFactorGraph.cpp
index d6836783b..438a06783 100644
--- a/gtsam/linear/GaussianFactorGraph.cpp
+++ b/gtsam/linear/GaussianFactorGraph.cpp
@@ -54,6 +54,20 @@ namespace gtsam {
     return keys;
   }
 
+  /* ************************************************************************* */
+  std::map GaussianFactorGraph::getKeyDimMap() const {
+    map spec;
+    BOOST_FOREACH ( const GaussianFactor::shared_ptr &gf, *this ) {
+      for ( GaussianFactor::const_iterator it = gf->begin() ; it != gf->end() ; it++ ) {
+        map::iterator it2 = spec.find(*it);
+        if ( it2 == spec.end() ) {
+          spec.insert(make_pair(*it, gf->getDim(it)));
+        }
+      }
+    }
+    return spec;
+  }
+
   /* ************************************************************************* */
     vector GaussianFactorGraph::getkeydim() const {
       // First find dimensions of each variable
diff --git a/gtsam/linear/GaussianFactorGraph.h b/gtsam/linear/GaussianFactorGraph.h
index 692310f26..910b25d1e 100644
--- a/gtsam/linear/GaussianFactorGraph.h
+++ b/gtsam/linear/GaussianFactorGraph.h
@@ -138,6 +138,9 @@ namespace gtsam {
     typedef FastSet Keys;
     Keys keys() const;
 
+    /* return a map of (Key, dimension) */
+    std::map getKeyDimMap() const;
+
     std::vector getkeydim() const;
 
     /** unnormalized error */
diff --git a/gtsam/linear/IterativeSolver.cpp b/gtsam/linear/IterativeSolver.cpp
index 2bba3e12b..ab3ccde13 100644
--- a/gtsam/linear/IterativeSolver.cpp
+++ b/gtsam/linear/IterativeSolver.cpp
@@ -6,25 +6,42 @@
  */
 
 #include 
+#include 
+#include 
 #include 
+#include 
 #include 
 
+using namespace std;
+
 namespace gtsam {
 
-std::string IterativeOptimizationParameters::getKernel() const { return kernelTranslator(kernel_); }
-std::string IterativeOptimizationParameters::getVerbosity() const { return verbosityTranslator(verbosity_); }
-void IterativeOptimizationParameters::setKernel(const std::string &src) { kernel_ = kernelTranslator(src); }
-void IterativeOptimizationParameters::setVerbosity(const std::string &src) { verbosity_ = verbosityTranslator(src); }
+/*****************************************************************************/
+string IterativeOptimizationParameters::getVerbosity() const { return verbosityTranslator(verbosity_); }
 
-IterativeOptimizationParameters::Kernel IterativeOptimizationParameters::kernelTranslator(const std::string &src) {
-  std::string s = src;  boost::algorithm::to_upper(s);
-  if (s == "CG") return IterativeOptimizationParameters::CG;
-  /* default is cg */
-  else return IterativeOptimizationParameters::CG;
+/*****************************************************************************/
+void IterativeOptimizationParameters::setVerbosity(const string &src) { verbosity_ = verbosityTranslator(src); }
+
+/*****************************************************************************/
+void IterativeOptimizationParameters::print() const {
+  print(cout);
 }
 
-IterativeOptimizationParameters::Verbosity IterativeOptimizationParameters::verbosityTranslator(const std::string &src)  {
-  std::string s = src;  boost::algorithm::to_upper(s);
+/*****************************************************************************/
+void IterativeOptimizationParameters::print(ostream &os) const {
+  os << "IterativeOptimizationParameters:" <<  endl
+     << "verbosity:     " << verbosityTranslator(verbosity_) <<  endl;
+}
+
+/*****************************************************************************/
+ ostream& operator<<(ostream &os, const IterativeOptimizationParameters &p) {
+  p.print(os);
+  return os;
+}
+
+/*****************************************************************************/
+IterativeOptimizationParameters::Verbosity IterativeOptimizationParameters::verbosityTranslator(const  string &src)  {
+   string s = src;  boost::algorithm::to_upper(s);
   if (s == "SILENT") return IterativeOptimizationParameters::SILENT;
   else if (s == "COMPLEXITY") return IterativeOptimizationParameters::COMPLEXITY;
   else if (s == "ERROR") return IterativeOptimizationParameters::ERROR;
@@ -32,18 +49,85 @@ IterativeOptimizationParameters::Verbosity IterativeOptimizationParameters::verb
   else return IterativeOptimizationParameters::SILENT;
 }
 
-std::string IterativeOptimizationParameters::kernelTranslator(IterativeOptimizationParameters::Kernel k)  {
-  if ( k == CG ) return "CG";
-  else return "UNKNOWN";
-}
-
-std::string IterativeOptimizationParameters::verbosityTranslator(IterativeOptimizationParameters::Verbosity verbosity)  {
+/*****************************************************************************/
+ string IterativeOptimizationParameters::verbosityTranslator(IterativeOptimizationParameters::Verbosity verbosity)  {
   if (verbosity == SILENT) return "SILENT";
   else if (verbosity == COMPLEXITY) return "COMPLEXITY";
   else if (verbosity == ERROR) return "ERROR";
   else return "UNKNOWN";
 }
 
+/*****************************************************************************/
+VectorValues IterativeSolver::optimize(
+    const GaussianFactorGraph &gfg,
+    boost::optional keyInfo,
+    boost::optional&> lambda)
+{
+  return optimize(
+           gfg,
+           keyInfo ? *keyInfo : KeyInfo(gfg),
+           lambda ? *lambda : std::map());
+}
+
+/*****************************************************************************/
+VectorValues IterativeSolver::optimize (
+  const GaussianFactorGraph &gfg,
+  const KeyInfo &keyInfo,
+  const std::map &lambda)
+{
+  return optimize(gfg, keyInfo, lambda, keyInfo.x0());
+}
+
+/****************************************************************************/
+KeyInfo::KeyInfo(const GaussianFactorGraph &fg, const Ordering &ordering)
+  : ordering_(ordering) {
+  initialize(fg);
+}
+
+/****************************************************************************/
+KeyInfo::KeyInfo(const GaussianFactorGraph &fg)
+  : ordering_(Ordering::Natural(fg)) {
+  initialize(fg);
+}
+
+/****************************************************************************/
+void KeyInfo::initialize(const GaussianFactorGraph &fg){
+  const map colspec = fg.getKeyDimMap();
+  const size_t n = ordering_.size();
+  size_t start = 0;
+
+  for ( size_t i = 0 ; i < n ; ++i ) {
+    const size_t key = ordering_[i];
+    const size_t dim = colspec.find(key)->second;
+    insert(make_pair(key, KeyInfoEntry(i, dim, start)));
+    start += dim ;
+  }
+  numCols_ = start;
+}
+
+/****************************************************************************/
+vector KeyInfo::colSpec() const {
+  std::vector result(size(), 0);
+  BOOST_FOREACH ( const gtsam::KeyInfo::value_type &item, *this ) {
+    result[item.second.index()] = item.second.dim();
+  }
+  return result;
+}
+
+/****************************************************************************/
+VectorValues KeyInfo::x0() const {
+  VectorValues result;
+  BOOST_FOREACH ( const gtsam::KeyInfo::value_type &item, *this ) {
+    result.insert(item.first, Vector::Zero(item.second.dim()));
+  }
+  return result;
+}
+
+/****************************************************************************/
+Vector KeyInfo::x0vector() const {
+  return Vector::Zero(numCols_);
+}
+
 
 }
 
diff --git a/gtsam/linear/IterativeSolver.h b/gtsam/linear/IterativeSolver.h
index 988293a4f..3c397c9e9 100644
--- a/gtsam/linear/IterativeSolver.h
+++ b/gtsam/linear/IterativeSolver.h
@@ -11,17 +11,27 @@
 
 #pragma once
 
+#include 
 #include 
-
+#include 
+#include 
+#include 
+#include 
+#include 
+#include