diff --git a/configure.ac b/configure.ac
index ea68efda2..3bbbff80f 100644
--- a/configure.ac
+++ b/configure.ac
@@ -26,30 +26,44 @@ AC_ARG_ENABLE([debug],
 AM_CONDITIONAL([DEBUG], [test x$debug = xtrue])
 
 # search for gsl
-AM_PATH_GSL(1.1)
+#AM_PATH_GSL(1.1)
      
 # enable using GSL for linalg
-AC_ARG_ENABLE([gsl],
-     [  --enable-gsl    Enable the GSL library],
-     [case "${enableval}" in
-       yes) gsl=true ;;
-       no)  gsl=false ;;
-       *) AC_MSG_ERROR([bad value ${enableval} for --enable-gsl]) ;;
-     esac],[gsl=false])
-
-AM_CONDITIONAL([GSL], [test x$gsl = xtrue])
+#AC_ARG_ENABLE([gsl],
+#     [  --enable-gsl    Enable the GSL library],
+#     [case "${enableval}" in
+#       yes) gsl=true ;;
+#       no)  gsl=false ;;
+#       *) AC_MSG_ERROR([bad value ${enableval} for --enable-gsl]) ;;
+#     esac],[gsl=false])
+#
+#AM_CONDITIONAL([GSL], [test x$gsl = xtrue])
 
      
 # enable using ATLAS for BLAS
-AC_ARG_ENABLE([atlas],
-     [  --enable-atlas    Enable ATLAS optimized BLAS],
-     [case "${enableval}" in
-       yes) atlas=true ;;
-       no)  atlas=false ;;
-       *) AC_MSG_ERROR([bad value ${enableval} for --enable-atlas]) ;;
-     esac],[atlas=false])
+#AC_ARG_ENABLE([atlas],
+#     [  --enable-atlas    Enable ATLAS optimized BLAS],
+#     [case "${enableval}" in
+#       yes) atlas=true ;;
+#       no)  atlas=false ;;
+#       *) AC_MSG_ERROR([bad value ${enableval} for --enable-atlas]) ;;
+#     esac],[atlas=false])
+#
+#AM_CONDITIONAL([ATLAS], [test x$atlas = xtrue])
 
-AM_CONDITIONAL([ATLAS], [test x$atlas = xtrue])
+# search for a blas implementation
+AX_BLAS()
+
+# enable BLAS with general purpose script
+AC_ARG_ENABLE([blas],
+     [  --enable-blas    Enable external BLAS library],
+     [case "${enableval}" in
+       yes) blas=true ;;
+       no)  blas=false ;;
+       *) AC_MSG_ERROR([bad value ${enableval} for --enable-blas]) ;;
+     esac],[blas=false])
+
+AM_CONDITIONAL([USE_BLAS], [test x$blas = xtrue])
 
 # Checks for programs.
 AC_PROG_CXX
diff --git a/cpp/Makefile.am b/cpp/Makefile.am
index a9f239532..94ced7c8f 100644
--- a/cpp/Makefile.am
+++ b/cpp/Makefile.am
@@ -293,27 +293,14 @@ include_HEADERS = $(headers)
 AM_CXXFLAGS += -I.. 
 AM_LDFLAGS = -L../CppUnitLite -lCppUnitLite $(BOOST_LDFLAGS) $(boost_serialization) #-L. -lgtsam
 
-# TESTING: adding cblas implementation with atlas
-#AM_CXXFLAGS += -DCBLAS 
-#libgtsam_la_CPPFLAGS += -DCBLAS
-#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
+# adding cblas implementation with atlas
+if USE_BLAS
+AM_CXXFLAGS += -DCBLAS 
+libgtsam_la_CPPFLAGS += -DCBLAS
+AM_LDFLAGS += $(BLAS_LIBS) $(LIBS) $(FLIBS) # -lcblas -latlas 
+libgtsam_la_LDFLAGS += $(BLAS_LIBS) $(LIBS) $(FLIBS) # -lcblas -latlas
+endif
 
-#######################
-# 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
 
 TESTS = $(check_PROGRAMS)
 CXXLINK = $(LIBTOOL) $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=link \
