diff --git a/cpp/Factorization.h b/cpp/Factorization.h
new file mode 100644
index 000000000..c290f741c
--- /dev/null
+++ b/cpp/Factorization.h
@@ -0,0 +1,49 @@
+/**
+ * @file    Factorization
+ * @brief   Template Linear solver class that uses a Gaussian Factor Graph
+ * @author  Kai Ni
+ * @author  Frank Dellaert
+ */ 
+
+// $Id: GaussianFactorGraph.h,v 1.24 2009/08/14 20:48:51 acunning Exp $
+
+#pragma once
+
+#include <stdexcept>
+#include "GaussianFactorGraph.h"
+
+namespace gtsam {
+
+	class Ordering;
+
+  /**
+   * A linear system solver using factorization
+   */
+  template <class NonlinearGraph, class Config>
+  class Factorization {
+  private:
+  	boost::shared_ptr<const Ordering> ordering_;
+  	bool useOldEliminate_;
+
+  public:
+  	Factorization(boost::shared_ptr<const Ordering> ordering, bool old=true)
+		: ordering_(ordering), useOldEliminate_(old) {
+  		if (!ordering) throw std::invalid_argument("Factorization constructor: ordering = NULL");
+  	}
+
+  	/**
+  	 * solve for the optimal displacement in the tangent space, and then solve
+  	 * the resulted linear system
+  	 */
+  	VectorConfig optimize(GaussianFactorGraph& fg) const {
+  		return fg.optimize(*ordering_, useOldEliminate_);
+  	}
+
+		/**
+		 * linearize the non-linear graph around the current config
+		 */
+  	boost::shared_ptr<GaussianFactorGraph> linearize(const NonlinearGraph& g, const Config& config) const {
+  		return g.linearize(config);
+  	}
+  };
+}
diff --git a/cpp/GaussianFactorGraph.h b/cpp/GaussianFactorGraph.h
index 18591a9e3..36f19d499 100644
--- a/cpp/GaussianFactorGraph.h
+++ b/cpp/GaussianFactorGraph.h
@@ -255,32 +255,4 @@ namespace gtsam {
 				size_t maxIterations = 0) const;
   };
 
-  /**
-   * A linear system solver using factorization
-   */
-  template <class NonlinearGraph, class Config>
-  class Factorization {
-  private:
-  	boost::shared_ptr<const Ordering> ordering_;
-  	bool useOldEliminate_;
-
-  public:
-  	Factorization(boost::shared_ptr<const Ordering> ordering, bool old=true)
-		: ordering_(ordering), useOldEliminate_(old) {}
-
-  	/**
-  	 * solve for the optimal displacement in the tangent space, and then solve
-  	 * the resulted linear system
-  	 */
-  	VectorConfig optimize(GaussianFactorGraph& fg) const {
-  		return fg.optimize(*ordering_, useOldEliminate_);
-  	}
-
-		/**
-		 * linearize the non-linear graph around the current config
-		 */
-  	boost::shared_ptr<GaussianFactorGraph> linearize(const NonlinearGraph& g, const Config& config) const {
-  		return g.linearize(config);
-  	}
-  };
-}
+} // namespace gtsam
diff --git a/cpp/Makefile.am b/cpp/Makefile.am
index bbc407a69..09eee333b 100644
--- a/cpp/Makefile.am
+++ b/cpp/Makefile.am
@@ -92,7 +92,7 @@ testBinaryBayesNet_SOURCES    = testBinaryBayesNet.cpp
 testBinaryBayesNet_LDADD      = libgtsam.la 
 
 # Gaussian inference 
-headers += GaussianFactorSet.h SharedGaussian.h SharedDiagonal.h VectorConfig.h
+headers += GaussianFactorSet.h SharedGaussian.h SharedDiagonal.h VectorConfig.h Factorization.h
 sources += NoiseModel.cpp Errors.cpp VectorMap.cpp VectorBTree.cpp GaussianFactor.cpp 
 sources += GaussianFactorGraph.cpp GaussianConditional.cpp GaussianBayesNet.cpp
 check_PROGRAMS += testVectorMap testVectorBTree testGaussianFactor testGaussianFactorGraph 
