diff --git a/.cproject b/.cproject
index 6664fce3b..05fb247ed 100644
--- a/.cproject
+++ b/.cproject
@@ -322,14 +322,6 @@
true
true
-
- make
- -j2
- testGaussianFactor.run
- true
- true
- true
-
make
-j2
@@ -356,6 +348,7 @@
make
+
tests/testBayesTree.run
true
false
@@ -363,6 +356,7 @@
make
+
testBinaryBayesNet.run
true
false
@@ -410,6 +404,7 @@
make
+
testSymbolicBayesNet.run
true
false
@@ -417,6 +412,7 @@
make
+
tests/testSymbolicFactor.run
true
false
@@ -424,6 +420,7 @@
make
+
testSymbolicFactorGraph.run
true
false
@@ -439,11 +436,20 @@
make
+
tests/testBayesTree
true
false
true
+
+ make
+ -j2
+ testGaussianFactor.run
+ true
+ true
+ true
+
make
-j2
@@ -478,7 +484,6 @@
make
-
testGraph.run
true
false
@@ -574,7 +579,6 @@
make
-
testInference.run
true
false
@@ -582,7 +586,6 @@
make
-
testGaussianBayesNet.run
true
false
@@ -590,7 +593,6 @@
make
-
testGaussianFactor.run
true
false
@@ -598,7 +600,6 @@
make
-
testJunctionTree.run
true
false
@@ -606,7 +607,6 @@
make
-
testSymbolicBayesNet.run
true
false
@@ -614,7 +614,6 @@
make
-
testSymbolicFactorGraph.run
true
false
@@ -692,6 +691,22 @@
true
true
+
+ make
+ -j2
+ check
+ true
+ true
+ true
+
+
+ make
+ -j2
+ check
+ true
+ true
+ true
+
make
-j2
@@ -748,14 +763,6 @@
true
true
-
- make
- -j2
- check
- true
- true
- true
-
make
-j2
@@ -764,6 +771,14 @@
true
true
+
+ make
+ -j2
+ check
+ true
+ true
+ true
+
make
-j2
@@ -820,22 +835,6 @@
true
true
-
- make
- -j2
- testGaussianFactorGraph
- true
- true
- true
-
-
- make
- -j2
- testGaussianJunctionTree
- true
- true
- true
-
make
-j2
@@ -876,6 +875,14 @@
true
true
+
+ make
+ -j2
+ testGaussianFactor.run
+ true
+ true
+ true
+
make
-j2
@@ -974,6 +981,7 @@
make
+
testErrors.run
true
false
@@ -1267,6 +1275,14 @@
true
true
+
+ make
+ -j2
+ check
+ true
+ true
+ true
+
make
-j2
@@ -1349,7 +1365,6 @@
make
-
testSimulated2DOriented.run
true
false
@@ -1389,7 +1404,6 @@
make
-
testSimulated2D.run
true
false
@@ -1397,7 +1411,6 @@
make
-
testSimulated3D.run
true
false
@@ -1427,6 +1440,62 @@
true
true
+
+ make
+ -j2
+ tests/testVector.run
+ true
+ true
+ true
+
+
+ make
+ -j2
+ tests/testMatrix.run
+ true
+ true
+ true
+
+
+ make
+ -j2
+ check
+ true
+ true
+ true
+
+
+ make
+ -j2
+ tests/testNumericalDerivative.run
+ true
+ true
+ true
+
+
+ make
+ -j2
+ tests/testBlockMatrices.run
+ true
+ true
+ true
+
+
+ make
+ -j2
+ clean
+ true
+ true
+ true
+
+
+ make
+ -j2
+ tests/testCholesky.run
+ true
+ true
+ true
+
make
-j2
@@ -1459,6 +1528,38 @@
true
true
+
+ make
+ -j2
+ check
+ true
+ true
+ true
+
+
+ make
+ -j2
+ tests/testGaussianConditional.run
+ true
+ true
+ true
+
+
+ make
+ -j2
+ tests/testGaussianFactorGraph.run
+ true
+ true
+ true
+
+
+ make
+ -j2
+ tests/testGaussianJunctionTree.run
+ true
+ true
+ true
+
make
-j2
@@ -1501,7 +1602,6 @@
make
-
tests/testGaussianISAM2
true
false
@@ -1523,6 +1623,46 @@
true
true
+
+ make
+ -j2
+ install
+ true
+ true
+ true
+
+
+ make
+ -j2
+ clean
+ true
+ true
+ true
+
+
+ make
+ -j2
+ check
+ true
+ true
+ true
+
+
+ make
+ -j2
+ all
+ true
+ true
+ true
+
+
+ make
+ -j2
+ dist
+ true
+ true
+ true
+
make
-j2
@@ -1619,54 +1759,6 @@
true
true
-
- make
- -j2
- install
- true
- true
- true
-
-
- make
- -j2
- clean
- true
- true
- true
-
-
- make
- -j2
- check
- true
- true
- true
-
-
- make
- -j2
- all
- true
- true
- true
-
-
- make
- -j2
- dist
- true
- true
- true
-
-
- make
- -j2
- check
- true
- true
- true
-
make
-j2
@@ -1699,6 +1791,14 @@
true
true
+
+ make
+ -j2
+ check
+ true
+ true
+ true
+
@@ -2021,14 +2121,6 @@
true
true
-
- make
- -j2
- testGaussianFactor.run
- true
- true
- true
-
make
-j2
@@ -2055,6 +2147,7 @@
make
+
tests/testBayesTree.run
true
false
@@ -2062,6 +2155,7 @@
make
+
testBinaryBayesNet.run
true
false
@@ -2109,6 +2203,7 @@
make
+
testSymbolicBayesNet.run
true
false
@@ -2116,6 +2211,7 @@
make
+
tests/testSymbolicFactor.run
true
false
@@ -2123,6 +2219,7 @@
make
+
testSymbolicFactorGraph.run
true
false
@@ -2138,11 +2235,20 @@
make
+
tests/testBayesTree
true
false
true
+
+ make
+ -j2
+ testGaussianFactor.run
+ true
+ true
+ true
+
make
-j2
@@ -2177,7 +2283,6 @@
make
-
testGraph.run
true
false
@@ -2273,7 +2378,6 @@
make
-
testInference.run
true
false
@@ -2281,7 +2385,6 @@
make
-
testGaussianBayesNet.run
true
false
@@ -2289,7 +2392,6 @@
make
-
testGaussianFactor.run
true
false
@@ -2297,7 +2399,6 @@
make
-
testJunctionTree.run
true
false
@@ -2305,7 +2406,6 @@
make
-
testSymbolicBayesNet.run
true
false
@@ -2313,7 +2413,6 @@
make
-
testSymbolicFactorGraph.run
true
false
@@ -2391,6 +2490,22 @@
true
true
+
+ make
+ -j2
+ check
+ true
+ true
+ true
+
+
+ make
+ -j2
+ check
+ true
+ true
+ true
+
make
-j2
@@ -2447,14 +2562,6 @@
true
true
-
- make
- -j2
- check
- true
- true
- true
-
make
-j2
@@ -2463,6 +2570,14 @@
true
true
+
+ make
+ -j2
+ check
+ true
+ true
+ true
+
make
-j2
@@ -2519,22 +2634,6 @@
true
true
-
- make
- -j2
- testGaussianFactorGraph
- true
- true
- true
-
-
- make
- -j2
- testGaussianJunctionTree
- true
- true
- true
-
make
-j2
@@ -2575,6 +2674,14 @@
true
true
+
+ make
+ -j2
+ testGaussianFactor.run
+ true
+ true
+ true
+
make
-j2
@@ -2673,6 +2780,7 @@
make
+
testErrors.run
true
false
@@ -2966,6 +3074,14 @@
true
true
+
+ make
+ -j2
+ check
+ true
+ true
+ true
+
make
-j2
@@ -3048,7 +3164,6 @@
make
-
testSimulated2DOriented.run
true
false
@@ -3088,7 +3203,6 @@
make
-
testSimulated2D.run
true
false
@@ -3096,7 +3210,6 @@
make
-
testSimulated3D.run
true
false
@@ -3126,6 +3239,62 @@
true
true
+
+ make
+ -j2
+ tests/testVector.run
+ true
+ true
+ true
+
+
+ make
+ -j2
+ tests/testMatrix.run
+ true
+ true
+ true
+
+
+ make
+ -j2
+ check
+ true
+ true
+ true
+
+
+ make
+ -j2
+ tests/testNumericalDerivative.run
+ true
+ true
+ true
+
+
+ make
+ -j2
+ tests/testBlockMatrices.run
+ true
+ true
+ true
+
+
+ make
+ -j2
+ clean
+ true
+ true
+ true
+
+
+ make
+ -j2
+ tests/testCholesky.run
+ true
+ true
+ true
+
make
-j2
@@ -3158,6 +3327,38 @@
true
true
+
+ make
+ -j2
+ check
+ true
+ true
+ true
+
+
+ make
+ -j2
+ tests/testGaussianConditional.run
+ true
+ true
+ true
+
+
+ make
+ -j2
+ tests/testGaussianFactorGraph.run
+ true
+ true
+ true
+
+
+ make
+ -j2
+ tests/testGaussianJunctionTree.run
+ true
+ true
+ true
+
make
-j2
@@ -3200,7 +3401,6 @@
make
-
tests/testGaussianISAM2
true
false
@@ -3222,6 +3422,46 @@
true
true
+
+ make
+ -j2
+ install
+ true
+ true
+ true
+
+
+ make
+ -j2
+ clean
+ true
+ true
+ true
+
+
+ make
+ -j2
+ check
+ true
+ true
+ true
+
+
+ make
+ -j2
+ all
+ true
+ true
+ true
+
+
+ make
+ -j2
+ dist
+ true
+ true
+ true
+
make
-j2
@@ -3318,54 +3558,6 @@
true
true
-
- make
- -j2
- install
- true
- true
- true
-
-
- make
- -j2
- clean
- true
- true
- true
-
-
- make
- -j2
- check
- true
- true
- true
-
-
- make
- -j2
- all
- true
- true
- true
-
-
- make
- -j2
- dist
- true
- true
- true
-
-
- make
- -j2
- check
- true
- true
- true
-
make
-j2
@@ -3398,6 +3590,14 @@
true
true
+
+ make
+ -j2
+ check
+ true
+ true
+ true
+
diff --git a/configure.ac b/configure.ac
index af3b269ae..2bbd60d12 100644
--- a/configure.ac
+++ b/configure.ac
@@ -13,9 +13,38 @@ AC_CONFIG_SRCDIR([gtsam/inference/SymbolicFactorGraph.cpp])
AC_CONFIG_SRCDIR([gtsam/linear/GaussianFactor.cpp])
AC_CONFIG_SRCDIR([gtsam/nonlinear/NonlinearOptimizer.cpp])
AC_CONFIG_SRCDIR([gtsam/slam/pose2SLAM.cpp])
-AC_CONFIG_SRCDIR([tests/testSQP.cpp])
+AC_CONFIG_SRCDIR([tests/testTupleValues.cpp])
AC_CONFIG_SRCDIR([examples/SimpleRotation.cpp])
+# pull in all subdirectories of Eigen - FIXME: find a better way of handling this
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/Cholesky/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/Core/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/Core/arch/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/Core/arch/AltiVec/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/Core/arch/Default/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/Core/arch/NEON/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/Core/arch/SSE/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/Core/products/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/Core/util/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/Eigen2Support/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/Eigen2Support/Geometry/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/Eigenvalues/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/Geometry/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/Geometry/arch/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/Householder/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/QR/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/SVD/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/Jacobi/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/misc/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/Sparse/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/LU/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/LU/arch/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/plugins/Makefile.am])
+AC_CONFIG_SRCDIR([gtsam/3rdparty/Eigen/src/StlSupport/Makefile.am])
+
# Check for OS
AC_CANONICAL_HOST # needs to be called at some point earlier
AM_CONDITIONAL([DARWIN], [case $host_os in darwin*) true;; *) false;; esac])
@@ -45,46 +74,6 @@ case $host_os in
;;
esac
-# 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=true])
-ak
-AM_CONDITIONAL([USE_BLAS], test x$blas = xtrue)
-AM_CONDITIONAL([USE_BLAS_MACOS], [test x$blas = xtrue && test x$ISMAC = xtrue])
-AM_CONDITIONAL([USE_BLAS_LINUX], [test x$blas = xtrue && test x$ISMAC = xfalse])
-
-# enable LAPACK
-AC_ARG_ENABLE([lapack],
- [ --enable-lapack Enable external LAPACK library],
- [case "${enableval}" in
- yes) lapack=true ;;
- no) lapack=false ;;
- *) AC_MSG_ERROR([bad value ${enableval} for --enable-lapack]) ;;
- esac],[lapack=true])
-
-AM_CONDITIONAL([USE_LAPACK], test x$lapack = xtrue)
-AM_CONDITIONAL([USE_LAPACK_MACOS], [test x$lapack = xtrue && test x$ISMAC = xtrue])
-AM_CONDITIONAL([USE_LAPACK_LINUX], [test x$lapack = xtrue && test x$ISMAC = xfalse])
-
-# On Mac, we use the Accelerate framework for BLAS/LAPACK
-AM_CONDITIONAL([USE_ACCELERATE_MACOS], [(test x$lapack = xtrue || test x$blas = xtrue) && test x$ISMAC = xtrue])
-
-#enable SparseQR for linear solving
-AC_ARG_ENABLE([spqr],
- [ --enable-spqr Enable SparseQR library support],
- [case "${enableval}" in
- yes) spqr=true ;;
- no) spqr=false ;;
- *) AC_MSG_ERROR([bad value ${enableval} for --enable-spqr]) ;;
- esac],[spqr=false])
-
-AM_CONDITIONAL([USE_SPQR], [test x$spqr = xtrue])
-
# enable profiling
AC_ARG_ENABLE([profiling],
[ --enable-profiling Enable profiling],
@@ -167,13 +156,42 @@ AC_ARG_WITH([ccolamd-lib],
[CCOLAMDLib=${HOME}/lib])
AC_SUBST([CCOLAMDLib])
-# For now we require blas, atlas, and lapack
-#AM_COND_IF([test x$ISMAC = xtrue],
-# [LINALG_CPPFLAGS="-I/System/Library/Frameworks/vecLib.framework/Headers ${CCOLAMDInc} -DGT_USE_LAPACK"],
-# [LINALG_CPPFLAGS="${CCOLAMDInc} -DGT_USE_LAPACK"])
-#AM_COND_IF([test x$ISMAC = xtrue],
-# [LINALG_LDFLAGS="-Wl,/System/Library/Frameworks/Accelerate.framework/Accelerate ${CCOLAMDLib}"],
-# [LINALG_LDFLAGS="-lcblas -latlas -llapack ${CCOLAMDLib}"])
-
-AC_CONFIG_FILES([CppUnitLite/Makefile gtsam/base/Makefile gtsam/geometry/Makefile gtsam/inference/Makefile gtsam/linear/Makefile gtsam/nonlinear/Makefile gtsam/slam/Makefile gtsam/Makefile tests/Makefile examples/Makefile examples/vSLAMexample/Makefile Makefile])
+AC_CONFIG_FILES([CppUnitLite/Makefile \
+gtsam/3rdparty/Eigen/src/Cholesky/Makefile \
+gtsam/3rdparty/Eigen/src/Core/Makefile \
+gtsam/3rdparty/Eigen/src/Core/arch/Makefile \
+gtsam/3rdparty/Eigen/src/Core/arch/AltiVec/Makefile \
+gtsam/3rdparty/Eigen/src/Core/arch/Default/Makefile \
+gtsam/3rdparty/Eigen/src/Core/arch/NEON/Makefile \
+gtsam/3rdparty/Eigen/src/Core/arch/SSE/Makefile \
+gtsam/3rdparty/Eigen/src/Core/products/Makefile \
+gtsam/3rdparty/Eigen/src/Core/util/Makefile \
+gtsam/3rdparty/Eigen/src/Eigenvalues/Makefile \
+gtsam/3rdparty/Eigen/src/Geometry/Makefile \
+gtsam/3rdparty/Eigen/src/Geometry/arch/Makefile \
+gtsam/3rdparty/Eigen/src/Householder/Makefile \
+gtsam/3rdparty/Eigen/src/QR/Makefile \
+gtsam/3rdparty/Eigen/src/SVD/Makefile \
+gtsam/3rdparty/Eigen/src/Jacobi/Makefile \
+gtsam/3rdparty/Eigen/src/misc/Makefile \
+gtsam/3rdparty/Eigen/src/Sparse/Makefile \
+gtsam/3rdparty/Eigen/src/LU/Makefile \
+gtsam/3rdparty/Eigen/src/LU/arch/Makefile \
+gtsam/3rdparty/Eigen/src/plugins/Makefile \
+gtsam/3rdparty/Eigen/src/StlSupport/Makefile \
+gtsam/3rdparty/Eigen/src/Eigen2Support/Makefile \
+gtsam/3rdparty/Eigen/src/Eigen2Support/Geometry/Makefile \
+gtsam/3rdparty/Eigen/src/Makefile \
+gtsam/3rdparty/Eigen/Makefile \
+gtsam/3rdparty/Makefile \
+gtsam/base/Makefile \
+gtsam/geometry/Makefile \
+gtsam/inference/Makefile \
+gtsam/linear/Makefile \
+gtsam/nonlinear/Makefile \
+gtsam/slam/Makefile gtsam/Makefile \
+tests/Makefile \
+examples/Makefile \
+examples/vSLAMexample/Makefile \
+Makefile])
AC_OUTPUT
diff --git a/examples/Makefile.am b/examples/Makefile.am
index 20c4a6650..3f7de2b33 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -27,22 +27,6 @@ AM_LDFLAGS = $(BOOST_LDFLAGS)
LDADD = ../gtsam/libgtsam.la
AM_DEFAULT_SOURCE_EXT = .cpp
-if USE_BLAS_LINUX
-AM_LDFLAGS += -lcblas -latlas
-endif
-
-if USE_LAPACK
-AM_CPPFLAGS += -DGT_USE_LAPACK
-endif
-
-if USE_LAPACK_LINUX
-AM_LDFLAGS += -llapack
-endif
-
-if USE_ACCELERATE_MACOS
-AM_LDFLAGS += -Wl,/System/Library/Frameworks/Accelerate.framework/Accelerate
-endif
-
# rule to run an executable
%.run: % $(LDADD)
./$^
diff --git a/examples/vSLAMexample/Makefile.am b/examples/vSLAMexample/Makefile.am
index 39fb706dc..8f951bb85 100644
--- a/examples/vSLAMexample/Makefile.am
+++ b/examples/vSLAMexample/Makefile.am
@@ -33,22 +33,6 @@ AM_LDFLAGS = $(BOOST_LDFLAGS)
LDADD = ../../gtsam/libgtsam.la
AM_DEFAULT_SOURCE_EXT = .cpp
-if USE_BLAS_LINUX
-AM_LDFLAGS += -lcblas -latlas
-endif
-
-if USE_LAPACK
-AM_CPPFLAGS += -DGT_USE_LAPACK
-endif
-
-if USE_LAPACK_LINUX
-AM_LDFLAGS += -llapack
-endif
-
-if USE_ACCELERATE_MACOS
-AM_LDFLAGS += -Wl,/System/Library/Frameworks/Accelerate.framework/Accelerate
-endif
-
# rule to run an executable
%.run: % $(LDADD)
./$^
diff --git a/gtsam/3rdparty/Eigen/Makefile.am b/gtsam/3rdparty/Eigen/Makefile.am
new file mode 100644
index 000000000..1902d0144
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/Makefile.am
@@ -0,0 +1,21 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS = src
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+# add primary headers
+headers = Array Cholesky Core Dense Eigen Eigen2Support Eigenvalues
+headers += Geometry Householder Jacobi LeastSquares LU QR QtAlignedMalloc
+headers += Sparse StdDeque StdList StdVector SVD
+
+Eigendir = $(pkgincludedir)/3rdparty/Eigen
+Eigen_includedir = $(includedir)/gtsam/3rdparty/Eigen
+Eigen_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/CMakeLists.txt b/gtsam/3rdparty/Eigen/src/CMakeLists.txt
deleted file mode 100644
index c326f374d..000000000
--- a/gtsam/3rdparty/Eigen/src/CMakeLists.txt
+++ /dev/null
@@ -1,7 +0,0 @@
-file(GLOB Eigen_src_subdirectories "*")
-escape_string_as_regex(ESCAPED_CMAKE_CURRENT_SOURCE_DIR "${CMAKE_CURRENT_SOURCE_DIR}")
-foreach(f ${Eigen_src_subdirectories})
- if(NOT f MATCHES "\\.txt" AND NOT f MATCHES "${ESCAPED_CMAKE_CURRENT_SOURCE_DIR}/[.].+" )
- add_subdirectory(${f})
- endif()
-endforeach()
diff --git a/gtsam/3rdparty/Eigen/src/Cholesky/CMakeLists.txt b/gtsam/3rdparty/Eigen/src/Cholesky/CMakeLists.txt
deleted file mode 100644
index d01488b41..000000000
--- a/gtsam/3rdparty/Eigen/src/Cholesky/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_Cholesky_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_Cholesky_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Cholesky COMPONENT Devel
- )
diff --git a/gtsam/3rdparty/Eigen/src/Cholesky/Makefile.am b/gtsam/3rdparty/Eigen/src/Cholesky/Makefile.am
new file mode 100644
index 000000000..83d37ea4a
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/Cholesky/Makefile.am
@@ -0,0 +1,21 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+# add headers from src folder for Eigen
+#Cholesky:
+headers = LDLT.h
+headers += LLT.h
+
+Choleskydir = $(pkgincludedir)/3rdparty/Eigen/src/Cholesky
+Cholesky_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/Cholesky
+Cholesky_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/Core/CMakeLists.txt b/gtsam/3rdparty/Eigen/src/Core/CMakeLists.txt
deleted file mode 100644
index 2346fc2bb..000000000
--- a/gtsam/3rdparty/Eigen/src/Core/CMakeLists.txt
+++ /dev/null
@@ -1,10 +0,0 @@
-FILE(GLOB Eigen_Core_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_Core_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Core COMPONENT Devel
- )
-
-ADD_SUBDIRECTORY(products)
-ADD_SUBDIRECTORY(util)
-ADD_SUBDIRECTORY(arch)
diff --git a/gtsam/3rdparty/Eigen/src/Core/Makefile.am b/gtsam/3rdparty/Eigen/src/Core/Makefile.am
new file mode 100644
index 000000000..bed40a038
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/Core/Makefile.am
@@ -0,0 +1,78 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS = arch products util
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+# add headers from src folder for Eigen
+headers =
+
+#Core:
+headers += ArrayBase.h
+headers += Array.h
+headers += ArrayWrapper.h
+headers += Assign.h
+headers += BandMatrix.h
+headers += Block.h
+headers += BooleanRedux.h
+headers += CommaInitializer.h
+headers += CwiseBinaryOp.h
+headers += CwiseNullaryOp.h
+headers += CwiseUnaryOp.h
+headers += CwiseUnaryView.h
+headers += DenseBase.h
+headers += DenseCoeffsBase.h
+headers += DenseStorage.h
+headers += Diagonal.h
+headers += DiagonalMatrix.h
+headers += DiagonalProduct.h
+headers += Dot.h
+headers += EigenBase.h
+headers += Flagged.h
+headers += ForceAlignedAccess.h
+headers += Functors.h
+headers += Fuzzy.h
+headers += GenericPacketMath.h
+headers += GlobalFunctions.h
+headers += IO.h
+headers += MapBase.h
+headers += Map.h
+headers += MathFunctions.h
+headers += MatrixBase.h
+headers += Matrix.h
+headers += NestByValue.h
+headers += NoAlias.h
+headers += NumTraits.h
+headers += PermutationMatrix.h
+headers += PlainObjectBase.h
+headers += ProductBase.h
+headers += Product.h
+headers += Random.h
+headers += Redux.h
+headers += Replicate.h
+headers += ReturnByValue.h
+headers += Reverse.h
+headers += Select.h
+headers += SelfAdjointView.h
+headers += SelfCwiseBinaryOp.h
+headers += SolveTriangular.h
+headers += StableNorm.h
+headers += Stride.h
+headers += Swap.h
+headers += Transpose.h
+headers += Transpositions.h
+headers += TriangularMatrix.h
+headers += VectorBlock.h
+headers += VectorwiseOp.h
+headers += Visitor.h
+
+Coredir = $(pkgincludedir)/3rdparty/Eigen/src/Core
+Core_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/Core
+Core_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/Core/arch/AltiVec/Makefile.am b/gtsam/3rdparty/Eigen/src/Core/arch/AltiVec/Makefile.am
new file mode 100644
index 000000000..2fea6bc92
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/Core/arch/AltiVec/Makefile.am
@@ -0,0 +1,23 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+headers =
+
+
+# Core/arch/AltiVec:
+headers += Complex.h
+headers += PacketMath.h
+
+AltiVecdir = $(pkgincludedir)/3rdparty/Eigen/src/Core/arch/AltiVec
+AltiVec_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/Core/arch/AltiVec
+AltiVec_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/Core/arch/CMakeLists.txt b/gtsam/3rdparty/Eigen/src/Core/arch/CMakeLists.txt
deleted file mode 100644
index 8456dec15..000000000
--- a/gtsam/3rdparty/Eigen/src/Core/arch/CMakeLists.txt
+++ /dev/null
@@ -1,4 +0,0 @@
-ADD_SUBDIRECTORY(SSE)
-ADD_SUBDIRECTORY(AltiVec)
-ADD_SUBDIRECTORY(NEON)
-ADD_SUBDIRECTORY(Default)
diff --git a/gtsam/3rdparty/Eigen/src/Core/arch/Default/Makefile.am b/gtsam/3rdparty/Eigen/src/Core/arch/Default/Makefile.am
new file mode 100644
index 000000000..d4d44be51
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/Core/arch/Default/Makefile.am
@@ -0,0 +1,21 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+headers =
+
+# Core/arch/Default:
+headers += Settings.h
+
+Defaultdir = $(pkgincludedir)/3rdparty/Eigen/src/Core/arch/Default
+Default_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/Core/arch/Default
+Default_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/Core/arch/Makefile.am b/gtsam/3rdparty/Eigen/src/Core/arch/Makefile.am
new file mode 100644
index 000000000..c831a22bd
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/Core/arch/Makefile.am
@@ -0,0 +1,4 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS = AltiVec Default NEON SSE
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/Core/arch/NEON/Makefile.am b/gtsam/3rdparty/Eigen/src/Core/arch/NEON/Makefile.am
new file mode 100644
index 000000000..950e2f98c
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/Core/arch/NEON/Makefile.am
@@ -0,0 +1,22 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+headers =
+
+# Core/arch/NEON:
+headers += Complex.h
+headers += PacketMath.h
+
+NEONdir = $(pkgincludedir)/3rdparty/Eigen/src/Core/arch/NEON
+NEON_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/Core/arch/NEON
+NEON_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/Core/arch/SSE/Makefile.am b/gtsam/3rdparty/Eigen/src/Core/arch/SSE/Makefile.am
new file mode 100644
index 000000000..7210a137d
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/Core/arch/SSE/Makefile.am
@@ -0,0 +1,22 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+headers =
+# Core/arch/SSE:
+headers += Complex.h
+headers += MathFunctions.h
+headers += PacketMath.h
+
+SSEdir = $(pkgincludedir)/3rdparty/Eigen/src/Core/arch/SSE
+SSE_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/Core/arch/SSE
+SSE_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/Core/products/Makefile.am b/gtsam/3rdparty/Eigen/src/Core/products/Makefile.am
new file mode 100644
index 000000000..33bf76927
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/Core/products/Makefile.am
@@ -0,0 +1,34 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+headers =
+
+# Core/products:
+headers += CoeffBasedProduct.h
+headers += GeneralBlockPanelKernel.h
+headers += GeneralMatrixMatrix.h
+headers += GeneralMatrixMatrixTriangular.h
+headers += GeneralMatrixVector.h
+headers += Parallelizer.h
+headers += SelfadjointMatrixMatrix.h
+headers += SelfadjointMatrixVector.h
+headers += SelfadjointProduct.h
+headers += SelfadjointRank2Update.h
+headers += TriangularMatrixMatrix.h
+headers += TriangularMatrixVector.h
+headers += TriangularSolverMatrix.h
+headers += TriangularSolverVector.h
+
+productsdir = $(pkgincludedir)/3rdparty/Eigen/src/Core/products
+products_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/Core/products
+products_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/Core/util/Makefile.am b/gtsam/3rdparty/Eigen/src/Core/util/Makefile.am
new file mode 100644
index 000000000..9fb84ed8b
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/Core/util/Makefile.am
@@ -0,0 +1,30 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+headers =
+
+# Core/util:
+headers += BlasUtil.h
+headers += Constants.h
+headers += DisableStupidWarnings.h
+headers += ForwardDeclarations.h
+headers += Macros.h
+headers += Memory.h
+headers += Meta.h
+headers += ReenableStupidWarnings.h
+headers += StaticAssert.h
+headers += XprHelper.h
+
+utildir = $(pkgincludedir)/3rdparty/Eigen/src/Core/util
+util_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/Core/util
+util_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/Eigen2Support/CMakeLists.txt b/gtsam/3rdparty/Eigen/src/Eigen2Support/CMakeLists.txt
deleted file mode 100644
index 7ae41b3cb..000000000
--- a/gtsam/3rdparty/Eigen/src/Eigen2Support/CMakeLists.txt
+++ /dev/null
@@ -1,8 +0,0 @@
-FILE(GLOB Eigen_Eigen2Support_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_Eigen2Support_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Eigen2Support COMPONENT Devel
- )
-
-ADD_SUBDIRECTORY(Geometry)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/Eigen2Support/Geometry/CMakeLists.txt b/gtsam/3rdparty/Eigen/src/Eigen2Support/Geometry/CMakeLists.txt
deleted file mode 100644
index c347a8f26..000000000
--- a/gtsam/3rdparty/Eigen/src/Eigen2Support/Geometry/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_Eigen2Support_Geometry_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_Eigen2Support_Geometry_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Eigen2Support/Geometry
- )
diff --git a/gtsam/3rdparty/Eigen/src/Eigen2Support/Geometry/Makefile.am b/gtsam/3rdparty/Eigen/src/Eigen2Support/Geometry/Makefile.am
new file mode 100644
index 000000000..2e351175e
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/Eigen2Support/Geometry/Makefile.am
@@ -0,0 +1,31 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+headers =
+
+# Eigen2SupportGeometry:
+headers += AlignedBox.h
+headers += All.h
+headers += AngleAxis.h
+headers += Hyperplane.h
+headers += ParametrizedLine.h
+headers += Quaternion.h
+headers += Rotation2D.h
+headers += RotationBase.h
+headers += Scaling.h
+headers += Transform.h
+headers += Translation.h
+
+Geometrydir = $(pkgincludedir)/3rdparty/Eigen/src/Eigen2Support/Geometry
+Geometry_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/Eigen2Support/Geometry
+Geometry_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/Eigen2Support/Makefile.am b/gtsam/3rdparty/Eigen/src/Eigen2Support/Makefile.am
new file mode 100644
index 000000000..fdef9f4ec
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/Eigen2Support/Makefile.am
@@ -0,0 +1,35 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS = Geometry
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+headers =
+
+# Eigen2Support:
+headers += Block.h
+headers += Cwise.h
+headers += CwiseOperators.h
+headers += Lazy.h
+headers += LeastSquares.h
+headers += LU.h
+headers += Macros.h
+headers += MathFunctions.h
+headers += Memory.h
+headers += Meta.h
+headers += Minor.h
+headers += QR.h
+headers += SVD.h
+headers += TriangularSolver.h
+headers += VectorBlock.h
+
+Eigen2Supportdir = $(pkgincludedir)/3rdparty/Eigen/src/Eigen2Support
+Eigen2Support_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/Eigen2Support
+Eigen2Support_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/Eigenvalues/Makefile.am b/gtsam/3rdparty/Eigen/src/Eigenvalues/Makefile.am
new file mode 100644
index 000000000..2f86a7c5e
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/Eigenvalues/Makefile.am
@@ -0,0 +1,30 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+headers =
+
+# Eigenvalues:
+headers += ComplexEigenSolver.h
+headers += ComplexSchur.h
+headers += EigenSolver.h
+headers += EigenvaluesCommon.h
+headers += GeneralizedSelfAdjointEigenSolver.h
+headers += HessenbergDecomposition.h
+headers += MatrixBaseEigenvalues.h
+headers += RealSchur.h
+headers += SelfAdjointEigenSolver.h
+headers += Tridiagonalization.h
+
+Eigenvaluesdir = $(pkgincludedir)/3rdparty/Eigen/src/Eigenvalues
+Eigenvalues_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/Eigenvalues
+Eigenvalues_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/Geometry/Makefile.am b/gtsam/3rdparty/Eigen/src/Geometry/Makefile.am
new file mode 100644
index 000000000..d1df4a34d
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/Geometry/Makefile.am
@@ -0,0 +1,34 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS = arch
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+headers =
+
+# Geometry:
+headers += AlignedBox.h
+headers += AngleAxis.h
+headers += EulerAngles.h
+headers += Homogeneous.h
+headers += Hyperplane.h
+headers += OrthoMethods.h
+headers += ParametrizedLine.h
+headers += Quaternion.h
+headers += Rotation2D.h
+headers += RotationBase.h
+headers += Scaling.h
+headers += Transform.h
+headers += Translation.h
+headers += Umeyama.h
+
+Geometrydir = $(pkgincludedir)/3rdparty/Eigen/src/Geometry
+Geometry_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/Geometry
+Geometry_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/Geometry/arch/CMakeLists.txt b/gtsam/3rdparty/Eigen/src/Geometry/arch/CMakeLists.txt
deleted file mode 100644
index 1267a79c7..000000000
--- a/gtsam/3rdparty/Eigen/src/Geometry/arch/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_Geometry_arch_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_Geometry_arch_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Geometry/arch COMPONENT Devel
- )
diff --git a/gtsam/3rdparty/Eigen/src/Geometry/arch/Makefile.am b/gtsam/3rdparty/Eigen/src/Geometry/arch/Makefile.am
new file mode 100644
index 000000000..8007c3f32
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/Geometry/arch/Makefile.am
@@ -0,0 +1,22 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+# add primary headers
+headers =
+
+# Geometry/arch:
+headers += Geometry_SSE.h
+
+archdir = $(pkgincludedir)/3rdparty/Eigen/src/Geometry/arch
+arch_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/Geometry/arch
+arch_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/Householder/CMakeLists.txt b/gtsam/3rdparty/Eigen/src/Householder/CMakeLists.txt
deleted file mode 100644
index ce4937db0..000000000
--- a/gtsam/3rdparty/Eigen/src/Householder/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_Householder_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_Householder_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Householder COMPONENT Devel
- )
diff --git a/gtsam/3rdparty/Eigen/src/Householder/Makefile.am b/gtsam/3rdparty/Eigen/src/Householder/Makefile.am
new file mode 100644
index 000000000..6b475d6dc
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/Householder/Makefile.am
@@ -0,0 +1,23 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+headers =
+
+# Householder:
+headers += BlockHouseholder.h
+headers += Householder.h
+headers += HouseholderSequence.h
+
+Householderdir = $(pkgincludedir)/3rdparty/Eigen/src/Householder
+Householder_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/Householder
+Householder_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/Jacobi/CMakeLists.txt b/gtsam/3rdparty/Eigen/src/Jacobi/CMakeLists.txt
deleted file mode 100644
index 490dac626..000000000
--- a/gtsam/3rdparty/Eigen/src/Jacobi/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_Jacobi_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_Jacobi_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Jacobi COMPONENT Devel
- )
diff --git a/gtsam/3rdparty/Eigen/src/Jacobi/Makefile.am b/gtsam/3rdparty/Eigen/src/Jacobi/Makefile.am
new file mode 100644
index 000000000..e5cdd5ea5
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/Jacobi/Makefile.am
@@ -0,0 +1,20 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+# add primary headers
+# Jacobi
+headers = Jacobi.h
+
+Jacobidir = $(pkgincludedir)/3rdparty/Eigen/src/Jacobi
+Jacobi_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/Jacobi
+Jacobi_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/LU/CMakeLists.txt b/gtsam/3rdparty/Eigen/src/LU/CMakeLists.txt
deleted file mode 100644
index e0d8d78c1..000000000
--- a/gtsam/3rdparty/Eigen/src/LU/CMakeLists.txt
+++ /dev/null
@@ -1,8 +0,0 @@
-FILE(GLOB Eigen_LU_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_LU_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/LU COMPONENT Devel
- )
-
-ADD_SUBDIRECTORY(arch)
diff --git a/gtsam/3rdparty/Eigen/src/LU/Makefile.am b/gtsam/3rdparty/Eigen/src/LU/Makefile.am
new file mode 100644
index 000000000..683fe4c65
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/LU/Makefile.am
@@ -0,0 +1,24 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS = arch
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+headers =
+
+# LU:
+headers += Determinant.h
+headers += FullPivLU.h
+headers += Inverse.h
+headers += PartialPivLU.h
+
+LUdir = $(pkgincludedir)/3rdparty/Eigen/src/LU
+LU_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/LU
+LU_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/LU/arch/CMakeLists.txt b/gtsam/3rdparty/Eigen/src/LU/arch/CMakeLists.txt
deleted file mode 100644
index f6b7ed9ec..000000000
--- a/gtsam/3rdparty/Eigen/src/LU/arch/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_LU_arch_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_LU_arch_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/LU/arch COMPONENT Devel
- )
diff --git a/gtsam/3rdparty/Eigen/src/LU/arch/Makefile.am b/gtsam/3rdparty/Eigen/src/LU/arch/Makefile.am
new file mode 100644
index 000000000..d94a1bb74
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/LU/arch/Makefile.am
@@ -0,0 +1,19 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+# LU/arch:
+headers = Inverse_SSE.h
+
+archdir = $(pkgincludedir)/3rdparty/Eigen/src/LU/arch
+arch_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/LU/arch
+arch_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/Makefile.am b/gtsam/3rdparty/Eigen/src/Makefile.am
new file mode 100644
index 000000000..6353a363e
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/Makefile.am
@@ -0,0 +1,9 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS = Cholesky Core Eigenvalues
+SUBDIRS += Geometry
+SUBDIRS += Jacobi misc Sparse
+SUBDIRS += LU plugins StlSupport
+SUBDIRS += Householder QR SVD
+SUBDIRS += Eigen2Support
diff --git a/gtsam/3rdparty/Eigen/src/QR/CMakeLists.txt b/gtsam/3rdparty/Eigen/src/QR/CMakeLists.txt
deleted file mode 100644
index 96f43d7f5..000000000
--- a/gtsam/3rdparty/Eigen/src/QR/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_QR_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_QR_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/QR COMPONENT Devel
- )
diff --git a/gtsam/3rdparty/Eigen/src/QR/Makefile.am b/gtsam/3rdparty/Eigen/src/QR/Makefile.am
new file mode 100644
index 000000000..eebc3de6f
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/QR/Makefile.am
@@ -0,0 +1,23 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+headers =
+
+# QR:
+headers += ColPivHouseholderQR.h
+headers += FullPivHouseholderQR.h
+headers += HouseholderQR.h
+
+QRdir = $(pkgincludedir)/3rdparty/Eigen/src/QR
+QR_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/QR
+QR_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/SVD/CMakeLists.txt b/gtsam/3rdparty/Eigen/src/SVD/CMakeLists.txt
deleted file mode 100644
index 55efc44b1..000000000
--- a/gtsam/3rdparty/Eigen/src/SVD/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_SVD_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_SVD_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SVD COMPONENT Devel
- )
diff --git a/gtsam/3rdparty/Eigen/src/SVD/Makefile.am b/gtsam/3rdparty/Eigen/src/SVD/Makefile.am
new file mode 100644
index 000000000..01f1c5f49
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/SVD/Makefile.am
@@ -0,0 +1,22 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+headers =
+
+ #SVD:
+headers += JacobiSVD.h
+headers += UpperBidiagonalization.h
+
+SVDdir = $(pkgincludedir)/3rdparty/Eigen/src/SVD
+SVD_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/SVD
+SVD_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/Sparse/CMakeLists.txt b/gtsam/3rdparty/Eigen/src/Sparse/CMakeLists.txt
deleted file mode 100644
index aa1468812..000000000
--- a/gtsam/3rdparty/Eigen/src/Sparse/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_Sparse_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_Sparse_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Sparse COMPONENT Devel
- )
diff --git a/gtsam/3rdparty/Eigen/src/Sparse/Makefile.am b/gtsam/3rdparty/Eigen/src/Sparse/Makefile.am
new file mode 100644
index 000000000..bd5f9d9e3
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/Sparse/Makefile.am
@@ -0,0 +1,45 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+headers =
+
+# Sparse:
+headers += AmbiVector.h
+headers += CompressedStorage.h
+headers += CoreIterators.h
+headers += DynamicSparseMatrix.h
+headers += MappedSparseMatrix.h
+headers += SparseAssign.h
+headers += SparseBlock.h
+headers += SparseCwiseBinaryOp.h
+headers += SparseCwiseUnaryOp.h
+headers += SparseDenseProduct.h
+headers += SparseDiagonalProduct.h
+headers += SparseDot.h
+headers += SparseFuzzy.h
+headers += SparseMatrixBase.h
+headers += SparseMatrix.h
+headers += SparseProduct.h
+headers += SparseRedux.h
+headers += SparseSelfAdjointView.h
+headers += SparseSparseProduct.h
+headers += SparseTranspose.h
+headers += SparseTriangularView.h
+headers += SparseUtil.h
+headers += SparseVector.h
+headers += SparseView.h
+headers += TriangularSolver.h
+
+Sparsedir = $(pkgincludedir)/3rdparty/Eigen/src/Sparse
+Sparse_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/Sparse
+Sparse_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/StlSupport/CMakeLists.txt b/gtsam/3rdparty/Eigen/src/StlSupport/CMakeLists.txt
deleted file mode 100644
index 0f094f637..000000000
--- a/gtsam/3rdparty/Eigen/src/StlSupport/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_StlSupport_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_StlSupport_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/StlSupport COMPONENT Devel
- )
diff --git a/gtsam/3rdparty/Eigen/src/StlSupport/Makefile.am b/gtsam/3rdparty/Eigen/src/StlSupport/Makefile.am
new file mode 100644
index 000000000..3711c44c5
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/StlSupport/Makefile.am
@@ -0,0 +1,24 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+headers =
+
+# StlSupport:
+headers += details.h
+headers += StdDeque.h
+headers += StdList.h
+headers += StdVector.h
+
+StlSupportdir = $(pkgincludedir)/3rdparty/Eigen/src/StlSupport
+StlSupport_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/StlSupport
+StlSupport_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/misc/CMakeLists.txt b/gtsam/3rdparty/Eigen/src/misc/CMakeLists.txt
deleted file mode 100644
index a58ffb745..000000000
--- a/gtsam/3rdparty/Eigen/src/misc/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_misc_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_misc_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/misc COMPONENT Devel
- )
diff --git a/gtsam/3rdparty/Eigen/src/misc/Makefile.am b/gtsam/3rdparty/Eigen/src/misc/Makefile.am
new file mode 100644
index 000000000..977ac554b
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/misc/Makefile.am
@@ -0,0 +1,23 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+headers =
+
+# misc:
+headers += Image.h
+headers += Kernel.h
+headers += Solve.h
+
+miscdir = $(pkgincludedir)/3rdparty/Eigen/src/misc
+misc_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/misc
+misc_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Eigen/src/plugins/CMakeLists.txt b/gtsam/3rdparty/Eigen/src/plugins/CMakeLists.txt
deleted file mode 100644
index 1a1d3ffbd..000000000
--- a/gtsam/3rdparty/Eigen/src/plugins/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_plugins_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_plugins_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/plugins COMPONENT Devel
- )
diff --git a/gtsam/3rdparty/Eigen/src/plugins/Makefile.am b/gtsam/3rdparty/Eigen/src/plugins/Makefile.am
new file mode 100644
index 000000000..3dde0d930
--- /dev/null
+++ b/gtsam/3rdparty/Eigen/src/plugins/Makefile.am
@@ -0,0 +1,27 @@
+# Eigen build/install - just need to copy header files into correct place
+
+# All the sub-directories that need to be built
+SUBDIRS =
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+SUBLIBS =
+
+# use nostdinc to turn off -I. and -I.., we do not need them because
+# header files are qualified so they can be included in external projects.
+AUTOMAKE_OPTIONS = nostdinc
+
+headers =
+
+# plugins:
+headers += ArrayCwiseBinaryOps.h
+headers += ArrayCwiseUnaryOps.h
+headers += BlockMethods.h
+headers += CommonCwiseBinaryOps.h
+headers += CommonCwiseUnaryOps.h
+headers += MatrixCwiseBinaryOps.h
+headers += MatrixCwiseUnaryOps.h
+
+pluginsdir = $(pkgincludedir)/3rdparty/Eigen/src/plugins
+plugins_includedir = $(includedir)/gtsam/3rdparty/Eigen/src/plugins
+plugins_HEADERS = $(headers)
\ No newline at end of file
diff --git a/gtsam/3rdparty/Makefile.am b/gtsam/3rdparty/Makefile.am
new file mode 100644
index 000000000..e30e3be6d
--- /dev/null
+++ b/gtsam/3rdparty/Makefile.am
@@ -0,0 +1,8 @@
+# 3rd Party libraries to be built and installed along with gtsam
+
+# All the sub-directories that need to be built
+SUBDIRS = Eigen
+
+# And the corresponding libraries produced
+# TODO add sublibs when something in 3rdparty actually needs compilation
+# SUBLIBS =
diff --git a/gtsam/Makefile.am b/gtsam/Makefile.am
index 6e4cfd769..0308bdc16 100644
--- a/gtsam/Makefile.am
+++ b/gtsam/Makefile.am
@@ -1,6 +1,6 @@
# All the sub-directories that need to be built
-SUBDIRS = base geometry inference linear nonlinear slam
+SUBDIRS = base geometry inference linear nonlinear slam 3rdparty
# And the corresponding libraries produced
SUBLIBS = base/libbase.la geometry/libgeometry.la inference/libinference.la \
@@ -12,7 +12,3 @@ libgtsam_la_SOURCES =
nodist_EXTRA_libgtsam_la_SOURCES = dummy.cxx
libgtsam_la_LIBADD = $(SUBLIBS) -L$(CCOLAMDLib) $(BOOST_LDFLAGS) -lccolamd
libgtsam_la_LDFLAGS = -no-undefined -version-info 0:0:0
-
-if USE_ACCELERATE_MACOS
-libgtsam_la_LDFLAGS += -Wl,/System/Library/Frameworks/Accelerate.framework/Accelerate
-endif
diff --git a/gtsam/base/DenseQRUtil.cpp b/gtsam/base/DenseQRUtil.cpp
index 1a086c3d2..e10167a67 100644
--- a/gtsam/base/DenseQRUtil.cpp
+++ b/gtsam/base/DenseQRUtil.cpp
@@ -19,149 +19,146 @@
#include
#include
-#include
-#include
-#include
#include // for memcpy and memset
using namespace std;
-namespace ublas = boost::numeric::ublas;
+//namespace ublas = boost::numeric::ublas;
-#ifdef GT_USE_LAPACK
-namespace gtsam {
-
- /* ************************************************************************* */
- int* MakeStairs(Matrix& A) {
-
- const int m = A.size1();
- const int n = A.size2();
- int* Stair = new int[n];
-
- // record the starting pointer of each row
- double* a[m];
- list remained;
- a[0] = A.data().begin();
- for(int i=0; i sorted;
- list::iterator itRemained;
- for(j = 0; j < n; ) {
- // remove the non-zero rows in the current column
- for(itRemained = remained.begin(); itRemained!=remained.end(); ) {
- if (*(a[*itRemained]) != 0) {
- sorted.push_back(*itRemained);
- itRemained = remained.erase(itRemained);
- } else {
- a[*itRemained]++;
- itRemained ++;
- }
- }
-
- // record the stair number
- Stair[j++] = m - remained.size();
-
- if(remained.empty()) break;
- }
-
- // all the remained columns have maximum stair
- for (; j::const_iterator itSorted;
- double* ptr1 = A.data().begin();
- double* ptr2 = A_new.data().begin();
- int row = 0, sizeOfRow = n * sizeof(double);
- for(itSorted=sorted.begin(); itSorted!=sorted.end(); itSorted++, row++)
- memcpy(ptr2+offset[row], ptr1+offset[*itSorted], sizeOfRow);
-
- A = A_new;
-
- return Stair;
- }
-
- void printColumnMajor(const double* a, const int m, const int n) {
- for(int i=0; icol");
- // convert from row major to column major
- ublas::matrix Acolwise(A);
- double *a = Acolwise.data().begin();
- toc("householder_denseqr: row->col");
-
- tic("householder_denseqr: denseqr_front");
- int npiv = min(m,n);
- int b = min(min(m,n),32);
- double W[b*(n+b)];
- DenseQR(m, n, npiv, a, Stair, W);
- toc("householder_denseqr: denseqr_front");
-
- tic("householder_denseqr: col->row");
- int k0 = 0;
- int j0;
- int k;
- memset(A.data().begin(), 0, m*n*sizeof(double));
- for(int j=0; jrow");
-
-
- if(allocedStair) delete[] Stair;
-
- toc("householder_denseqr");
- }
-
- void householder_denseqr_colmajor(ublas::matrix& A, int *Stair) {
- int m = A.size1();
- int n = A.size2();
-
- assert(Stair != NULL);
-
- tic(1, "householder_denseqr");
- int npiv = min(m,n);
- int b = min(min(m,n),32);
- double W[b*(n+b)];
- DenseQR(m, n, npiv, A.data().begin(), Stair, W);
- toc(1, "householder_denseqr");
- }
-
-} // namespace gtsam
-#endif
+//#ifdef GT_USE_LAPACK
+//namespace gtsam {
+//
+// /* ************************************************************************* */
+// int* MakeStairs(Matrix& A) {
+//
+// const int m = A.rows();
+// const int n = A.cols();
+// int* Stair = new int[n];
+//
+// // record the starting pointer of each row
+// double* a[m];
+// list remained;
+// a[0] = A.data();
+// for(int i=0; i sorted;
+// list::iterator itRemained;
+// for(j = 0; j < n; ) {
+// // remove the non-zero rows in the current column
+// for(itRemained = remained.begin(); itRemained!=remained.end(); ) {
+// if (*(a[*itRemained]) != 0) {
+// sorted.push_back(*itRemained);
+// itRemained = remained.erase(itRemained);
+// } else {
+// a[*itRemained]++;
+// itRemained ++;
+// }
+// }
+//
+// // record the stair number
+// Stair[j++] = m - remained.size();
+//
+// if(remained.empty()) break;
+// }
+//
+// // all the remained columns have maximum stair
+// for (; j::const_iterator itSorted;
+// double* ptr1 = A.data();
+// double* ptr2 = A_new.data();
+// int row = 0, sizeOfRow = n * sizeof(double);
+// for(itSorted=sorted.begin(); itSorted!=sorted.end(); itSorted++, row++)
+// memcpy(ptr2+offset[row], ptr1+offset[*itSorted], sizeOfRow);
+//
+// A = A_new;
+//
+// return Stair;
+// }
+//
+// void printColumnMajor(const double* a, const int m, const int n) {
+// for(int i=0; icol");
+// // convert from row major to column major
+// MatrixColMajor Acolwise(A);
+// double *a = Acolwise.data();
+// toc("householder_denseqr: row->col");
+//
+// tic("householder_denseqr: denseqr_front");
+// int npiv = min(m,n);
+// int b = min(min(m,n),32);
+// double W[b*(n+b)];
+// DenseQR(m, n, npiv, a, Stair, W);
+// toc("householder_denseqr: denseqr_front");
+//
+// tic("householder_denseqr: col->row");
+// int k0 = 0;
+// int j0;
+// int k;
+// memset(A.data(), 0, m*n*sizeof(double));
+// for(int j=0; jrow");
+//
+//
+// if(allocedStair) delete[] Stair;
+//
+// toc("householder_denseqr");
+// }
+//
+// void householder_denseqr_colmajor(MatrixColMajor& A, int *Stair) {
+// int m = A.rows();
+// int n = A.cols();
+//
+// assert(Stair != NULL);
+//
+// tic(1, "householder_denseqr");
+// int npiv = min(m,n);
+// int b = min(min(m,n),32);
+// double W[b*(n+b)];
+// DenseQR(m, n, npiv, A.data(), Stair, W);
+// toc(1, "householder_denseqr");
+// }
+//
+//} // namespace gtsam
+//#endif
diff --git a/gtsam/base/DenseQRUtil.h b/gtsam/base/DenseQRUtil.h
index 8f29012b7..86d946b12 100644
--- a/gtsam/base/DenseQRUtil.h
+++ b/gtsam/base/DenseQRUtil.h
@@ -21,17 +21,17 @@
#include
-#ifdef GT_USE_LAPACK
-#include
-
-namespace gtsam {
-
- /** make stairs and speed up householder_denseqr. Stair is defined as the row index of where zero entries start in each column */
- int* MakeStairs(Matrix &A);
-
- /** Householder tranformation, zeros below diagonal */
- void householder_denseqr(Matrix &A, int* Stair = NULL);
-
- void householder_denseqr_colmajor(boost::numeric::ublas::matrix& A, int *Stair);
-}
-#endif
+//#ifdef GT_USE_LAPACK
+//#include
+//
+//namespace gtsam {
+//
+// /** make stairs and speed up householder_denseqr. Stair is defined as the row index of where zero entries start in each column */
+// int* MakeStairs(Matrix &A);
+//
+// /** Householder tranformation, zeros below diagonal */
+// void householder_denseqr(Matrix &A, int* Stair = NULL);
+//
+// void householder_denseqr_colmajor(MatrixColMajor& A, int *Stair);
+//}
+//#endif
diff --git a/gtsam/base/FixedVector.h b/gtsam/base/FixedVector.h
index 5bb7f75b1..ab4a36552 100644
--- a/gtsam/base/FixedVector.h
+++ b/gtsam/base/FixedVector.h
@@ -28,10 +28,10 @@ namespace gtsam {
* size checking.
*/
template
-class FixedVector : public boost::numeric::ublas::bounded_vector ,
+class FixedVector : public Eigen::Matrix,
public Testable > {
public:
- typedef boost::numeric::ublas::bounded_vector Base;
+ typedef Eigen::Matrix Base;
/** default constructor */
FixedVector() {}
@@ -44,7 +44,7 @@ public:
/** Initialize with a C-style array */
FixedVector(const double* values) {
- std::copy(values, values+N, this->data().begin());
+ std::copy(values, values+N, this->data());
}
/**
@@ -69,10 +69,11 @@ public:
* @param constant value
*/
inline static FixedVector repeat(double value) {
- FixedVector v;
- for (size_t i=0; i, Testable {
+
+ /** default constructor - should be unnecessary */
+ LieVector() {}
+
+ /** initialize from a normal vector */
+ LieVector(const Vector& v) : Vector(v) {}
+
+ /** wrap a double */
+ LieVector(double d) : Vector(Vector_(1, d)) {}
+
+ /** constructor with size and initial data, row order ! */
+ LieVector(size_t m, const double* const data);
+
+ /** Specify arguments directly, as in Vector_() - always force these to be doubles */
+ LieVector(size_t m, ...);
+
+ /** get the underlying vector */
+ inline Vector vector() const {
+ return static_cast(*this);
+ }
+
+ /** print @param s optional string naming the object */
+ inline void print(const std::string& name="") const {
+ gtsam::print(vector(), name);
+ }
+
+ /** equality up to tolerance */
+ inline bool equals(const LieVector& expected, double tol=1e-5) const {
+ return gtsam::equal(vector(), expected.vector(), tol);
+ }
+
/**
- * LieVector is a wrapper around vector to allow it to be a Lie type
+ * Returns dimensionality of the tangent space
*/
- struct LieVector : public Vector, public Lie, Testable {
+ inline size_t dim() const { return this->size(); }
- /** default constructor - should be unnecessary */
- LieVector() {}
+ /**
+ * Returns Exponential map update of T
+ * Default implementation calls global binary function
+ */
+ inline LieVector expmap(const Vector& v) const { return LieVector(vector() + v); }
- /** initialize from a normal vector */
- LieVector(const Vector& v) : Vector(v) {}
+ /** expmap around identity */
+ static inline LieVector Expmap(const Vector& v) { return LieVector(v); }
- /** wrap a double */
- LieVector(double d) : Vector(Vector_(1, d)) {}
+ /**
+ * Returns Log map
+ * Default Implementation calls global binary function
+ */
+ inline Vector logmap(const LieVector& lp) const {
+ return lp.vector() - vector();
+ }
- /** constructor with size and initial data, row order ! */
- LieVector(size_t m, const double* const data);
+ /** Logmap around identity - just returns with default cast back */
+ static inline Vector Logmap(const LieVector& p) { return p; }
- /** Specify arguments directly, as in Vector_() - always force these to be doubles */
- LieVector(size_t m, ...);
+ /** compose with another object */
+ inline LieVector compose(const LieVector& p) const {
+ return LieVector(vector() + p);
+ }
- /** get the underlying vector */
- inline Vector vector() const {
- return static_cast(*this);
- }
+ /** between operation */
+ inline LieVector between(const LieVector& l2,
+ boost::optional H1=boost::none,
+ boost::optional H2=boost::none) const {
+ if(H1) *H1 = -eye(dim());
+ if(H2) *H2 = eye(l2.dim());
+ return LieVector(l2.vector() - vector());
+ }
- /** print @param s optional string naming the object */
- inline void print(const std::string& name="") const {
- gtsam::print(vector(), name);
- }
-
- /** equality up to tolerance */
- inline bool equals(const LieVector& expected, double tol=1e-5) const {
- return gtsam::equal(vector(), expected.vector(), tol);
- }
-
- /**
- * Returns dimensionality of the tangent space
- */
- inline size_t dim() const { return this->size(); }
-
- /**
- * Returns Exponential map update of T
- * Default implementation calls global binary function
- */
- inline LieVector expmap(const Vector& v) const { return LieVector(vector() + v); }
-
- /** expmap around identity */
- static inline LieVector Expmap(const Vector& v) { return LieVector(v); }
-
- /**
- * Returns Log map
- * Default Implementation calls global binary function
- */
- inline Vector logmap(const LieVector& lp) const {
- return lp.vector() - vector();
- }
-
- /** Logmap around identity - just returns with default cast back */
- static inline Vector Logmap(const LieVector& p) { return p; }
-
- /** compose with another object */
- inline LieVector compose(const LieVector& p) const {
- return LieVector(vector() + p);
- }
-
- /** between operation */
- inline LieVector between(const LieVector& l2,
- boost::optional H1=boost::none,
- boost::optional H2=boost::none) const {
- if(H1) *H1 = -eye(dim());
- if(H2) *H2 = eye(l2.dim());
- return LieVector(l2.vector() - vector());
- }
-
- /** invert the object and yield a new one */
- inline LieVector inverse() const {
- return LieVector(-1.0 * vector());
- }
- };
+ /** invert the object and yield a new one */
+ inline LieVector inverse() const {
+ return LieVector(-1.0 * vector());
+ }
+};
} // \namespace gtsam
diff --git a/gtsam/base/Makefile.am b/gtsam/base/Makefile.am
index b9da4a7f5..84ec41381 100644
--- a/gtsam/base/Makefile.am
+++ b/gtsam/base/Makefile.am
@@ -13,15 +13,15 @@ check_PROGRAMS =
# base Math
-headers += FixedVector.h types.h blockMatrices.h Matrix-inl.h lapack.h
-sources += Vector.cpp Matrix.cpp cholesky.cpp
-check_PROGRAMS += tests/testFixedVector tests/testVector tests/testMatrix tests/testCholesky
+headers += FixedVector.h types.h blockMatrices.h Matrix-inl.h
+sources += Vector.cpp Matrix.cpp
+sources += cholesky.cpp
+check_PROGRAMS += tests/testFixedVector tests/testVector tests/testMatrix tests/testBlockMatrices
+check_PROGRAMS += tests/testCholesky
check_PROGRAMS += tests/testNumericalDerivative
-if USE_LAPACK
-sources += DenseQR.cpp DenseQRUtil.cpp
-check_PROGRAMS += tests/testDenseQRUtil
-endif
+#sources += DenseQR.cpp DenseQRUtil.cpp
+#check_PROGRAMS += tests/testDenseQRUtil
# Testing
headers += Testable.h TestableAssertions.h numericalDerivative.h
@@ -39,7 +39,8 @@ sources += DSFVector.cpp
check_PROGRAMS += tests/testBTree tests/testDSF tests/testDSFVector
# Timing tests
-noinst_PROGRAMS = tests/timeMatrix tests/timeublas tests/timeVirtual tests/timeTest
+noinst_PROGRAMS = tests/timeMatrix tests/timeVirtual tests/timeTest
+#noinst_PROGRAMS += tests/timeublas
#----------------------------------------------------------------------------------------------------
# Create a libtool library that is not installed
@@ -54,18 +55,9 @@ libbase_la_SOURCES = $(sources)
AM_CPPFLAGS =
-# On Mac, we compile using the BLAS/LAPACK headers in the Accelerate framework
-if USE_ACCELERATE_MACOS
-AM_CPPFLAGS += -I/System/Library/Frameworks/vecLib.framework/Headers
-endif
-
AM_CPPFLAGS += $(BOOST_CPPFLAGS) -I$(CCOLAMDInc) -I$(top_srcdir)
AM_LDFLAGS = $(BOOST_LDFLAGS)
-if USE_BLAS
-AM_CPPFLAGS += -DGT_USE_CBLAS
-endif
-
#----------------------------------------------------------------------------------------------------
# rules to build local programs
#----------------------------------------------------------------------------------------------------
@@ -74,22 +66,6 @@ AM_DEFAULT_SOURCE_EXT = .cpp
AM_LDFLAGS += $(boost_serialization) -L$(CCOLAMDLib) -lccolamd
LDADD = libbase.la ../../CppUnitLite/libCppUnitLite.a
-if USE_BLAS_LINUX
-AM_LDFLAGS += -lcblas -latlas
-endif
-
-if USE_LAPACK
-AM_CPPFLAGS += -DGT_USE_LAPACK
-endif
-
-if USE_LAPACK_LINUX
-AM_LDFLAGS += -llapack
-endif
-
-if USE_ACCELERATE_MACOS
-AM_LDFLAGS += -Wl,/System/Library/Frameworks/Accelerate.framework/Accelerate
-endif
-
# rule to run an executable
%.run: % $(LDADD)
./$^
diff --git a/gtsam/base/Matrix-inl.h b/gtsam/base/Matrix-inl.h
index a038f28d1..da59874ca 100644
--- a/gtsam/base/Matrix-inl.h
+++ b/gtsam/base/Matrix-inl.h
@@ -26,84 +26,82 @@ namespace gtsam {
using namespace std;
-/* ************************************************************************* */
-/**
- * backSubstitute U*x=b
- * @param U an upper triangular matrix
- * @param b an RHS vector
- * @return the solution x of U*x=b
- */
-template
-Vector backSubstituteUpper(const boost::numeric::ublas::matrix_expression& _U,
- const boost::numeric::ublas::vector_expression& _b, bool unit=false) {
- const MATRIX& U(*static_cast(&_U));
- const VECTOR& b(*static_cast(&_b));
- size_t n = U.size2();
-#ifndef NDEBUG
- size_t m = U.size1();
- if (m!=n)
- throw invalid_argument("backSubstituteUpper: U must be square");
-#endif
-
- Vector result(n);
- for (size_t i = n; i > 0; i--) {
- double zi = b(i-1);
- for (size_t j = i+1; j <= n; j++)
- zi -= U(i-1,j-1) * result(j-1);
-#ifndef NDEBUG
- if(!unit && fabs(U(i-1,i-1)) <= numeric_limits::epsilon()) {
- stringstream ss;
- ss << "backSubstituteUpper: U is singular,\n";
- print(U, "U: ", ss);
- throw invalid_argument(ss.str());
- }
-#endif
- if (!unit) zi /= U(i-1,i-1);
- result(i-1) = zi;
- }
-
- return result;
-}
-
-/* ************************************************************************* */
-/**
- * backSubstitute x'*U=b'
- * @param U an upper triangular matrix
- * @param b an RHS vector
- * @param unit, set true if unit triangular
- * @return the solution x of x'*U=b'
- * TODO: use boost
- */
-template
-Vector backSubstituteUpper(const boost::numeric::ublas::vector_expression& _b,
- const boost::numeric::ublas::matrix_expression& _U, bool unit=false) {
- const VECTOR& b(*static_cast(&_b));
- const MATRIX& U(*static_cast(&_U));
- size_t n = U.size2();
-#ifndef NDEBUG
- size_t m = U.size1();
- if (m!=n)
- throw invalid_argument("backSubstituteUpper: U must be square");
-#endif
-
- Vector result(n);
- for (size_t i = 1; i <= n; i++) {
- double zi = b(i-1);
- for (size_t j = 1; j < i; j++)
- zi -= U(j-1,i-1) * result(j-1);
-#ifndef NDEBUG
- if(!unit && fabs(U(i-1,i-1)) <= numeric_limits::epsilon()) {
- stringstream ss;
- ss << "backSubstituteUpper: U is singular,\n";
- print(U, "U: ", ss);
- throw invalid_argument(ss.str());
- }
-#endif
- if (!unit) zi /= U(i-1,i-1);
- result(i-1) = zi;
- }
-
- return result;
-}
+///* ************************************************************************* */
+///**
+// * backSubstitute U*x=b
+// * @param U an upper triangular matrix
+// * @param b an RHS vector
+// * @return the solution x of U*x=b
+// */
+//template
+//Vector backSubstituteUpper(const MATRIX& U, const VECTOR& b, bool unit=false) {
+//// const MATRIX& U(*static_cast(&_U));
+//// const VECTOR& b(*static_cast(&_b));
+// size_t n = U.cols();
+//#ifndef NDEBUG
+// size_t m = U.rows();
+// if (m!=n)
+// throw invalid_argument("backSubstituteUpper: U must be square");
+//#endif
+//
+// Vector result(n);
+// for (size_t i = n; i > 0; i--) {
+// double zi = b(i-1);
+// for (size_t j = i+1; j <= n; j++)
+// zi -= U(i-1,j-1) * result(j-1);
+//#ifndef NDEBUG
+// if(!unit && fabs(U(i-1,i-1)) <= numeric_limits::epsilon()) {
+// stringstream ss;
+// ss << "backSubstituteUpper: U is singular,\n";
+// print(U, "U: ", ss);
+// throw invalid_argument(ss.str());
+// }
+//#endif
+// if (!unit) zi /= U(i-1,i-1);
+// result(i-1) = zi;
+// }
+//
+// return result;
+//}
+//
+///* ************************************************************************* */
+///**
+// * backSubstitute x'*U=b'
+// * @param U an upper triangular matrix
+// * @param b an RHS vector
+// * @param unit, set true if unit triangular
+// * @return the solution x of x'*U=b'
+// * TODO: use boost
+// */
+//template
+//Vector backSubstituteUpper(const VECTOR& b, const MATRIX& U, bool unit=false) {
+//// const VECTOR& b(*static_cast(&_b));
+//// const MATRIX& U(*static_cast(&_U));
+// size_t n = U.cols();
+//#ifndef NDEBUG
+// size_t m = U.rows();
+// if (m!=n)
+// throw invalid_argument("backSubstituteUpper: U must be square");
+//#endif
+//
+// Vector result(n);
+// for (size_t i = 1; i <= n; i++) {
+// double zi = b(i-1);
+// for (size_t j = 1; j < i; j++)
+// zi -= U(j-1,i-1) * result(j-1);
+//#ifndef NDEBUG
+// if(!unit && fabs(U(i-1,i-1)) <= numeric_limits::epsilon()) {
+// stringstream ss;
+// ss << "backSubstituteUpper: U is singular,\n";
+// print(U, "U: ", ss);
+// throw invalid_argument(ss.str());
+// }
+//#endif
+// if (!unit) zi /= U(i-1,i-1);
+// result(i-1) = zi;
+// }
+//
+// return result;
+//}
}
diff --git a/gtsam/base/Matrix.cpp b/gtsam/base/Matrix.cpp
index 330338a49..f9a604290 100644
--- a/gtsam/base/Matrix.cpp
+++ b/gtsam/base/Matrix.cpp
@@ -22,143 +22,108 @@
#include
#include
-#ifdef GT_USE_CBLAS
-extern "C" {
-#include
-}
-#endif
-
-#ifdef GT_USE_LAPACK
-extern "C" {
-#include
-}
-#endif
-
-#include
#include
-#include
-#include
-#include
-#include
-#include
-#include
#include
+#include
+#include
+
+#include
#include
#include
using namespace std;
-namespace ublas = boost::numeric::ublas;
namespace gtsam {
/** Explicit instantiations of template functions for standard types */
-template Vector backSubstituteUpper(const boost::numeric::ublas::matrix_expression& U,
- const boost::numeric::ublas::vector_expression& b, bool unit);
-template Vector backSubstituteUpper(const boost::numeric::ublas::vector_expression& b,
- const boost::numeric::ublas::matrix_expression& U, bool unit);
+//template Vector backSubstituteUpper(const Matrix& U, const Vector& b, bool unit);
+//template Vector backSubstituteUpper(const Vector& b, const Matrix& U, bool unit);
/* ************************************************************************* */
Matrix Matrix_( size_t m, size_t n, const double* const data) {
- Matrix A(m,n);
- copy(data, data+m*n, A.data().begin());
- return A;
+ Matrix A(m,n);
+ copy(data, data+m*n, A.data());
+ return A;
}
/* ************************************************************************* */
Matrix Matrix_( size_t m, size_t n, const Vector& v)
{
- Matrix A(m,n);
- // column-wise copy
- for( size_t j = 0, k=0 ; j < n ; j++)
- for( size_t i = 0; i < m ; i++,k++)
- A(i,j) = v(k);
- return A;
+ Matrix A(m,n);
+ // column-wise copy
+ for( size_t j = 0, k=0 ; j < n ; j++)
+ for( size_t i = 0; i < m ; i++,k++)
+ A(i,j) = v(k);
+ return A;
}
/* ************************************************************************* */
Matrix Matrix_(size_t m, size_t n, ...) {
- Matrix A(m,n);
- va_list ap;
- va_start(ap, n);
- for( size_t i = 0 ; i < m ; i++)
- for( size_t j = 0 ; j < n ; j++) {
- double value = va_arg(ap, double);
- A(i,j) = value;
- }
- va_end(ap);
- return A;
+ Matrix A(m,n);
+ va_list ap;
+ va_start(ap, n);
+ for( size_t i = 0 ; i < m ; i++)
+ for( size_t j = 0 ; j < n ; j++) {
+ double value = va_arg(ap, double);
+ A(i,j) = value;
+ }
+ va_end(ap);
+ return A;
}
/* ************************************************************************* */
-/** create a matrix with value zero */
-/* ************************************************************************* */
-Matrix zeros( size_t m, size_t n )
-{
- Matrix A(m,n, 0.0);
- return A;
+Matrix zeros( size_t m, size_t n ) {
+ return Matrix::Zero(m,n);
}
-/**
- * Identity matrix
- */
-Matrix eye( size_t m, size_t n){
- Matrix A = zeros(m,n);
- for(size_t i = 0; i tol)
- equal = false;
- }
-
-
-
- return equal;
+ return v.asDiagonal();
}
/* ************************************************************************* */
bool assert_equal(const Matrix& expected, const Matrix& actual, double tol) {
- if (equal_with_abs_tol(expected,actual,tol)) return true;
+ if (equal_with_abs_tol(expected,actual,tol)) return true;
- size_t n1 = expected.size2(), m1 = expected.size1();
- size_t n2 = actual.size2(), m2 = actual.size1();
+ size_t n1 = expected.cols(), m1 = expected.rows();
+ size_t n2 = actual.cols(), m2 = actual.rows();
- cout << "not equal:" << endl;
- print(expected,"expected = ");
- print(actual,"actual = ");
- if(m1!=m2 || n1!=n2)
- cout << m1 << "," << n1 << " != " << m2 << "," << n2 << endl;
- else
- print(actual-expected, "actual - expected = ");
- return false;
+ cout << "not equal:" << endl;
+ print(expected,"expected = ");
+ print(actual,"actual = ");
+ if(m1!=m2 || n1!=n2)
+ cout << m1 << "," << n1 << " != " << m2 << "," << n2 << endl;
+ else {
+ Matrix diff = actual-expected;
+ print(diff, "actual - expected = ");
+ }
+ return false;
+}
+
+/* ************************************************************************* */
+bool assert_equal(const MatrixColMajor& expected, const MatrixColMajor& actual, double tol) {
+ if (equal_with_abs_tol(expected,actual,tol)) return true;
+
+ size_t n1 = expected.cols(), m1 = expected.rows();
+ size_t n2 = actual.cols(), m2 = actual.rows();
+
+ cout << "not equal:" << endl;
+ print(expected,"expected = ");
+ print(actual,"actual = ");
+ if(m1!=m2 || n1!=n2)
+ cout << m1 << "," << n1 << " != " << m2 << "," << n2 << endl;
+ else {
+ Matrix diff = actual-expected;
+ print(diff, "actual - expected = ");
+ }
+ return false;
}
/* ************************************************************************* */
@@ -176,217 +141,226 @@ bool assert_equal(const std::list& As, const std::list& Bs, doub
/* ************************************************************************* */
static bool is_linear_dependent(const Matrix& A, const Matrix& B, double tol) {
- // This local static function is used by linear_independent and
- // linear_dependent just below.
- size_t n1 = A.size2(), m1 = A.size1();
- size_t n2 = B.size2(), m2 = B.size1();
+ // This local static function is used by linear_independent and
+ // linear_dependent just below.
+ size_t n1 = A.cols(), m1 = A.rows();
+ size_t n2 = B.cols(), m2 = B.rows();
- bool dependent = true;
- if(m1!=m2 || n1!=n2) dependent = false;
+ bool dependent = true;
+ if(m1!=m2 || n1!=n2) dependent = false;
- for(size_t i=0; dependent && i=A.size1())
- throw invalid_argument("Row index out of bounds!");
-
- const double * Aptr = A.data().begin() + A.size2() * i;
- return Vector_(A.size2(), Aptr);
+ size_t m = A.rows(), n = A.cols();
+ Vector v(m*n);
+ for( size_t j = 0, k=0 ; j < n ; j++)
+ for( size_t i = 0; i < m ; i++,k++)
+ v(k) = A(i,j);
+ return v;
}
/* ************************************************************************* */
void print(const Matrix& A, const string &s, ostream& stream) {
- size_t m = A.size1(), n = A.size2();
+ size_t m = A.rows(), n = A.cols();
- // print out all elements
- stream << s << "[\n";
- for( size_t i = 0 ; i < m ; i++) {
- for( size_t j = 0 ; j < n ; j++) {
- double aij = A(i,j);
- if(aij != 0.0)
- stream << setw(12) << setprecision(9) << aij << ",\t";
- else
- stream << " 0.0,\t";
- }
- stream << endl;
- }
- stream << "];" << endl;
+ // print out all elements
+ stream << s << "[\n";
+ for( size_t i = 0 ; i < m ; i++) {
+ for( size_t j = 0 ; j < n ; j++) {
+ double aij = A(i,j);
+ if(aij != 0.0)
+ stream << setw(12) << setprecision(9) << aij << ",\t";
+ else
+ stream << " 0.0,\t";
+ }
+ stream << endl;
+ }
+ stream << "];" << endl;
+}
+
+/* ************************************************************************* */
+void print(const MatrixColMajor& A, const string &s, ostream& stream) {
+ size_t m = A.rows(), n = A.cols();
+
+ // print out all elements
+ stream << s << "[\n";
+ for( size_t i = 0 ; i < m ; i++) {
+ for( size_t j = 0 ; j < n ; j++) {
+ double aij = A(i,j);
+ if(aij != 0.0)
+ stream << setw(12) << setprecision(9) << aij << ",\t";
+ else
+ stream << " 0.0,\t";
+ }
+ stream << endl;
+ }
+ stream << "];" << endl;
}
/* ************************************************************************* */
@@ -396,50 +370,26 @@ void save(const Matrix& A, const string &s, const string& filename) {
stream.close();
}
-/* ************************************************************************* */
-Matrix sub(const Matrix& A, size_t i1, size_t i2, size_t j1, size_t j2) {
-
- size_t m=i2-i1, n=j2-j1;
- Matrix B(m,n);
- for (size_t i=i1,k=0;i colproxy(A, j);
- colproxy = col;
+ A.col(j) = col;
}
/* ************************************************************************* */
void insertColumn(Matrix& A, const Vector& col, size_t i, size_t j) {
- ublas::matrix_column colproxy(A, j);
- ublas::vector_range > colsubproxy(colproxy,
- ublas::range (i, i+col.size()));
- colsubproxy = col;
+ A.col(j).segment(i, col.size()) = col;
}
/* ************************************************************************* */
-
Vector columnNormSquare(const MatrixColMajor &A) {
- Vector v (A.size2()) ;
- for ( size_t i = 0 ; i < A.size2() ; ++i ) {
- ublas::matrix_column mc (A, i);
- v[i] = dot(mc, mc) ;
+ Vector v (A.cols()) ;
+ for ( size_t i = 0 ; i < (size_t) A.cols() ; ++i ) {
+ v[i] = A.col(i).dot(A.col(i));
}
return v ;
}
@@ -447,67 +397,70 @@ Vector columnNormSquare(const MatrixColMajor &A) {
/* ************************************************************************* */
void solve(Matrix& A, Matrix& B)
{
- typedef ublas::permutation_matrix pmatrix;
+ // Eigen version - untested
+ Eigen::FullPivLU lu;
+ lu.compute(A);
+ B = lu.solve(B);
+ A = lu.matrixLU();
+
+// typedef ublas::permutation_matrix pmatrix;
// create a working copy of the input
- Matrix A_(A);
+// Matrix A_(A);
// create a permutation matrix for the LU-factorization
- pmatrix pm(A_.size1());
+// pmatrix pm(A_.rows());
// perform LU-factorization
- size_t res = lu_factorize(A_,pm);
- if( res != 0 ) throw runtime_error ("Matrix::solve: lu_factorize failed!");
+ // FIXME: add back with Eigen functionality
+ // size_t res = lu_factorize(A_,pm);
+ // if( res != 0 ) throw runtime_error ("Matrix::solve: lu_factorize failed!");
// backsubstitute to get the inverse
- lu_substitute(A_, pm, B);
+// lu_substitute(A_, pm, B);
}
/* ************************************************************************* */
-Matrix inverse(const Matrix& originalA)
+Matrix inverse(const Matrix& A)
{
- Matrix A(originalA);
- Matrix B = eye(A.size2());
- solve(A,B);
- return B;
+ return A.inverse();
}
/* ************************************************************************* */
/** Householder QR factorization, Golub & Van Loan p 224, explicit version */
/* ************************************************************************* */
pair qr(const Matrix& A) {
+ const size_t m = A.rows(), n = A.cols(), kprime = min(m,n);
- const size_t m = A.size1(), n = A.size2(), kprime = min(m,n);
-
- Matrix Q=eye(m,m),R(A);
- Vector v(m);
+ Matrix Q=eye(m,m),R(A);
+ Vector v(m);
- // loop over the kprime first columns
- for(size_t j=0; j < kprime; j++){
+ // loop over the kprime first columns
+ for(size_t j=0; j < kprime; j++){
- // we now work on the matrix (m-j)*(n-j) matrix A(j:end,j:end)
- const size_t mm=m-j;
+ // we now work on the matrix (m-j)*(n-j) matrix A(j:end,j:end)
+ const size_t mm=m-j;
- // copy column from matrix to xjm, i.e. x(j:m) = A(j:m,j)
- Vector xjm(mm);
- for(size_t k = 0 ; k < mm; k++)
- xjm(k) = R(j+k, j);
-
- // calculate the Householder vector v
- double beta; Vector vjm;
- boost::tie(beta,vjm) = house(xjm);
+ // copy column from matrix to xjm, i.e. x(j:m) = A(j:m,j)
+ Vector xjm(mm);
+ for(size_t k = 0 ; k < mm; k++)
+ xjm(k) = R(j+k, j);
- // pad with zeros to get m-dimensional vector v
- for(size_t k = 0 ; k < m; k++)
- v(k) = k qr(const Matrix& A) {
* on a number of different matrices for which all columns change.
*/
/* ************************************************************************* */
-inline void householder_update_manual(Matrix &A, size_t j, double beta, const Vector& vjm) {
- const size_t m = A.size1(), n = A.size2();
-
- Vector w(n);
- for( size_t c = 0; c < n; c++) {
- w(c) = 0.0;
- // dangerous as relies on row-major scheme
- const double *a = &A(j,c), * const v = &vjm(0);
- for( size_t r=j, s=0 ; r < m ; r++, s++, a+=n )
-
- w(c) += (*a) * v[s];
- w(c) *= beta;
- }
-
- // rank 1 update A(j:m,:) -= v(j:m)*w'
- for( size_t c = 0 ; c < n; c++) {
- double wc = w(c);
- double *a = &A(j,c); const double * const v =&vjm(0);
- for( size_t r=j, s=0 ; r < m ; r++, s++, a+=n )
-
- (*a) -= v[s] * wc;
- }
-}
-
-void householder_update(Matrix &A, size_t j, double beta, const Vector& vjm) {
-#if defined GT_USE_CBLAS
-
- // CBLAS version not working, using manual approach
- householder_update_manual(A,j,beta,vjm);
-
-
-
-#else
- householder_update_manual(A,j,beta,vjm);
-#endif
-}
-
-/* ************************************************************************* */
-// update A, b
-// A' \define A_{S}-ar and b'\define b-ad
-// __attribute__ ((noinline)) // uncomment to prevent inlining when profiling
-inline void updateAb_manual(Matrix& A, Vector& b, size_t j, const Vector& a,
- const Vector& r, double d) {
- const size_t m = A.size1(), n = A.size2();
- for (size_t i = 0; i < m; i++) { // update all rows
- double ai = a(i);
- b(i) -= ai * d;
- double *Aij = A.data().begin() + i * n + j + 1;
- const double *rptr = r.data().begin() + j + 1;
-
- for (size_t j2 = j + 1; j2 < n; j2++, Aij++, rptr++)
- *Aij -= ai * (*rptr);
- }
-}
+//void householder_update(Matrix &A, size_t j, double beta, const Vector& vjm) {
+// const size_t m = A.rows(), n = A.cols();
+//
+// Vector w(n);
+// for( size_t c = 0; c < n; c++) {
+// w(c) = 0.0;
+// // dangerous as relies on row-major scheme
+// const double *a = &A(j,c), * const v = &vjm(0);
+// for( size_t r=j, s=0 ; r < m ; r++, s++, a+=n )
+//
+// w(c) += (*a) * v[s];
+// w(c) *= beta;
+// }
+//
+// // rank 1 update A(j:m,:) -= v(j:m)*w'
+// for( size_t c = 0 ; c < n; c++) {
+// double wc = w(c);
+// double *a = &A(j,c); const double * const v =&vjm(0);
+// for( size_t r=j, s=0 ; r < m ; r++, s++, a+=n )
+// (*a) -= v[s] * wc;
+// }
+//}
/**
* Perform updates of system matrices
*/
static void updateAb(Matrix& A, Vector& b, size_t j, const Vector& a,
const Vector& r, double d) {
- // TODO: reimplement using BLAS
- updateAb_manual(A,b,j,a,r,d);
+ const size_t m = A.rows(), n = A.cols();
+ for (size_t i = 0; i < m; i++) { // update all rows
+ double ai = a(i);
+ b(i) -= ai * d;
+ double *Aij = A.data() + i * n + j + 1;
+ const double *rptr = r.data() + j + 1;
+
+ for (size_t j2 = j + 1; j2 < n; j2++, Aij++, rptr++)
+ *Aij -= ai * (*rptr);
+ }
}
/* ************************************************************************* */
list >
weighted_eliminate(Matrix& A, Vector& b, const Vector& sigmas) {
- size_t m = A.size1(), n = A.size2(); // get size(A)
+ size_t m = A.rows(), n = A.cols(); // get size(A)
size_t maxRank = min(m,n);
// create list
@@ -601,8 +530,8 @@ weighted_eliminate(Matrix& A, Vector& b, const Vector& sigmas) {
// Then update A and b by substituting x with d-rS, zero-ing out x's column.
for (size_t j=0; j(A, j2)); // TODO: don't use ublas
+ r(j2) = pseudo.dot(A.col(j2));
// create the rhs
- double d = inner_prod(pseudo, b);
+ double d = pseudo.dot(b);
// construct solution (r, d, sigma)
// TODO: avoid sqrt, store precision or at least variance
@@ -638,105 +567,105 @@ weighted_eliminate(Matrix& A, Vector& b, const Vector& sigmas) {
* version with Householder vectors below diagonal, as in GVL
*/
/* ************************************************************************* */
-inline void householder_manual(Matrix &A, size_t k) {
- const size_t m = A.size1(), n = A.size2(), kprime = min(k,min(m,n));
+void householder_(Matrix &A, size_t k)
+{
+ const size_t m = A.rows(), n = A.cols(), kprime = min(k,min(m,n));
// loop over the kprime first columns
for(size_t j=0; j < kprime; j++){
- // below, the indices r,c always refer to original A
-
- // copy column from matrix to xjm, i.e. x(j:m) = A(j:m,j)
- Vector xjm(m-j);
- for(size_t r = j ; r < m; r++)
- xjm(r-j) = A(r,j);
+ // copy column from matrix to vjm, i.e. v(j:m) = A(j:m,j)
+ Vector vjm = A.col(j).segment(j, m-j);
// calculate the Householder vector, in place
- double beta = houseInPlace(xjm);
- Vector& vjm = xjm;
+ double beta = houseInPlace(vjm);
- // do outer product update A = (I-beta vv')*A = A - v*(beta*A'*v)' = A - v*w'
- householder_update(A, j, beta, vjm);
+ // do outer product update A(j:m,:) = (I-beta vv')*A = A - v*(beta*A'*v)' = A - v*w'
+ Vector w = beta * A.middleRows(j,m-j).transpose() * vjm;
+ A.middleRows(j,m-j) -= vjm * w.transpose();
// the Householder vector is copied in the zeroed out part
- for( size_t r = j+1 ; r < m ; r++ )
- A(r,j) = vjm(r-j);
+ A.col(j).segment(j+1, m-(j+1)) = vjm.segment(1, m-(j+1));
+ } // column j
+// const size_t m = A.rows(), n = A.cols(), kprime = min(k,min(m,n));
+// // loop over the kprime first columns
+// for(size_t j=0; j < kprime; j++){
+// // below, the indices r,c always refer to original A
+//
+// // copy column from matrix to xjm, i.e. x(j:m) = A(j:m,j)
+// Vector xjm(m-j);
+// for(size_t r = j ; r < m; r++)
+// xjm(r-j) = A(r,j);
+//
+// // calculate the Householder vector, in place
+// double beta = houseInPlace(xjm);
+// Vector& vjm = xjm;
+//
+// // do outer product update A = (I-beta vv')*A = A - v*(beta*A'*v)' = A - v*w'
+// householder_update(A, j, beta, vjm);
+//
+// // the Householder vector is copied in the zeroed out part
+// for( size_t r = j+1 ; r < m ; r++ )
+// A(r,j) = vjm(r-j);
+//
+// } // column j
+}
+
+/* ************************************************************************* */
+void householder(Matrix &A, size_t k) {
+ // version with zeros below diagonal
+ householder_(A,k);
+ const size_t m = A.rows(), n = A.cols(), kprime = min(k,min(m,n));
+ for(size_t j=0; j < kprime; j++)
+ A.col(j).segment(j+1, m-(j+1)).setZero();
+// for( size_t i = j+1 ; i < m ; i++ )
+// A(i,j) = 0.0;
+}
+
+/* ************************************************************************* */
+void householder_(MatrixColMajor& A, size_t k, bool copy_vectors) {
+ const size_t m = A.rows(), n = A.cols(), kprime = min(k,min(m,n));
+ // loop over the kprime first columns
+ for(size_t j=0; j < kprime; j++) {
+ // copy column from matrix to vjm, i.e. v(j:m) = A(j:m,j)
+ Vector vjm = A.col(j).segment(j, m-j);
+
+ // calculate the Householder vector, in place
+ double beta = houseInPlace(vjm);
+
+ // do outer product update A(j:m,:) = (I-beta vv')*A = A - v*(beta*A'*v)' = A - v*w'
+ tic(1, "householder_update"); // bottleneck for system
+ // don't touch old columns
+ Vector w = beta * A.block(j,j,m-j,n-j).transpose() * vjm;
+ A.block(j,j,m-j,n-j) -= vjm * w.transpose();
+ toc(1, "householder_update");
+
+ // the Householder vector is copied in the zeroed out part
+ if (copy_vectors) {
+ tic(2, "householder_vector_copy");
+ A.col(j).segment(j+1, m-(j+1)) = vjm.segment(1, m-(j+1));
+ toc(2, "householder_vector_copy");
+ }
} // column j
}
-void householder_(Matrix &A, size_t k)
-{
- householder_manual(A, k);
+/* ************************************************************************* */
+void householder(MatrixColMajor& A, size_t k) {
+ // version with zeros below diagonal
+ tic(1, "householder_");
+ householder_(A,k,false);
+ toc(1, "householder_");
+// tic(2, "householder_zero_fill");
+// const size_t m = A.rows(), n = A.cols(), kprime = min(k,min(m,n));
+// for(size_t j=0; j < kprime; j++)
+// A.col(j).segment(j+1, m-(j+1)).setZero();
+// toc(2, "householder_zero_fill");
}
/* ************************************************************************* */
-/** version with zeros below diagonal */
-/* ************************************************************************* */
-void householder(Matrix &A, size_t k) {
- householder_(A,k);
- const size_t m = A.size1(), n = A.size2(), kprime = min(k,min(m,n));
- for(size_t j=0; j < kprime; j++)
- for( size_t i = j+1 ; i < m ; i++ )
- A(i,j) = 0.0;
-}
-
-/* ************************************************************************* */
-/** in-place householder */
-/* ************************************************************************* */
-#ifdef GT_USE_LAPACK
-#ifdef USE_LAPACK_QR
-void householder(Matrix &A) {
- __CLPK_integer m = A.size1();
- __CLPK_integer n = A.size2();
-
- // convert from row major to column major
- double a[m*n]; size_t k = 0;
- for(size_t j=0; j::epsilon()) {
- stringstream ss;
- ss << "backSubstituteUpper: L is singular,\n";
- print(L, "L: ", ss);
- throw invalid_argument(ss.str());
+ stringstream ss;
+ ss << "backSubstituteUpper: L is singular,\n";
+ print(L, "L: ", ss);
+ throw invalid_argument(ss.str());
}
#endif
if (!unit) zi /= L(i-1,i-1);
@@ -761,31 +690,92 @@ Vector backSubstituteLower(const Matrix& L, const Vector& b, bool unit) {
return result;
}
+/* ************************************************************************* */
+Vector backSubstituteUpper(const Matrix& U, const Vector& b, bool unit) {
+ // @return the solution x of U*x=b
+ size_t n = U.cols();
+#ifndef NDEBUG
+ size_t m = U.rows();
+ if (m!=n)
+ throw invalid_argument("backSubstituteUpper: U must be square");
+#endif
+
+ Vector result(n);
+ for (size_t i = n; i > 0; i--) {
+ double zi = b(i-1);
+ for (size_t j = i+1; j <= n; j++)
+ zi -= U(i-1,j-1) * result(j-1);
+#ifndef NDEBUG
+ if(!unit && fabs(U(i-1,i-1)) <= numeric_limits::epsilon()) {
+ stringstream ss;
+ ss << "backSubstituteUpper: U is singular,\n";
+ print(U, "U: ", ss);
+ throw invalid_argument(ss.str());
+ }
+#endif
+ if (!unit) zi /= U(i-1,i-1);
+ result(i-1) = zi;
+ }
+
+ return result;
+}
+
+/* ************************************************************************* */
+Vector backSubstituteUpper(const Vector& b, const Matrix& U, bool unit) {
+ // @return the solution x of x'*U=b'
+
+ size_t n = U.cols();
+#ifndef NDEBUG
+ size_t m = U.rows();
+ if (m!=n)
+ throw invalid_argument("backSubstituteUpper: U must be square");
+#endif
+
+ Vector result(n);
+ for (size_t i = 1; i <= n; i++) {
+ double zi = b(i-1);
+ for (size_t j = 1; j < i; j++)
+ zi -= U(j-1,i-1) * result(j-1);
+#ifndef NDEBUG
+ if(!unit && fabs(U(i-1,i-1)) <= numeric_limits::epsilon()) {
+ stringstream ss;
+ ss << "backSubstituteUpper: U is singular,\n";
+ print(U, "U: ", ss);
+ throw invalid_argument(ss.str());
+ }
+#endif
+ if (!unit) zi /= U(i-1,i-1);
+ result(i-1) = zi;
+ }
+
+ return result;
+}
+
/* ************************************************************************* */
Matrix stack(size_t nrMatrices, ...)
{
- size_t dimA1 = 0;
- size_t dimA2 = 0;
- va_list ap;
- va_start(ap, nrMatrices);
- for(size_t i = 0 ; i < nrMatrices ; i++) {
- Matrix *M = va_arg(ap, Matrix *);
- dimA1 += M->size1();
- dimA2 = M->size2(); // TODO: should check if all the same !
- }
- va_end(ap);
- va_start(ap, nrMatrices);
- Matrix A(dimA1, dimA2);
- size_t vindex = 0;
- for( size_t i = 0 ; i < nrMatrices ; i++) {
- Matrix *M = va_arg(ap, Matrix *);
- for(size_t d1 = 0; d1 < M->size1(); d1++)
- for(size_t d2 = 0; d2 < M->size2(); d2++)
- A(vindex+d1, d2) = (*M)(d1, d2);
- vindex += M->size1();
- }
+ size_t dimA1 = 0;
+ size_t dimA2 = 0;
+ va_list ap;
+ va_start(ap, nrMatrices);
+ for(size_t i = 0 ; i < nrMatrices ; i++) {
+ Matrix *M = va_arg(ap, Matrix *);
+ dimA1 += M->rows();
+ dimA2 = M->cols(); // TODO: should check if all the same !
+ }
+ va_end(ap);
+ va_start(ap, nrMatrices);
+ Matrix A(dimA1, dimA2);
+ size_t vindex = 0;
+ for( size_t i = 0 ; i < nrMatrices ; i++) {
+ Matrix *M = va_arg(ap, Matrix *);
+ for(size_t d1 = 0; d1 < (size_t) M->rows(); d1++)
+ for(size_t d2 = 0; d2 < (size_t) M->cols(); d2++)
+ A(vindex+d1, d2) = (*M)(d1, d2);
+ vindex += M->rows();
+ }
- return A;
+ return A;
}
/* ************************************************************************* */
@@ -796,30 +786,17 @@ Matrix collect(const std::vector& matrices, size_t m, size_t n)
size_t dimA2 = n*matrices.size();
if (!m && !n) {
BOOST_FOREACH(const Matrix* M, matrices) {
- dimA1 = M->size1(); // TODO: should check if all the same !
- dimA2 += M->size2();
+ dimA1 = M->rows(); // TODO: should check if all the same !
+ dimA2 += M->cols();
}
}
- // memcpy version
+ // stl::copy version
Matrix A(dimA1, dimA2);
- double * Aptr = A.data().begin();
size_t hindex = 0;
BOOST_FOREACH(const Matrix* M, matrices) {
- size_t row_len = M->size2();
-
- // find the size of the row to copy
- size_t row_size = sizeof(double) * row_len;
-
- // loop over rows
- for(size_t d1 = 0; d1 < M->size1(); ++d1) { // rows
- // get a pointer to the start of the row in each matrix
- double * Arow = Aptr + d1*dimA2 + hindex;
- double * Mrow = const_cast (M->data().begin() + d1*row_len);
-
- // do direct memory copy to move the row over
- memcpy(Arow, Mrow, row_size);
- }
+ size_t row_len = M->cols();
+ A.block(0, hindex, dimA1, row_len) = *M;
hindex += row_len;
}
@@ -829,23 +806,23 @@ Matrix collect(const std::vector& matrices, size_t m, size_t n)
/* ************************************************************************* */
Matrix collect(size_t nrMatrices, ...)
{
- vector matrices;
- va_list ap;
- va_start(ap, nrMatrices);
- for( size_t i = 0 ; i < nrMatrices ; i++) {
- Matrix *M = va_arg(ap, Matrix *);
- matrices.push_back(M);
- }
-return collect(matrices);
+ vector matrices;
+ va_list ap;
+ va_start(ap, nrMatrices);
+ for( size_t i = 0 ; i < nrMatrices ; i++) {
+ Matrix *M = va_arg(ap, Matrix *);
+ matrices.push_back(M);
+ }
+ return collect(matrices);
}
/* ************************************************************************* */
// row scaling, in-place
void vector_scale_inplace(const Vector& v, Matrix& A) {
- size_t m = A.size1(); size_t n = A.size2();
+ size_t m = A.rows(); size_t n = A.cols();
for (size_t i=0; i 0); // Rank failure
- double l_i_i = sqrt(p);
- L(i,i) = l_i_i;
-
- BNU::matrix_column l_i(L, i);
- project( l_i, BNU::range(i+1, A.size1()) )
- = ( BNU::project( BNU::column(A, i), BNU::range(i+1, A.size1()) )
- - BNU::prod( BNU::project(L, BNU::range(i+1, A.size1()), BNU::range(0, i)),
- BNU::project(BNU::row(L, i), BNU::range(0, i) ) ) ) / l_i_i;
- }
- return L;
+ Matrix L = zeros(A.rows(), A.rows());
+ Eigen::LLT llt;
+ llt.compute(A);
+ return llt.matrixL();
}
Matrix RtR(const Matrix &A)
{
- return trans(LLt(A));
+ return LLt(A).transpose();
}
/*
@@ -952,18 +910,20 @@ Matrix RtR(const Matrix &A)
*/
Matrix cholesky_inverse(const Matrix &A)
{
- Matrix L = LLt(A);
- Matrix inv(boost::numeric::ublas::identity_matrix(A.size1()));
- inplace_solve (L, inv, BNU::lower_tag ());
- return BNU::prod(trans(inv), inv);
-}
+ // FIXME: replace with real algorithm
+ return A.inverse();
+// Matrix L = LLt(A);
+// Matrix inv(eye(A.rows()));
+// inplace_solve (L, inv, BNU::lower_tag ());
+// return BNU::prod(trans(inv), inv);
+}
#if 0
/* ************************************************************************* */
// TODO, would be faster with Cholesky
Matrix inverse_square_root(const Matrix& A) {
- size_t m = A.size2(), n = A.size1();
+ size_t m = A.cols(), n = A.rows();
if (m!=n)
throw invalid_argument("inverse_square_root: A must be square");
@@ -974,8 +934,8 @@ Matrix inverse_square_root(const Matrix& A) {
// invert and sqrt diagonal of S
// We also arbitrarily choose sign to make result have positive signs
- for(size_t i = 0; i(A.size1()));
- inplace_solve (R, inv, BNU::upper_tag ());
+ Matrix inv = eye(A.rows());
+// inplace_solve(R, inv, BNU::upper_tag ());
+ R.triangularView().solveInPlace(inv);
return inv;
}
+/* ************************************************************************* */
+void svd(const Matrix& A, Matrix& U, Vector& S, Matrix& V) {
+ const size_t m = A.rows(), n = A.cols();
+ if (m < n) {
+ V = trans(A);
+ svd(V, S, U); // A'=V*diag(s)*U'
+ } else {
+ U = A; // copy
+ svd(U, S, V); // call in-place version
+ }
+}
+
+/* ************************************************************************* */
+void svd(Matrix& A, Vector& S, Matrix& V) {
+ Eigen::JacobiSVD svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV);
+ S = svd.singularValues();
+ A = svd.matrixU() * -1.0; // sign issues in tests - still valid, though
+ V = svd.matrixV() * -1.0;
+}
+
+/* ************************************************************************* */
+boost::tuple DLT(const Matrix& A, double rank_tol) {
+
+ // Check size of A
+ int m = A.rows(), n = A.cols();
+ if (m < n) throw invalid_argument(
+ "DLT: m rank_tol) rank++;
+ // Find minimum singular value and corresponding column index
+ int min_j = n - 1;
+ double min_S = S(min_j);
+ for (int j = 0; j < n - 1; j++)
+ if (S(j) < min_S) {
+ min_j = j;
+ min_S = S(j);
+ }
+
+ // Return rank, minimum singular value, and corresponding column of V
+ return boost::tuple(rank, min_S, Vector(column(V, min_j)));
+}
/* ************************************************************************* */
Matrix expm(const Matrix& A, size_t K) {
- Matrix E = eye(A.size1()), A_k = eye(A.size1());
+ Matrix E = eye(A.rows()), A_k = eye(A.rows());
for(size_t k=1;k<=K;k++) {
A_k = A_k*A/k;
E = E + A_k;
diff --git a/gtsam/base/Matrix.h b/gtsam/base/Matrix.h
index 7219e1f17..a0f204ce4 100644
--- a/gtsam/base/Matrix.h
+++ b/gtsam/base/Matrix.h
@@ -11,10 +11,11 @@
/**
* @file Matrix.h
- * @brief typedef and functions to augment Boost's ublas::matrix
+ * @brief typedef and functions to augment Eigen's MatrixXd
* @author Christian Potthast
* @author Kai Ni
* @author Frank Dellaert
+ * @author Alex Cunningham
*/
// \callgraph
@@ -23,20 +24,28 @@
#include
#include
-#include
#include
+#include
#include
/**
- * Vector is a *global* typedef
+ * Matrix is a *global* typedef
* wrap-matlab does this typedef as well
* we use the default < double,row_major,unbounded_array >
*/
-#if ! defined (MEX_H)
-typedef boost::numeric::ublas::matrix Matrix;
-typedef boost::numeric::ublas::matrix MatrixColMajor;
-#endif
+// FIXME: replace to handle matlab wrapper
+//#if ! defined (MEX_H)
+//typedef boost::numeric::ublas::matrix Matrix;
+//typedef boost::numeric::ublas::matrix MatrixColMajor;
+//#endif
+
+typedef Eigen::Matrix Matrix;
+typedef Eigen::Matrix MatrixColMajor;
+
+// Matrix expressions for accessing parts of matrices
+typedef Eigen::Block SubMatrix;
+typedef Eigen::Block ConstSubMatrix;
namespace gtsam {
@@ -72,7 +81,23 @@ Matrix diag(const Vector& v);
/**
* equals with an tolerance
*/
-bool equal_with_abs_tol(const Matrix& A, const Matrix& B, double tol = 1e-9);
+template
+bool equal_with_abs_tol(const Eigen::DenseBase& A, const Eigen::DenseBase& B, double tol = 1e-9) {
+
+ const size_t n1 = A.cols(), m1 = A.rows();
+ const size_t n2 = B.cols(), m2 = B.rows();
+
+ if(m1!=m2 || n1!=n2) return false;
+
+ for(size_t i=0; i tol)
+ return false;
+ }
+ return true;
+}
/**
* equality is just equal_with_abs_tol 1e-9
@@ -91,7 +116,9 @@ inline bool operator!=(const Matrix& A, const Matrix& B) {
/**
* equals with an tolerance, prints out message if unequal
*/
+// FIXME: make better use of templates to test these properly
bool assert_equal(const Matrix& A, const Matrix& B, double tol = 1e-9);
+bool assert_equal(const MatrixColMajor& A, const MatrixColMajor& B, double tol = 1e-9);
/**
* equals with an tolerance, prints out message if unequal
@@ -108,11 +135,6 @@ bool linear_independent(const Matrix& A, const Matrix& B, double tol = 1e-9);
*/
bool linear_dependent(const Matrix& A, const Matrix& B, double tol = 1e-9);
-/**
- * overload * for matrix-vector multiplication (as BOOST does not)
- */
-inline Vector operator*(const Matrix& A, const Vector & v) { return prod(A,v);}
-
/**
* BLAS Level-2 style e <- e + alpha*A*x
*/
@@ -144,15 +166,12 @@ void transposeMultiplyAdd(const Matrix& A, const Vector& e, Vector& x);
*/
void transposeMultiplyAdd(double alpha, const Matrix& A, const Vector& e, SubVector x);
-/**
- * overload * for vector*matrix multiplication (as BOOST does not)
- */
-inline Vector operator*(const Vector & v, const Matrix& A) { return prod(v,A);}
-
-/**
- * overload * for matrix multiplication (as BOOST does not)
- */
-inline Matrix operator*(const Matrix& A, const Matrix& B) { return prod(A,B);}
+/** products using old-style format to improve compatibility */
+template
+inline MATRIX prod(const MATRIX& A, const MATRIX&B) {
+ MATRIX result = A * B;
+ return result;
+}
/**
* convert to column vector, column order !!!
@@ -163,6 +182,7 @@ Vector Vector_(const Matrix& A);
* print a matrix
*/
void print(const Matrix& A, const std::string& s = "", std::ostream& stream = std::cout);
+void print(const MatrixColMajor& A, const std::string& s = "", std::ostream& stream = std::cout);
/**
* save a matrix to file, which can be loaded by matlab
@@ -178,7 +198,11 @@ void save(const Matrix& A, const std::string &s, const std::string& filename);
* @param j2 last col index + 1
* @return submatrix A(i1:i2-1,j1:j2-1)
*/
-Matrix sub(const Matrix& A, size_t i1, size_t i2, size_t j1, size_t j2);
+template
+Eigen::Block sub(const MATRIX& A, size_t i1, size_t i2, size_t j1, size_t j2) {
+ size_t m=i2-i1, n=j2-j1;
+ return A.block(i1,j1,m,n);
+}
/**
* insert a submatrix IN PLACE at a specified location in a larger matrix
@@ -191,13 +215,26 @@ Matrix sub(const Matrix& A, size_t i1, size_t i2, size_t j1, size_t j2);
void insertSub(Matrix& big, const Matrix& small, size_t i, size_t j);
/**
- * extracts a column from a matrix
- * NOTE: using this without the underscore is the ublas version!
+ * Extracts a column view from a matrix that avoids a copy
* @param matrix to extract column from
* @param index of the column
- * @return the column in vector form
+ * @return a const view of the matrix
*/
-Vector column_(const Matrix& A, size_t j);
+template
+const typename MATRIX::ConstColXpr column(const MATRIX& A, size_t j) {
+ return A.col(j);
+}
+
+/**
+ * Extracts a row view from a matrix that avoids a copy
+ * @param matrix to extract row from
+ * @param index of the row
+ * @return a const view of the matrix
+ */
+template
+const typename MATRIX::ConstRowXpr row(const MATRIX& A, size_t j) {
+ return A.row(j);
+}
/**
* inserts a column into a matrix IN PLACE
@@ -210,15 +247,25 @@ Vector column_(const Matrix& A, size_t j);
void insertColumn(Matrix& A, const Vector& col, size_t j);
void insertColumn(Matrix& A, const Vector& col, size_t i, size_t j);
-Vector columnNormSquare(const MatrixColMajor &A) ;
+Vector columnNormSquare(const MatrixColMajor &A);
/**
- * extracts a row from a matrix
- * @param matrix to extract row from
- * @param index of the row
- * @return the row in vector form
+ * Zeros all of the elements below the diagonal of a matrix, in place
+ * @param A is a matrix, to be modified in place
+ * @param cols is the number of columns to zero, use zero for all columns
*/
-Vector row_(const Matrix& A, size_t j);
+template
+void zeroBelowDiagonal(MATRIX& A, size_t cols=0) {
+ const size_t m = A.rows(), n = A.cols();
+ const size_t k = (cols) ? std::min(cols, std::min(m,n)) : std::min(m,n);
+ for (size_t j=0; j qr(const Matrix& A);
+/**
+ * QR factorization using Eigen's internal block QR algorithm
+ * @param A is the input matrix, and is the output
+ * @param clear_below_diagonal enables zeroing out below diagonal
+ */
+template
+void inplace_QR(MATRIX& A, bool clear_below_diagonal=true) {
+ size_t rows = A.rows();
+ size_t cols = A.cols();
+ size_t size = std::min(rows,cols);
+
+ typedef Eigen::internal::plain_diag_type::type HCoeffsType;
+ typedef Eigen::internal::plain_row_type::type RowVectorType;
+ HCoeffsType hCoeffs(size);
+ RowVectorType temp(cols);
+
+ Eigen::internal::householder_qr_inplace_blocked(A, hCoeffs, 48, temp.data());
+
+ zeroBelowDiagonal(A);
+}
+
/**
* Imperative version of Householder rank 1 update
*/
-void householder_update(Matrix &A, size_t j, double beta, const Vector& vjm);
+//void householder_update(Matrix &A, size_t j, double beta, const Vector& vjm);
/**
* Imperative algorithm for in-place full elimination with
@@ -271,24 +339,43 @@ void householder_(Matrix& A, size_t k);
*/
void householder(Matrix& A, size_t k);
-#ifdef GT_USE_LAPACK
-#ifdef USE_LAPACK_QR
/**
- * Householder tranformation, zeros below diagonal
+ * Householder tranformation, Householder vectors below diagonal
+ * Column Major Version
+ * @param k number of columns to zero out below diagonal
+ * @param A matrix
+ * @param copy_vectors - true to copy Householder vectors below diagonal
* @return nothing: in place !!!
*/
-void householder(Matrix& A);
+void householder_(MatrixColMajor& A, size_t k, bool copy_vectors=true);
/**
- * Householder tranformation directly on a column-major matrix. Does not zero
- * below the diagonal, so it will contain Householder vectors.
+ * Householder tranformation, zeros below diagonal
+ * Column Major version
+ * @param k number of columns to zero out below diagonal
+ * @param A matrix
* @return nothing: in place !!!
- * FIXME: this uses the LAPACK QR function, so it cannot be used on Ubuntu (without
- * lots of hacks, at least)
*/
-void householderColMajor(MatrixColMajor& A);
-#endif
-#endif
+void householder(MatrixColMajor& A, size_t k);
+
+//#ifdef GT_USE_LAPACK
+//#ifdef USE_LAPACK_QR
+///**
+// * Householder tranformation, zeros below diagonal
+// * @return nothing: in place !!!
+// */
+//void householder(Matrix& A);
+//
+///**
+// * Householder tranformation directly on a column-major matrix. Does not zero
+// * below the diagonal, so it will contain Householder vectors.
+// * @return nothing: in place !!!
+// * FIXME: this uses the LAPACK QR function, so it cannot be used on Ubuntu (without
+// * lots of hacks, at least)
+// */
+//void householderColMajor(MatrixColMajor& A);
+//#endif
+//#endif
/**
* backSubstitute U*x=b
@@ -296,11 +383,11 @@ void householderColMajor(MatrixColMajor& A);
* @param b an RHS vector
* @param unit, set true if unit triangular
* @return the solution x of U*x=b
- * TODO: use boost
*/
-template
-Vector backSubstituteUpper(const boost::numeric::ublas::matrix_expression& U,
- const boost::numeric::ublas::vector_expression& b, bool unit=false);
+//FIXME: add back expression form
+//template
+//Vector backSubstituteUpper(const MATRIX& U, const VECTOR& b, bool unit=false);
+Vector backSubstituteUpper(const Matrix& U, const Vector& b, bool unit=false);
/**
* backSubstitute x'*U=b'
@@ -308,11 +395,11 @@ Vector backSubstituteUpper(const boost::numeric::ublas::matrix_expression
-Vector backSubstituteUpper(const boost::numeric::ublas::vector_expression& b,
- const boost::numeric::ublas::matrix_expression& U, bool unit=false);
+//FIXME: add back expression form
+//template
+//Vector backSubstituteUpper(const VECTOR& b, const MATRIX& U, bool unit=false);
+Vector backSubstituteUpper(const Vector& b, const Matrix& U, bool unit=false);
/**
* backSubstitute L*x=b
@@ -320,7 +407,6 @@ Vector backSubstituteUpper(const boost::numeric::ublas::vector_expression n then U*S*V' = (m*n)*(n*n)*(n*n) (the m-n columns of U are of no use)
+ * if m < n then U*S*V' = (m*m)*(m*m)*(m*n) (the n-m columns of V are of no use)
+ * Careful! The dimensions above reflect V', not V, which is n*m if m n then U*S*V' = (m*n)*(n*n)*(n*n) (the m-n columns of U are of no use)
+ * You can just pass empty matrix V and vector S, they will be re-allocated.
+ */
+void svd(Matrix& A, Vector& S, Matrix& V);
+
+/**
+ * Direct linear transform algorithm that calls svd
+ * to find a vector v that minimizes the algebraic error A*v
+ * @param A of size m*n, where m>=n (pad with zero rows if not!)
+ * Returns rank of A, minimum error (singular value),
+ * and corresponding eigenvector (column of V, with A=U*S*V')
+ */
+boost::tuple
+DLT(const Matrix& A, double rank_tol = 1e-9);
+
/**
* Numerical exponential map, naive approach, not industrial strength !!!
* @param A matrix to exponentiate
@@ -392,3 +512,38 @@ Matrix expm(const Matrix& A, size_t K=7);
result_.addFailure(Failure(name_, __FILE__, __LINE__, #expected, #actual)); }
} // namespace gtsam
+
+#include
+#include
+
+namespace boost {
+namespace serialization {
+
+// split version - sends sizes ahead
+template
+void save(Archive & ar, const Matrix & m, unsigned int version)
+{
+ const int rows = m.rows(), cols = m.cols(), elements = rows*cols;
+ std::vector raw_data(elements);
+ std::copy(m.data(), m.data()+elements, raw_data.begin());
+ ar << make_nvp("rows", rows);
+ ar << make_nvp("cols", cols);
+ ar << make_nvp("data", raw_data);
+}
+template
+void load(Archive & ar, Matrix & m, unsigned int version)
+{
+ size_t rows, cols;
+ std::vector raw_data;
+ ar >> make_nvp("rows", rows);
+ ar >> make_nvp("cols", cols);
+ ar >> make_nvp("data", raw_data);
+ m = Matrix(rows, cols);
+ std::copy(raw_data.begin(), raw_data.end(), m.data());
+}
+
+} // namespace serialization
+} // namespace boost
+
+BOOST_SERIALIZATION_SPLIT_FREE(Matrix)
+
diff --git a/gtsam/base/Vector.cpp b/gtsam/base/Vector.cpp
index 835e9402a..7c526a626 100644
--- a/gtsam/base/Vector.cpp
+++ b/gtsam/base/Vector.cpp
@@ -10,11 +10,11 @@
* -------------------------------------------------------------------------- */
/**
-* @file Vector.cpp
-* @brief typedef and functions to augment Boost's ublas::vector
-* @author Kai Ni
-* @author Frank Dellaert
-*/
+ * @file Vector.cpp
+ * @brief typedef and functions to augment Boost's ublas::vector
+ * @author Kai Ni
+ * @author Frank Dellaert
+ */
#include
#include
@@ -31,460 +31,456 @@
#include
#endif
-#include
-#include
#include
#include
#include
using namespace std;
-namespace ublas = boost::numeric::ublas;
namespace gtsam {
-
- /* ************************************************************************* */
- void odprintf_(const char *format, ostream& stream, ...) {
- char buf[4096], *p = buf;
- int n;
-
- va_list args;
- va_start(args, stream);
- #ifdef WIN32
- n = _vsnprintf(p, sizeof buf - 3, format, args); // buf-3 is room for CR/LF/NUL
- #else
- n = vsnprintf(p, sizeof buf - 3, format, args); // buf-3 is room for CR/LF/NUL
- #endif
- va_end(args);
-
- #ifdef WIN32
- OutputDebugString(buf);
- #else
- stream << buf;
- #endif
- }
- /* ************************************************************************* */
-
- void odprintf(const char *format, ...)
- {
- char buf[4096], *p = buf;
- int n;
+/* ************************************************************************* */
+void odprintf_(const char *format, ostream& stream, ...) {
+ char buf[4096], *p = buf;
+ int n;
- va_list args;
- va_start(args, format);
- #ifdef WIN32
- n = _vsnprintf(p, sizeof buf - 3, format, args); // buf-3 is room for CR/LF/NUL
- #else
- n = vsnprintf(p, sizeof buf - 3, format, args); // buf-3 is room for CR/LF/NUL
- #endif
- va_end(args);
+ va_list args;
+ va_start(args, stream);
+#ifdef WIN32
+ n = _vsnprintf(p, sizeof buf - 3, format, args); // buf-3 is room for CR/LF/NUL
+#else
+ n = vsnprintf(p, sizeof buf - 3, format, args); // buf-3 is room for CR/LF/NUL
+#endif
+ va_end(args);
- #ifdef WIN32
- OutputDebugString(buf);
- #else
- cout << buf;
- #endif
- }
+#ifdef WIN32
+ OutputDebugString(buf);
+#else
+ stream << buf;
+#endif
+}
- /* ************************************************************************* */
- Vector Vector_( size_t m, const double* const data) {
- Vector v(m);
- copy(data, data+m, v.data().begin());
- return v;
- }
-
- /* ************************************************************************* */
- Vector Vector_(size_t m, ...) {
- Vector v(m);
- va_list ap;
- va_start(ap, m);
- for( size_t i = 0 ; i < m ; i++) {
- double value = va_arg(ap, double);
- v(i) = value;
- }
- va_end(ap);
- return v;
- }
-
- /* ************************************************************************* */
- bool zero(const Vector& v) {
- bool result = true;
- size_t n = v.size();
- for( size_t j = 0 ; j < n ; j++)
- result = result && (v(j) == 0.0);
- return result;
- }
-
- /* ************************************************************************* */
- Vector repeat(size_t n, double value) {
- Vector v(n, value);
- return v;
- }
+/* ************************************************************************* */
- /* ************************************************************************* */
- Vector delta(size_t n, size_t i, double value) {
- Vector v = zero(n);
- v(i) = value;
- return v;
- }
+void odprintf(const char *format, ...)
+{
+ char buf[4096], *p = buf;
+ int n;
- /* ************************************************************************* */
- void print(const Vector& v, const string& s, ostream& stream) {
- size_t n = v.size();
- odprintf_("%s [", stream, s.c_str());
- for(size_t i=0; i= it2[i]))
- return false;
- return true;
- }
+#ifdef WIN32
+ OutputDebugString(buf);
+#else
+ cout << buf;
+#endif
+}
- /* ************************************************************************* */
- bool equal_with_abs_tol(const Vector& vec1, const Vector& vec2, double tol) {
- Vector::const_iterator it1 = vec1.begin();
- Vector::const_iterator it2 = vec2.begin();
- if (vec1.size()!=vec2.size()) return false;
- for(size_t i=0; i tol)
- return false;
- }
- return true;
- }
-
- /* ************************************************************************* */
- bool assert_equal(const Vector& expected, const Vector& actual, double tol) {
- if (equal_with_abs_tol(expected,actual,tol)) return true;
- cout << "not equal:" << endl;
- print(expected, "expected");
- print(actual, "actual");
- return false;
- }
+/* ************************************************************************* */
+Vector Vector_( size_t m, const double* const data) {
+// Vector v(m);
+// for (size_t i=0; i scale;
- for(size_t i=0; itol&&fabs(it2[i])tol))
- return false;
- if(it1[i] == 0 && it2[i] == 0)
- continue;
- if (!scale)
- scale = it1[i] / it2[i];
- else if (fabs(it1[i] - it2[i] * (*scale)) > tol)
- return false;
- }
- return scale.is_initialized();
- }
+/* ************************************************************************* */
+Vector Vector_(const std::vector& d) {
+ Vector v(d.size());
+ copy(d.begin(), d.end(), v.data());
+ return v;
+}
- /* ************************************************************************* */
- Vector sub(const Vector &v, size_t i1, size_t i2) {
- size_t n = i2-i1;
- Vector v_return(n);
- for( size_t i = 0; i < n; i++ )
- v_return(i) = v(i1 + i);
- return v_return;
- }
-
- /* ************************************************************************* */
- void subInsert(Vector& big, const Vector& small, size_t i) {
- ublas::vector_range colsubproxy(big, ublas::range (i, i+small.size()));
- colsubproxy = small;
- }
+/* ************************************************************************* */
+bool zero(const Vector& v) {
+ bool result = true;
+ size_t n = v.size();
+ for( size_t j = 0 ; j < n ; j++)
+ result = result && (v(j) == 0.0);
+ return result;
+}
- /* ************************************************************************* */
- Vector emul(const Vector &a, const Vector &b) {
- size_t n = a.size();
- assert (b.size()==n);
- Vector c(n);
- for( size_t i = 0; i < n; i++ )
- c(i) = a(i)*b(i);
- return c;
+/* ************************************************************************* */
+Vector repeat(size_t n, double value) {
+ return Vector::Constant(n, value);
+}
+
+/* ************************************************************************* */
+Vector delta(size_t n, size_t i, double value) {
+ return Vector::Unit(n, i) * value;
+// Vector v = zero(n);
+// v(i) = value;
+// return v;
+}
+
+/* ************************************************************************* */
+void print(const Vector& v, const string& s, ostream& stream) {
+ size_t n = v.size();
+ odprintf_("%s [", stream, s.c_str());
+ for(size_t i=0; i= vec2[i]))
+ return false;
+ return true;
+}
+
+/* ************************************************************************* */
+bool equal_with_abs_tol(const Vector& vec1, const Vector& vec2, double tol) {
+ if (vec1.size()!=vec2.size()) return false;
+ size_t m = vec1.size();
+ for(size_t i=0; i tol)
+ return false;
+ }
+ return true;
+}
+
+/* ************************************************************************* */
+bool equal_with_abs_tol(const SubVector& vec1, const SubVector& vec2, double tol) {
+ if (vec1.size()!=vec2.size()) return false;
+ size_t m = vec1.size();
+ for(size_t i=0; i tol)
+ return false;
+ }
+ return true;
+}
+
+/* ************************************************************************* */
+bool assert_equal(const Vector& expected, const Vector& actual, double tol) {
+ if (equal_with_abs_tol(expected,actual,tol)) return true;
+ cout << "not equal:" << endl;
+ print(expected, "expected");
+ print(actual, "actual");
+ return false;
+}
+
+/* ************************************************************************* */
+bool assert_equal(SubVector expected, SubVector actual, double tol) {
+ if (equal_with_abs_tol(expected,actual,tol)) return true;
+ cout << "not equal:" << endl;
+ print(expected, "expected");
+ print(actual, "actual");
+ return false;
+}
+
+///* ************************************************************************* */
+//bool assert_equal(ConstSubVector expected, ConstSubVector actual, double tol) {
+// if (equal_with_abs_tol(Vector(expected),Vector(actual),tol)) return true;
+// cout << "not equal:" << endl;
+// print(Vector(expected), "expected");
+// print(Vector(actual), "actual");
+// return false;
+//}
+
+/* ************************************************************************* */
+bool linear_dependent(const Vector& vec1, const Vector& vec2, double tol) {
+ if (vec1.size()!=vec2.size()) return false;
+ boost::optional scale;
+ size_t m = vec1.size();
+ for(size_t i=0; itol&&fabs(vec2[i])tol))
+ return false;
+ if(vec1[i] == 0 && vec2[i] == 0)
+ continue;
+ if (!scale)
+ scale = vec1[i] / vec2[i];
+ else if (fabs(vec1[i] - vec2[i] * (*scale)) > tol)
+ return false;
+ }
+ return scale.is_initialized();
+}
+
+/* ************************************************************************* */
+ConstSubVector sub(const Vector &v, size_t i1, size_t i2) {
+ return v.segment(i1,i2-i1);
+}
+
+/* ************************************************************************* */
+void subInsert(Vector& big, const Vector& small, size_t i) {
+ big.segment(i, small.size()) = small;
+}
+
+/* ************************************************************************* */
+Vector emul(const Vector &a, const Vector &b) {
+ assert (b.size()==a.size());
+ return a.cwiseProduct(b);
+}
+
+/* ************************************************************************* */
+Vector ediv(const Vector &a, const Vector &b) {
+ assert (b.size()==a.size());
+ return a.cwiseQuotient(b);
+}
+
+/* ************************************************************************* */
+Vector ediv_(const Vector &a, const Vector &b) {
+ size_t n = a.size();
+ assert (b.size()==a.size());
+ Vector c(n);
+ for( size_t i = 0; i < n; i++ ) {
+ double ai = a(i), bi = b(i);
+ c(i) = (bi==0.0 && ai==0.0) ? 0.0 : a(i)/b(i);
+ }
+ return c;
+}
+
+/* ************************************************************************* */
+double sum(const Vector &a) {
+ return a.sum();
+// double result = 0.0;
+// size_t n = a.size();
+// for( size_t i = 0; i < n; i++ )
+// result += a(i);
+// return result;
+}
+
+/* ************************************************************************* */
+double norm_2(const Vector& v) {
+ return v.norm();
+}
+
+/* ************************************************************************* */
+Vector reciprocal(const Vector &a) {
+ size_t n = a.size();
+ Vector b(n);
+ for( size_t i = 0; i < n; i++ )
+ b(i) = 1.0/a(i);
+ return b;
+}
+
+/* ************************************************************************* */
+Vector esqrt(const Vector& v) {
+ return v.cwiseSqrt();
+// Vector s(v.size());
+// transform(v.begin(), v.end(), s.begin(), ::sqrt);
+// return s;
+}
+
+/* ************************************************************************* */
+Vector abs(const Vector& v) {
+ return v.cwiseAbs();
+}
+
+/* ************************************************************************* */
+double max(const Vector &a) {
+ return a.maxCoeff();
+}
+
+/* ************************************************************************* */
+double dot(const Vector& a, const Vector& b) {
+ assert (b.size()==a.size());
+ return a.dot(b);
+}
+
+/* ************************************************************************* */
+void scal(double alpha, Vector& x) {
+ x *= alpha;
+}
+
+/* ************************************************************************* */
+void axpy(double alpha, const Vector& x, Vector& y) {
+ assert (y.size()==x.size());
+ y += alpha * x;
+// size_t n = x.size();
+// for (size_t i = 0; i < n; i++)
+// y[i] += alpha * x[i];
+}
+
+/* ************************************************************************* */
+void axpy(double alpha, const Vector& x, SubVector y) {
+ assert (y.size()==x.size());
+ y += alpha * x;
+// size_t n = x.size();
+// for (size_t i = 0; i < n; i++)
+// y[i] += alpha * x[i];
+}
+
+/* ************************************************************************* */
+// imperative version, pass in x
+double houseInPlace(Vector &v) {
+ const double x0 = v(0);
+ const double x02 = x0*x0;
+
+ // old code - GSL verison was actually a bit slower
+ const double sigma = v.squaredNorm() - x02;
+
+ v(0) = 1.0;
+
+ if( sigma == 0.0 )
+ return 0.0;
+ else {
+ double mu = sqrt(x02 + sigma);
+ if( x0 <= 0.0 )
+ v(0) = x0 - mu;
+ else
+ v(0) = -sigma / (x0 + mu);
+
+ const double v02 = v(0)*v(0);
+ v = v / v(0);
+ return 2.0 * v02 / (sigma + v02);
+ }
+}
+
+/* ************************************************************************* */
+pair house(const Vector &x) {
+ Vector v(x);
+ double beta = houseInPlace(v);
+ return make_pair(beta, v);
+}
+
+/* ************************************************************************* */
+// Fast version *no error checking* !
+// Pass in initialized vector of size m or will crash !
+double weightedPseudoinverse(const Vector& a, const Vector& weights,
+ Vector& pseudo) {
+
+ size_t m = weights.size();
+ static const double inf = std::numeric_limits::infinity();
+
+ // Check once for zero entries of a. TODO: really needed ?
+ vector isZero;
+ for (size_t i = 0; i < m; ++i) isZero.push_back(fabs(a[i]) < 1e-9);
+
+ // If there is a valid (a!=0) constraint (sigma==0) return the first one
+ for (size_t i = 0; i < m; ++i)
+ if (weights[i] == inf && !isZero[i]) {
+ pseudo = delta(m, i, 1 / a[i]);
+ return inf;
}
- /* ************************************************************************* */
- Vector ediv(const Vector &a, const Vector &b) {
- size_t n = a.size();
- assert (b.size()==n);
- Vector c(n);
- for( size_t i = 0; i < n; i++ )
- c(i) = a(i)/b(i);
- return c;
- }
-
- /* ************************************************************************* */
- Vector ediv_(const Vector &a, const Vector &b) {
- size_t n = a.size();
- assert (b.size()==n);
- Vector c(n);
- for( size_t i = 0; i < n; i++ ) {
- double ai = a(i), bi = b(i);
- c(i) = (bi==0.0 && ai==0.0) ? 0.0 : a(i)/b(i);
- }
- return c;
- }
-
- /* ************************************************************************* */
- double sum(const Vector &a) {
- double result = 0.0;
- size_t n = a.size();
- for( size_t i = 0; i < n; i++ )
- result += a(i);
- return result;
- }
-
- /* ************************************************************************* */
- Vector reciprocal(const Vector &a) {
- size_t n = a.size();
- Vector b(n);
- for( size_t i = 0; i < n; i++ )
- b(i) = 1.0/a(i);
- return b;
- }
-
- /* ************************************************************************* */
- Vector esqrt(const Vector& v) {
- Vector s(v.size());
- transform(v.begin(), v.end(), s.begin(), ::sqrt);
- return s;
+ // Form psuedo-inverse inv(a'inv(Sigma)a)a'inv(Sigma)
+ // For diagonal Sigma, inv(Sigma) = diag(precisions)
+ double precision = 0;
+ for (size_t i = 0; i < m; i++) {
+ double ai = a[i];
+ if (!isZero[i]) // also catches remaining sigma==0 rows
+ precision += weights[i] * (ai * ai);
}
- /* ************************************************************************* */
- Vector abs(const Vector& v) {
- Vector s(v.size());
- transform(v.begin(), v.end(), s.begin(), ::fabs);
- return s;
+ // precision = a'inv(Sigma)a
+ if (precision < 1e-9)
+ for (size_t i = 0; i < m; i++)
+ pseudo[i] = 0;
+ else {
+ // emul(precisions,a)/precision
+ double variance = 1.0 / precision;
+ for (size_t i = 0; i < m; i++)
+ pseudo[i] = isZero[i] ? 0 : variance * weights[i] * a[i];
+ }
+ return precision; // sum of weights
+}
+
+/* ************************************************************************* */
+// Slow version with error checking
+pair
+weightedPseudoinverse(const Vector& a, const Vector& weights) {
+ int m = weights.size();
+ if (a.size() != m)
+ throw invalid_argument("a and weights have different sizes!");
+ Vector pseudo(m);
+ double precision = weightedPseudoinverse(a, weights, pseudo);
+ return make_pair(pseudo, precision);
+}
+
+/* ************************************************************************* */
+Vector concatVectors(const std::list& vs)
+{
+ size_t dim = 0;
+ BOOST_FOREACH(Vector v, vs)
+ dim += v.size();
+
+ Vector A(dim);
+ size_t index = 0;
+ BOOST_FOREACH(Vector v, vs) {
+ for(int d = 0; d < v.size(); d++)
+ A(d+index) = v(d);
+ index += v.size();
}
- /* ************************************************************************* */
- double max(const Vector &a) {
- return *(std::max_element(a.begin(), a.end()));
+ return A;
+}
+
+/* ************************************************************************* */
+Vector concatVectors(size_t nrVectors, ...)
+{
+ va_list ap;
+ list vs;
+ va_start(ap, nrVectors);
+ for( size_t i = 0 ; i < nrVectors ; i++) {
+ Vector* V = va_arg(ap, Vector*);
+ vs.push_back(*V);
}
+ va_end(ap);
+ return concatVectors(vs);
+}
- /* ************************************************************************* */
- double dot(const Vector& a, const Vector& b) {
- size_t n = a.size();
- assert (b.size()==n);
- double result = 0.0;
- for (size_t i = 0; i < n; i++)
- result += a[i] * b[i];
- return result;
- }
+/* ************************************************************************* */
+Vector rand_vector_norm(size_t dim, double mean, double sigma)
+{
+ boost::normal_distribution norm_dist(mean, sigma);
+ boost::variate_generator > norm(generator, norm_dist);
- /* ************************************************************************* */
- void scal(double alpha, Vector& x) {
- size_t n = x.size();
- for (size_t i = 0; i < n; i++)
- x[i] *= alpha;
- }
+ Vector v(dim);
+ for(int i = 0; i house(const Vector &x) {
- Vector v(x);
- double beta = houseInPlace(v);
- return make_pair(beta, v);
- }
-
- /* ************************************************************************* */
- // Fast version *no error checking* !
- // Pass in initialized vector of size m or will crash !
- double weightedPseudoinverse(const Vector& a, const Vector& weights,
- Vector& pseudo) {
-
- size_t m = weights.size();
- static const double inf = std::numeric_limits::infinity();
-
- // Check once for zero entries of a. TODO: really needed ?
- vector isZero;
- for (size_t i = 0; i < m; ++i) isZero.push_back(fabs(a[i]) < 1e-9);
-
- // If there is a valid (a!=0) constraint (sigma==0) return the first one
- for (size_t i = 0; i < m; ++i)
- if (weights[i] == inf && !isZero[i]) {
- pseudo = delta(m, i, 1 / a[i]);
- return inf;
- }
-
- // Form psuedo-inverse inv(a'inv(Sigma)a)a'inv(Sigma)
- // For diagonal Sigma, inv(Sigma) = diag(precisions)
- double precision = 0;
- for (size_t i = 0; i < m; i++) {
- double ai = a[i];
- if (!isZero[i]) // also catches remaining sigma==0 rows
- precision += weights[i] * (ai * ai);
- }
-
- // precision = a'inv(Sigma)a
- if (precision < 1e-9)
- for (size_t i = 0; i < m; i++)
- pseudo[i] = 0;
- else {
- // emul(precisions,a)/precision
- double variance = 1.0 / precision;
- for (size_t i = 0; i < m; i++)
- pseudo[i] = isZero[i] ? 0 : variance * weights[i] * a[i];
- }
- return precision; // sum of weights
- }
-
- /* ************************************************************************* */
- // Slow version with error checking
- pair
- weightedPseudoinverse(const Vector& a, const Vector& weights) {
- size_t m = weights.size();
- if (a.size() != m)
- throw invalid_argument("a and weights have different sizes!");
- Vector pseudo(m);
- double precision = weightedPseudoinverse(a, weights, pseudo);
- return make_pair(pseudo, precision);
- }
-
- /* ************************************************************************* */
- Vector concatVectors(const std::list& vs)
- {
- size_t dim = 0;
- BOOST_FOREACH(Vector v, vs)
- dim += v.size();
-
- Vector A(dim);
- size_t index = 0;
- BOOST_FOREACH(Vector v, vs) {
- for(size_t d = 0; d < v.size(); d++)
- A(d+index) = v(d);
- index += v.size();
- }
-
- return A;
- }
-
- /* ************************************************************************* */
- Vector concatVectors(size_t nrVectors, ...)
- {
- va_list ap;
- list vs;
- va_start(ap, nrVectors);
- for( size_t i = 0 ; i < nrVectors ; i++) {
- Vector* V = va_arg(ap, Vector*);
- vs.push_back(*V);
- }
- va_end(ap);
- return concatVectors(vs);
- }
-
- /* ************************************************************************* */
- Vector rand_vector_norm(size_t dim, double mean, double sigma)
- {
- boost::normal_distribution norm_dist(mean, sigma);
- boost::variate_generator > norm(generator, norm_dist);
-
- Vector v(dim);
- Vector::iterator it_v;
- for(it_v=v.begin(); it_v!=v.end(); it_v++)
- *it_v = norm();
-
- return v;
- }
-
- /* ************************************************************************* */
-
} // namespace gtsam
diff --git a/gtsam/base/Vector.h b/gtsam/base/Vector.h
index 707eeeec1..6005384f0 100644
--- a/gtsam/base/Vector.h
+++ b/gtsam/base/Vector.h
@@ -21,17 +21,22 @@
#pragma once
#include
-#include
-#include
+#include
+#include
#include
// Vector is a *global* typedef
// wrap-matlab does this typedef as well
-#if ! defined (MEX_H)
-typedef boost::numeric::ublas::vector Vector;
-#endif
-typedef boost::numeric::ublas::vector_range SubVector;
-typedef boost::numeric::ublas::vector_range ConstSubVector;
+// TODO: fix matlab wrapper
+//#if ! defined (MEX_H)
+//typedef boost::numeric::ublas::vector Vector;
+//#endif
+
+// Typedef arbitary length vector
+typedef Eigen::VectorXd Vector;
+
+typedef Eigen::VectorBlock SubVector;
+typedef Eigen::VectorBlock ConstSubVector;
namespace gtsam {
@@ -51,6 +56,11 @@ Vector Vector_( size_t m, const double* const data);
*/
Vector Vector_(size_t m, ...);
+/**
+ * Create a numeric vector from an STL vector of doubles
+ */
+Vector Vector_(const std::vector& data);
+
/**
* Create vector initialized to a constant value
* @param size
@@ -81,13 +91,13 @@ inline Vector basis(size_t n, size_t i) { return delta(n, i, 1.0); }
* Create zero vector
* @param size
*/
-inline Vector zero(size_t n) { return repeat(n,0.0);}
+inline Vector zero(size_t n) { return Vector::Zero(n);}
/**
* Create vector initialized to ones
* @param size
*/
-inline Vector ones(size_t n) { return repeat(n,1.0);}
+inline Vector ones(size_t n) { return Vector::Ones(n); }
/**
* check if all zero
@@ -125,6 +135,7 @@ bool greaterThanOrEqual(const Vector& v1, const Vector& v2);
* VecA == VecB up to tolerance
*/
bool equal_with_abs_tol(const Vector& vec1, const Vector& vec2, double tol=1e-9);
+bool equal_with_abs_tol(const SubVector& vec1, const SubVector& vec2, double tol=1e-9);
/**
* Override of equal in Lie.h
@@ -157,7 +168,7 @@ bool assert_equal(const Vector& vec1, const Vector& vec2, double tol=1e-9);
* @return bool
*/
bool assert_equal(SubVector vec1, SubVector vec2, double tol=1e-9);
-bool assert_equal(ConstSubVector vec1, ConstSubVector vec2, double tol=1e-9);
+//bool assert_equal(ConstSubVector vec1, ConstSubVector vec2, double tol=1e-9);
/**
* check whether two vectors are linearly dependent
@@ -175,7 +186,7 @@ bool linear_dependent(const Vector& vec1, const Vector& vec2, double tol=1e-9);
* @param i2 last row index + 1
* @return subvector v(i1:i2)
*/
-Vector sub(const Vector &v, size_t i1, size_t i2);
+ConstSubVector sub(const Vector &v, size_t i1, size_t i2);
/**
* Inserts a subvector into a vector IN PLACE
@@ -216,6 +227,14 @@ Vector ediv_(const Vector &a, const Vector &b);
*/
double sum(const Vector &a);
+/**
+ * Calculates L2 norm for a vector
+ * modeled after boost.ublas for compatibility
+ * @param vector
+ * @return the L2 norm
+ */
+double norm_2(const Vector& v);
+
/**
* elementwise reciprocal of vector elements
* @param a vector
@@ -247,6 +266,7 @@ double max(const Vector &a);
/** Dot product */
double dot(const Vector &a, const Vector& b);
+// TODO: remove simple blas functions - these are one-liners with Eigen
/**
* BLAS Level 1 scal: x <- alpha*x
*/
@@ -258,11 +278,6 @@ void scal(double alpha, Vector& x);
void axpy(double alpha, const Vector& x, Vector& y);
void axpy(double alpha, const Vector& x, SubVector y);
-/**
- * Divide every element of a Vector into a scalar
- */
-Vector operator/(double s, const Vector& v);
-
/**
* house(x,j) computes HouseHolder vector v and scaling factor beta
* from x, such that the corresponding Householder reflection zeroes out
@@ -310,5 +325,34 @@ Vector rand_vector_norm(size_t dim, double mean = 0, double sigma = 1);
} // namespace gtsam
+// FIXME: make this go away - use the Sampler class instead
static boost::minstd_rand generator(42u);
+
+#include
+#include
+
+namespace boost {
+namespace serialization {
+
+// split version - copies into an STL vector for serialization
+template
+void save(Archive & ar, const Vector & v, unsigned int version)
+{
+ const size_t n = v.size();
+ std::vector raw_data(n);
+ copy(v.data(), v.data()+n, raw_data.begin());
+ ar << make_nvp("data", raw_data);
+}
+template
+void load(Archive & ar, Vector & v, unsigned int version)
+{
+ std::vector raw_data;
+ ar >> make_nvp("data", raw_data);
+ v = gtsam::Vector_(raw_data);
+}
+
+} // namespace serialization
+} // namespace boost
+
+BOOST_SERIALIZATION_SPLIT_FREE(Vector)
diff --git a/gtsam/base/blockMatrices.h b/gtsam/base/blockMatrices.h
index 8bf8876a2..3523a84d6 100644
--- a/gtsam/base/blockMatrices.h
+++ b/gtsam/base/blockMatrices.h
@@ -18,60 +18,12 @@
#pragma once
#include
-#include
+#include
+#include
+#include
namespace gtsam {
-/** This is a wrapper around ublas::matrix_column that stores a copy of a
- * ublas::matrix_range. This does not copy the matrix data itself. The
- * purpose of this class is to be allow a column-of-a-range to be returned
- * from a function, given that a standard column-of-a-range just stores a
- * reference to the range. The stored range stores a reference to the original
- * matrix.
- */
-template
-class BlockColumn : public boost::numeric::ublas::vector_expression > {
-protected:
- typedef boost::numeric::ublas::matrix_range Range;
- typedef boost::numeric::ublas::matrix_column Base;
- Range range_;
- Base base_;
-public:
- typedef BlockColumn Self;
- typedef typename Base::matrix_type matrix_type;
- typedef typename Base::size_type size_type;
- typedef typename Base::difference_type difference_type;
- typedef typename Base::value_type value_type;
- typedef typename Base::const_reference const_reference;
- typedef typename Base::reference reference;
- typedef typename Base::storage_category storage_category;
- typedef Self closure_type;
- typedef const Self const_closure_type;
- typedef typename Base::iterator iterator;
- typedef typename Base::const_iterator const_iterator;
-
- BlockColumn(const boost::numeric::ublas::matrix_range& block, size_t column) :
- range_(block), base_(range_, column) {}
- BlockColumn(const BlockColumn& rhs) :
- range_(rhs.range_), base_(rhs.base_) {}
- BlockColumn& operator=(const BlockColumn& rhs) { base_.operator=(rhs.base_); return *this; }
- template BlockColumn& operator=(const boost::numeric::ublas::vector_expression& rhs) { base_.operator=(rhs); return *this; }
- typename Base::size_type size() const { return base_.size(); }
- const typename Base::matrix_closure_type& data() const { return base_.data(); }
- typename Base::matrix_closure_type& data() { return base_.data(); }
- typename Base::const_reference operator()(typename Base::size_type i) const { return base_(i); }
- typename Base::reference operator()(typename Base::size_type i) { return base_(i); }
- BlockColumn& assign_temporary(BlockColumn& rhs) { base_.assign_temporary(rhs.base_); return *this; }
- BlockColumn& assign_temporary(Base& rhs) { base_.assign_temporary(rhs); return *this; }
- bool same_closure(const BlockColumn& rhs) { return base_.same_closure(rhs.base_); }
- bool same_closure(const Base& rhs) { return base_.same_closure(rhs); }
- template BlockColumn& assign(const boost::numeric::ublas::vector_expression& rhs) { base_.assign(rhs); return *this; }
- iterator begin() { return base_.begin(); }
- const_iterator begin() const { return base_.begin(); }
- iterator end() { return base_.end(); }
- const_iterator end() const { return base_.end(); }
-};
-
template class SymmetricBlockView;
/**
@@ -95,10 +47,12 @@ template
class VerticalBlockView {
public:
typedef MATRIX FullMatrix;
- typedef typename boost::numeric::ublas::matrix_range Block;
- typedef typename boost::numeric::ublas::matrix_range constBlock;
- typedef BlockColumn Column;
- typedef BlockColumn constColumn;
+ typedef Eigen::Block Block;
+ typedef Eigen::Block constBlock;
+
+ // columns of blocks
+ typedef Eigen::VectorBlock Column;
+ typedef Eigen::VectorBlock constColumn;
protected:
FullMatrix& matrix_; // the reference to the full matrix
@@ -112,18 +66,19 @@ protected:
public:
/** Construct from an empty matrix (asserts that the matrix is empty) */
VerticalBlockView(FullMatrix& matrix) :
- matrix_(matrix), rowStart_(0), rowEnd_(matrix_.size1()), blockStart_(0) {
+ matrix_(matrix), rowStart_(0), rowEnd_(matrix_.rows()), blockStart_(0) {
fillOffsets((size_t*)0, (size_t*)0);
assertInvariants();
}
/**
* Construct from a non-empty matrix and copy the block structure from
- * another block view. */
+ * another block view.
+ */
template
VerticalBlockView(FullMatrix& matrix, const RHS& rhs) :
matrix_(matrix) {
- if(matrix_.size1() != rhs.size1() || matrix_.size2() != rhs.size2())
+ if((size_t) matrix_.rows() != rhs.rows() || (size_t) matrix_.cols() != rhs.cols())
throw std::invalid_argument(
"In VerticalBlockView<>(FullMatrix& matrix, const RHS& rhs), matrix and rhs must\n"
"already be of the same size. If not, construct the VerticalBlockView from an\n"
@@ -136,12 +91,13 @@ public:
/** Construct from iterators over the sizes of each vertical block */
template
VerticalBlockView(FullMatrix& matrix, ITERATOR firstBlockDim, ITERATOR lastBlockDim) :
- matrix_(matrix), rowStart_(0), rowEnd_(matrix_.size1()), blockStart_(0) {
+ matrix_(matrix), rowStart_(0), rowEnd_(matrix_.rows()), blockStart_(0) {
fillOffsets(firstBlockDim, lastBlockDim);
assertInvariants();
}
- /** Construct from a vector of the sizes of each vertical block, resize the
+ /**
+ * Construct from a vector of the sizes of each vertical block, resize the
* matrix so that its height is matrixNewHeight and its width fits the given
* block dimensions.
*/
@@ -149,17 +105,17 @@ public:
VerticalBlockView(FullMatrix& matrix, ITERATOR firstBlockDim, ITERATOR lastBlockDim, size_t matrixNewHeight) :
matrix_(matrix), rowStart_(0), rowEnd_(matrixNewHeight), blockStart_(0) {
fillOffsets(firstBlockDim, lastBlockDim);
- matrix_.resize(matrixNewHeight, variableColOffsets_.back(), false);
+ matrix_.resize(matrixNewHeight, variableColOffsets_.back());
assertInvariants();
}
/** Row size
*/
- size_t size1() const { assertInvariants(); return rowEnd_ - rowStart_; }
+ size_t rows() const { assertInvariants(); return rowEnd_ - rowStart_; }
/** Column size
*/
- size_t size2() const { assertInvariants(); return variableColOffsets_.back() - variableColOffsets_[blockStart_]; }
+ size_t cols() const { assertInvariants(); return variableColOffsets_.back() - variableColOffsets_[blockStart_]; }
/** Block count
@@ -167,60 +123,62 @@ public:
size_t nBlocks() const { assertInvariants(); return variableColOffsets_.size() - 1 - blockStart_; }
- Block operator()(size_t block) {
+ /** Access a single block in the underlying matrix with read/write access */
+ inline Block operator()(size_t block) {
return range(block, block+1);
}
- constBlock operator()(size_t block) const {
+ /** Access a const block view */
+ inline const constBlock operator()(size_t block) const {
return range(block, block+1);
}
- Block range(size_t startBlock, size_t endBlock) {
+ /** access ranges of blocks at a time */
+ inline Block range(size_t startBlock, size_t endBlock) {
assertInvariants();
size_t actualStartBlock = startBlock + blockStart_;
size_t actualEndBlock = endBlock + blockStart_;
checkBlock(actualStartBlock);
assert(actualEndBlock < variableColOffsets_.size());
- return Block(matrix_,
- boost::numeric::ublas::range(rowStart_, rowEnd_),
- boost::numeric::ublas::range(variableColOffsets_[actualStartBlock], variableColOffsets_[actualEndBlock]));
+ const size_t& startCol = variableColOffsets_[actualStartBlock];
+ const size_t& endCol = variableColOffsets_[actualEndBlock];
+ return matrix_.block(rowStart_, startCol, rowEnd_-rowStart_, endCol-startCol);
}
- constBlock range(size_t startBlock, size_t endBlock) const {
+ inline const constBlock range(size_t startBlock, size_t endBlock) const {
assertInvariants();
size_t actualStartBlock = startBlock + blockStart_;
size_t actualEndBlock = endBlock + blockStart_;
checkBlock(actualStartBlock);
assert(actualEndBlock < variableColOffsets_.size());
- return constBlock(matrix_,
- boost::numeric::ublas::range(rowStart_, rowEnd_),
- boost::numeric::ublas::range(variableColOffsets_[actualStartBlock], variableColOffsets_[actualEndBlock]));
+ const size_t& startCol = variableColOffsets_[actualStartBlock];
+ const size_t& endCol = variableColOffsets_[actualEndBlock];
+ return ((const FullMatrix&)matrix_).block(rowStart_, startCol, rowEnd_-rowStart_, endCol-startCol);
}
- Block full() {
+ inline Block full() {
return range(0,nBlocks());
}
- constBlock full() const {
+ inline const constBlock full() const {
return range(0,nBlocks());
}
+ /** get a single column out of a block */
Column column(size_t block, size_t columnOffset) {
assertInvariants();
size_t actualBlock = block + blockStart_;
checkBlock(actualBlock);
assert(variableColOffsets_[actualBlock] + columnOffset < variableColOffsets_[actualBlock+1]);
- Block blockMat(operator()(block));
- return Column(blockMat, columnOffset);
+ return matrix_.col(variableColOffsets_[actualBlock] + columnOffset).segment(rowStart_, rowEnd_-rowStart_);
}
- constColumn column(size_t block, size_t columnOffset) const {
+ const constColumn column(size_t block, size_t columnOffset) const {
assertInvariants();
size_t actualBlock = block + blockStart_;
checkBlock(actualBlock);
- assert(variableColOffsets_[actualBlock] + columnOffset < matrix_.size2());
- constBlock blockMat(operator()(block));
- return constColumn(blockMat, columnOffset);
+ assert(variableColOffsets_[actualBlock] + columnOffset < (size_t) matrix_.cols());
+ return ((const FullMatrix&)matrix_).col(variableColOffsets_[actualBlock] + columnOffset).segment(rowStart_, rowEnd_-rowStart_);
}
size_t offset(size_t block) const {
@@ -237,17 +195,21 @@ public:
size_t rowEnd() const { return rowEnd_; }
size_t firstBlock() const { return blockStart_; }
- /** Copy the block structure and resize the underlying matrix, but do not
+ /** access to full matrix */
+ const FullMatrix& fullMatrix() const { return matrix_; }
+
+ /**
+ * Copy the block structure and resize the underlying matrix, but do not
* copy the matrix data. If blockStart(), rowStart(), and/or rowEnd() have
* been modified, this copies the structure of the corresponding matrix view.
* In the destination VerticalBlockView, blockStart() and rowStart() will
- * thus be 0, rowEnd() will be size2() of the source VerticalBlockView, and
+ * thus be 0, rowEnd() will be cols() of the source VerticalBlockView, and
* the underlying matrix will be the size of the view of the source matrix.
*/
template
void copyStructureFrom(const RHS& rhs) {
- if(matrix_.size1() != rhs.size1() || matrix_.size2() != rhs.size2())
- matrix_.resize(rhs.size1(), rhs.size2(), false);
+ if((size_t) matrix_.rows() != (size_t) rhs.rows() || (size_t) matrix_.cols() != (size_t) rhs.cols())
+ matrix_.resize(rhs.rows(), rhs.cols());
if(rhs.blockStart_ == 0)
variableColOffsets_ = rhs.variableColOffsets_;
else {
@@ -261,7 +223,7 @@ public:
}
}
rowStart_ = 0;
- rowEnd_ = matrix_.size1();
+ rowEnd_ = matrix_.rows();
blockStart_ = 0;
assertInvariants();
}
@@ -275,13 +237,14 @@ public:
* If blockStart(), rowStart(), and/or rowEnd() have been modified, this
* copies the structure of the corresponding matrix view. In the destination
* VerticalBlockView, blockStart() and rowStart() will thus be 0, rowEnd()
- * will be size2() of the source VerticalBlockView, and the underlying matrix
+ * will be cols() of the source VerticalBlockView, and the underlying matrix
* will be the size of the view of the source matrix.
*/
template
VerticalBlockView