diff --git a/cpp/Matrix.cpp b/cpp/Matrix.cpp
index 3f60e4cad..c06053067 100644
--- a/cpp/Matrix.cpp
+++ b/cpp/Matrix.cpp
@@ -34,10 +34,6 @@
 #include "Vector.h"
 #include "svdcmp.h"
 
-// use for switching quickly between GSL and nonGSL versions without reconfigure
-//#define REVERTGSL
-
-
 using namespace std;
 namespace ublas = boost::numeric::ublas;
 
@@ -158,7 +154,6 @@ bool assert_equal(const std::list<Matrix>& As, const std::list<Matrix>& Bs, doub
 /* ************************************************************************* */
 void multiplyAdd(double alpha, const Matrix& A, const Vector& x, Vector& e) {
 #if defined CBLAS
-
 	// get sizes
 	const size_t m = A.size1(), n = A.size2();
 
@@ -227,15 +222,15 @@ void transposeMultiplyAdd(double alpha, const Matrix& A, const Vector& e, Vector
 
 	// get pointers
 	const double * Aptr = A.data().begin();
-	const double * Xptr = e.data().begin();
-	double * Eptr = x.data().begin();
+	const double * Eptr = e.data().begin();
+	double * Xptr = x.data().begin();
 
 	// fill in parameters
 	const double beta = 1.0;
 	const int incx = 1, incy = 1, ida = n;
 
 	// execute blas call
-	cblas_dgemv(CblasRowMajor, CblasTrans, m, n, alpha, Aptr, ida, Xptr, incx, beta, Eptr, incy);
+	cblas_dgemv(CblasRowMajor, CblasTrans, m, n, alpha, Aptr, ida, Eptr, incx, beta, Xptr, incy);
 
 #elif defined GSL
 	gsl_vector_const_view eg = gsl_vector_const_view_array(e.data().begin(), e.size());
@@ -482,28 +477,60 @@ inline void householder_update_manual(Matrix &A, int j, double beta, const Vecto
 }
 
 void householder_update(Matrix &A, int j, double beta, const Vector& vjm) {
-	// TODO: SWAP IN ATLAS VERSION OF THE SYSTEM
+	const size_t m = A.size1(), n = A.size2();
+#if defined CBLAS
+
+	// CBLAS version not working, using manual approach
+	householder_update_manual(A,j,beta,vjm);
+
 //	// straight atlas version
-//	const size_t m = A.size1(), n = A.size2(), mj = m-j;
+//	const size_t mj = m-j;
 //
 //	// find pointers to the data
 //	const double * vptr = vjm.data().begin(); // mj long
 //	double * Aptr = A.data().begin() + n*j; // mj x n - note that this starts at row j
 //
 //	// first step: get w = beta*trans(A(j:m,:))*vjm
-//	Vector w(n);
+//	Vector w = zero(n);
 //	double * wptr = w.data().begin();
 //
+//	// DEBUG: create a duplicate version of the problem to solve simultaneously
+//	Matrix aA(A); Vector avjm(vjm);
+//
 //	// execute w generation
 //	cblas_dgemv(CblasRowMajor, CblasTrans, mj, n, beta, Aptr, n, vptr, 1, 0.0, wptr, 1);
 //
+//	// Execute w generation with alternate code
+//	Vector aw(n);
+//	for( size_t c = 0; c < n; c++) {
+//		aw(c) = 0.0;
+//		// dangerous as relies on row-major scheme
+//		const double *a = &aA(j,c), * const v = &avjm(0);
+//		for( size_t r=j, s=0 ; r < m ; r++, s++, a+=n )
+//			// w(c) += A(r,c) * vjm(r-j)
+//			aw(c) += (*a) * v[s];
+//		aw(c) *= beta;
+//	}
+//
+//	print(w, "CBLAS w");
+//	print(aw, "Alternate w");
+//
 //	// second step: rank 1 update A(j:m,:) = v(j:m)*w' + A(j:m,:)
-//	cblas_dger(CblasRowMajor, mj, n, 1.0, vptr, 1, wptr, 1, Aptr, n);
+//	cblas_dger(CblasRowMajor, mj, n, 1.0, vptr, 1, wptr, 1, Aptr, n); // not correct
+//
+//	// Execute second step using alternate code
+//	for( size_t c = 0 ; c < n; c++) {
+//		double wc = aw(c);
+//		double *a = &aA(j,c); const double * const v =&avjm(0);
+//		for( size_t r=j, s=0 ; r < m ; r++, s++, a+=n )
+//			// A(r,c) -= vjm(r-j) * wjn(c-j);
+//			(*a) -= v[s] * wc;
+//	}
+//
+//	// copy in alternate results, which should be correct
+//	A = aA;
 
-
-#ifdef GSL
-#ifndef REVERTGSL
-	const size_t m = A.size1(), n = A.size2();
+#elif defined GSL
 	// use GSL version
 	gsl_vector_const_view v = gsl_vector_const_view_array(vjm.data().begin(), m-j);
 	gsl_matrix_view Ag = gsl_matrix_view_array(A.data().begin(), m, n);
@@ -512,10 +539,6 @@ void householder_update(Matrix &A, int j, double beta, const Vector& vjm) {
 #else
 	householder_update_manual(A,j,beta,vjm);
 #endif
-
-#else
-	householder_update_manual(A,j,beta,vjm);
-#endif
 }
 
 /* ************************************************************************* */
diff --git a/cpp/NonlinearOptimizer-inl.h b/cpp/NonlinearOptimizer-inl.h
index b0e22496f..55ce44c36 100644
--- a/cpp/NonlinearOptimizer-inl.h
+++ b/cpp/NonlinearOptimizer-inl.h
@@ -149,11 +149,12 @@ namespace gtsam {
 
 		// update config
 		shared_config newConfig(new C(expmap(*config_,delta))); // TODO: updateConfig
-		if (verbosity >= TRYCONFIG)
-			newConfig->print("config");
+//		if (verbosity >= TRYCONFIG)
+//			newConfig->print("config");
 
 		// create new optimization state with more adventurous lambda
 		NonlinearOptimizer next(graph_, newConfig, solver_, lambda_ / factor);
+		cout << "next error = " << next.error_ << endl;
 
 		if(lambdaMode >= CAUTIOUS) {
 			throw runtime_error("CAUTIOUS mode not working yet, please use BOUNDED.");
@@ -171,8 +172,9 @@ namespace gtsam {
 
 			// Either we're not cautious, or we are but the adventerous lambda is better than the same one.
 			return next;
-
-		} else {
+		} else if (lambda_ > 1e+10) // if lambda gets too big, something is broken
+			throw runtime_error("Lambda has grown too large!");
+		else {
 
 			// A more adventerous lambda was worse.  If we're cautious, try the same lambda.
 			if(lambdaMode == CAUTIOUS) {
@@ -184,6 +186,7 @@ namespace gtsam {
 			// Either we're not cautious, or the same lambda was worse than the current error.
 			// The more adventerous lambda was worse too, so make lambda more conservative
 			// and keep the same config.
+
 			// TODO: can we avoid copying the config ?
 			if(lambdaMode >= BOUNDED && lambda_ >= 1.0e5) {
 				return NonlinearOptimizer(graph_, newConfig, solver_, lambda_);;
@@ -216,6 +219,7 @@ namespace gtsam {
 			linear->print("linear");
 
 		// try lambda steps with successively larger lambda until we achieve descent
+		cout << "Trying Lambda for the first time" << endl;
 		return try_lambda(*linear, verbosity, lambdaFactor, lambdaMode);
 	}
 
diff --git a/myconfigure b/myconfigure
index 6bec8b101..a6ae991a2 100755
--- a/myconfigure
+++ b/myconfigure
@@ -1 +1 @@
-./configure --prefix=$HOME --with-toolbox=$HOME/toolbox/ --with-boost=/opt/local/include/ --enable-gsl=no --enable-atlas=no CXXFLAGS=" -g -O2 -march=core2 -DNDEBUG" --disable-static
+./configure --prefix=$HOME --with-toolbox=$HOME/toolbox/ --with-boost=/opt/local/include/ --enable-blas --with-blas=atlas CXXFLAGS=" -g -O2 -march=core2 -DNDEBUG" --disable-static