@@ -296,24 +296,24 @@ AM_LDFLAGS = -L../CppUnitLite -lCppUnitLite $(BOOST_LDFLAGS) $(boost_serializati
 # TESTING: adding cblas implementation with atlas
 #AM_CXXFLAGS += -DCBLAS 
 #libgtsam_la_CPPFLAGS += -DCBLAS
-#AM_LDFLAGS += -lcblas -latlas 
-#libgtsam_la_LDFLAGS += -lcblas -latlas
+#AM_LDFLAGS += -L/opt/local/lib -lcblas -latlas -llapack -lptcblas -lptf77blas -lf77blas
+#libgtsam_la_LDFLAGS += -L/opt/local/lib -lcblas -latlas -llapack -lptcblas -lptf77blas -lf77blas
 
 #######################
 # GSL/ATLAS inclusion
 #######################
 #GSL using ATLAS
-if GSL
-AM_CXXFLAGS += -DGSL $(GSL_CFLAGS)
-libgtsam_la_CPPFLAGS += -DGSL $(GSL_CFLAGS)
-if ATLAS
-AM_LDFLAGS += $(GSL_LIBS_NO_CBLAS) -lcblas -latlas 
-libgtsam_la_LDFLAGS += $(GSL_LIBS_NO_CBLAS) -lcblas -latlas
-else 
-AM_LDFLAGS += $(GSL_LIBS) 
-libgtsam_la_LDFLAGS += $(GSL_LIBS) 
-endif
-endif
+#if GSL
+#AM_CXXFLAGS += -DGSL $(GSL_CFLAGS)
+#libgtsam_la_CPPFLAGS += -DGSL $(GSL_CFLAGS)
+#if ATLAS
+#AM_LDFLAGS += $(GSL_LIBS_NO_CBLAS) -lcblas -latlas 
+#libgtsam_la_LDFLAGS += $(GSL_LIBS_NO_CBLAS) -lcblas -latlas
+#else 
+#AM_LDFLAGS += $(GSL_LIBS) 
+#libgtsam_la_LDFLAGS += $(GSL_LIBS) 
+#endif
+#endif
 
 TESTS = $(check_PROGRAMS)
 CXXLINK = $(LIBTOOL) $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=link \
diff --git a/cpp/NonlinearOptimizer-inl.h b/cpp/NonlinearOptimizer-inl.h
index a555da1cd..285b4a127 100644
--- a/cpp/NonlinearOptimizer-inl.h
+++ b/cpp/NonlinearOptimizer-inl.h
@@ -43,14 +43,18 @@ namespace gtsam {
 		return converged;
 	}
 
-	/* ************************************************************************* */
-	// Constructors without the solver
 	/* ************************************************************************* */
 	template<class G, class C, class L, class S, class W>
 	NonlinearOptimizer<G, C, L, S, W>::NonlinearOptimizer(shared_graph graph,
 			shared_config config, shared_solver solver, double lambda) :
-		graph_(graph), config_(config), error_(graph->error(
-				*config)), lambda_(lambda), solver_(solver) {
+		graph_(graph), config_(config), lambda_(lambda), solver_(solver) {
+		if (!graph) throw std::invalid_argument(
+				"NonlinearOptimizer constructor: graph = NULL");
+		if (!config) throw std::invalid_argument(
+				"NonlinearOptimizer constructor: config = NULL");
+		if (!solver) throw std::invalid_argument(
+				"NonlinearOptimizer constructor: solver = NULL");
+		error_ = graph->error(*config);
 	}
 
 	/* ************************************************************************* */
diff --git a/cpp/NonlinearOptimizer.h b/cpp/NonlinearOptimizer.h
index ab6a88b73..cbcdc81c7 100644
--- a/cpp/NonlinearOptimizer.h
+++ b/cpp/NonlinearOptimizer.h
@@ -9,8 +9,9 @@
 #define NONLINEAROPTIMIZER_H_
 
 #include <boost/shared_ptr.hpp>
-#include "NonlinearFactorGraph.h"
 #include "VectorConfig.h"
+#include "NonlinearFactorGraph.h"
+#include "Factorization.h"
 
 namespace gtsam {
 
@@ -69,7 +70,7 @@ namespace gtsam {
 		// keep a configuration and its error
 		// These typically change once per iteration (in a functional way)
 		const shared_config config_;
-		const double error_;
+		double error_; // TODO FD: no more const because in constructor I need to set it after checking :-(
 
 		// keep current lambda for use within LM only
 		// TODO: red flag, should we have an LM class ?