add lietorch src code and eigen src code, update readme

eigen-master/Eigen/SparseCore

+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_SPARSECORE_MODULE_H
+#define EIGEN_SPARSECORE_MODULE_H
+
+#include "Core"
+
+#include "src/Core/util/DisableStupidWarnings.h"
+
+#include <vector>
+#include <map>
+#include <cstdlib>
+#include <cstring>
+#include <algorithm>
+#include <numeric>
+
+/**
+ * \defgroup SparseCore_Module SparseCore module
+ *
+ * This module provides a sparse matrix representation, and basic associated matrix manipulations
+ * and operations.
+ *
+ * See the \ref TutorialSparse "Sparse tutorial"
+ *
+ * \code
+ * #include <Eigen/SparseCore>
+ * \endcode
+ *
+ * This module depends on: Core.
+ */
+
+// IWYU pragma: begin_exports
+#include "src/SparseCore/SparseUtil.h"
+#include "src/SparseCore/SparseMatrixBase.h"
+#include "src/SparseCore/SparseAssign.h"
+#include "src/SparseCore/CompressedStorage.h"
+#include "src/SparseCore/AmbiVector.h"
+#include "src/SparseCore/SparseCompressedBase.h"
+#include "src/SparseCore/SparseMatrix.h"
+#include "src/SparseCore/SparseMap.h"
+#include "src/SparseCore/SparseVector.h"
+#include "src/SparseCore/SparseRef.h"
+#include "src/SparseCore/SparseCwiseUnaryOp.h"
+#include "src/SparseCore/SparseCwiseBinaryOp.h"
+#include "src/SparseCore/SparseTranspose.h"
+#include "src/SparseCore/SparseBlock.h"
+#include "src/SparseCore/SparseDot.h"
+#include "src/SparseCore/SparseRedux.h"
+#include "src/SparseCore/SparseView.h"
+#include "src/SparseCore/SparseDiagonalProduct.h"
+#include "src/SparseCore/ConservativeSparseSparseProduct.h"
+#include "src/SparseCore/SparseSparseProductWithPruning.h"
+#include "src/SparseCore/SparseProduct.h"
+#include "src/SparseCore/SparseDenseProduct.h"
+#include "src/SparseCore/SparseSelfAdjointView.h"
+#include "src/SparseCore/SparseTriangularView.h"
+#include "src/SparseCore/TriangularSolver.h"
+#include "src/SparseCore/SparsePermutation.h"
+#include "src/SparseCore/SparseFuzzy.h"
+#include "src/SparseCore/SparseSolverBase.h"
+// IWYU pragma: end_exports
+
+#include "src/Core/util/ReenableStupidWarnings.h"
+
+#endif  // EIGEN_SPARSECORE_MODULE_H

eigen-master/Eigen/SparseLU

+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
+// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_SPARSELU_MODULE_H
+#define EIGEN_SPARSELU_MODULE_H
+
+#include "SparseCore"
+
+/**
+ * \defgroup SparseLU_Module SparseLU module
+ * This module defines a supernodal factorization of general sparse matrices.
+ * The code is fully optimized for supernode-panel updates with specialized kernels.
+ * Please, see the documentation of the SparseLU class for more details.
+ */
+
+// Ordering interface
+#include "OrderingMethods"
+
+#include "src/Core/util/DisableStupidWarnings.h"
+
+// IWYU pragma: begin_exports
+#include "src/SparseLU/SparseLU_Structs.h"
+#include "src/SparseLU/SparseLU_SupernodalMatrix.h"
+#include "src/SparseLU/SparseLUImpl.h"
+#include "src/SparseCore/SparseColEtree.h"
+#include "src/SparseLU/SparseLU_Memory.h"
+#include "src/SparseLU/SparseLU_heap_relax_snode.h"
+#include "src/SparseLU/SparseLU_relax_snode.h"
+#include "src/SparseLU/SparseLU_pivotL.h"
+#include "src/SparseLU/SparseLU_panel_dfs.h"
+#include "src/SparseLU/SparseLU_kernel_bmod.h"
+#include "src/SparseLU/SparseLU_panel_bmod.h"
+#include "src/SparseLU/SparseLU_column_dfs.h"
+#include "src/SparseLU/SparseLU_column_bmod.h"
+#include "src/SparseLU/SparseLU_copy_to_ucol.h"
+#include "src/SparseLU/SparseLU_pruneL.h"
+#include "src/SparseLU/SparseLU_Utils.h"
+#include "src/SparseLU/SparseLU.h"
+// IWYU pragma: end_exports
+
+#include "src/Core/util/ReenableStupidWarnings.h"
+
+#endif  // EIGEN_SPARSELU_MODULE_H

eigen-master/Eigen/SparseQR

+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_SPARSEQR_MODULE_H
+#define EIGEN_SPARSEQR_MODULE_H
+
+#include "SparseCore"
+#include "OrderingMethods"
+#include "src/Core/util/DisableStupidWarnings.h"
+
+/** \defgroup SparseQR_Module SparseQR module
+ * \brief Provides QR decomposition for sparse matrices
+ *
+ * This module provides a simplicial version of the left-looking Sparse QR decomposition.
+ * The columns of the input matrix should be reordered to limit the fill-in during the
+ * decomposition. Built-in methods (COLAMD, AMD) or external  methods (METIS) can be used to this end.
+ * See the \link OrderingMethods_Module OrderingMethods\endlink module for the list
+ * of built-in and external ordering methods.
+ *
+ * \code
+ * #include <Eigen/SparseQR>
+ * \endcode
+ *
+ *
+ */
+
+// IWYU pragma: begin_exports
+#include "src/SparseCore/SparseColEtree.h"
+#include "src/SparseQR/SparseQR.h"
+// IWYU pragma: end_exports
+
+#include "src/Core/util/ReenableStupidWarnings.h"
+
+#endif

eigen-master/Eigen/StdDeque

+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2009 Hauke Heibel <hauke.heibel@googlemail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_STDDEQUE_MODULE_H
+#define EIGEN_STDDEQUE_MODULE_H
+
+#include "Core"
+#include <deque>
+
+#if EIGEN_COMP_MSVC && EIGEN_OS_WIN64 && \
+    (EIGEN_MAX_STATIC_ALIGN_BYTES <= 16) /* MSVC auto aligns up to 16 bytes in 64 bit builds */
+
+#define EIGEN_DEFINE_STL_DEQUE_SPECIALIZATION(...)
+
+#else
+
+// IWYU pragma: begin_exports
+#include "src/StlSupport/StdDeque.h"
+// IWYU pragma: end_exports
+
+#endif
+
+#endif  // EIGEN_STDDEQUE_MODULE_H

eigen-master/Eigen/StdList

+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Hauke Heibel <hauke.heibel@googlemail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_STDLIST_MODULE_H
+#define EIGEN_STDLIST_MODULE_H
+
+#include "Core"
+#include <list>
+
+#if EIGEN_COMP_MSVC && EIGEN_OS_WIN64 && \
+    (EIGEN_MAX_STATIC_ALIGN_BYTES <= 16) /* MSVC auto aligns up to 16 bytes in 64 bit builds */
+
+#define EIGEN_DEFINE_STL_LIST_SPECIALIZATION(...)
+
+#else
+
+// IWYU pragma: begin_exports
+#include "src/StlSupport/StdList.h"
+// IWYU pragma: end_exports
+
+#endif
+
+#endif  // EIGEN_STDLIST_MODULE_H

eigen-master/Eigen/StdVector

+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2009 Hauke Heibel <hauke.heibel@googlemail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_STDVECTOR_MODULE_H
+#define EIGEN_STDVECTOR_MODULE_H
+
+#include "Core"
+#include <vector>
+
+#if EIGEN_COMP_MSVC && EIGEN_OS_WIN64 && \
+    (EIGEN_MAX_STATIC_ALIGN_BYTES <= 16) /* MSVC auto aligns up to 16 bytes in 64 bit builds */
+
+#define EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(...)
+
+#else
+
+// IWYU pragma: begin_exports
+#include "src/StlSupport/StdVector.h"
+// IWYU pragma: end_exports
+
+#endif
+
+#endif  // EIGEN_STDVECTOR_MODULE_H

eigen-master/Eigen/SuperLUSupport

+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_SUPERLUSUPPORT_MODULE_H
+#define EIGEN_SUPERLUSUPPORT_MODULE_H
+
+#include "SparseCore"
+
+#include "src/Core/util/DisableStupidWarnings.h"
+
+#ifdef EMPTY
+#define EIGEN_EMPTY_WAS_ALREADY_DEFINED
+#endif
+
+typedef int int_t;
+#include <slu_Cnames.h>
+#include <supermatrix.h>
+#include <slu_util.h>
+
+// slu_util.h defines a preprocessor token named EMPTY which is really polluting,
+// so we remove it in favor of a SUPERLU_EMPTY token.
+// If EMPTY was already defined then we don't undef it.
+
+#if defined(EIGEN_EMPTY_WAS_ALREADY_DEFINED)
+#undef EIGEN_EMPTY_WAS_ALREADY_DEFINED
+#elif defined(EMPTY)
+#undef EMPTY
+#endif
+
+#define SUPERLU_EMPTY (-1)
+
+namespace Eigen {
+struct SluMatrix;
+}
+
+/** \ingroup Support_modules
+ * \defgroup SuperLUSupport_Module SuperLUSupport module
+ *
+ * This module provides an interface to the <a href="http://crd-legacy.lbl.gov/~xiaoye/SuperLU/">SuperLU</a> library.
+ * It provides the following factorization class:
+ * - class SuperLU: a supernodal sequential LU factorization.
+ * - class SuperILU: a supernodal sequential incomplete LU factorization (to be used as a preconditioner for iterative
+ * methods).
+ *
+ * \warning This wrapper requires at least versions 4.0 of SuperLU. The 3.x versions are not supported.
+ *
+ * \warning When including this module, you have to use SUPERLU_EMPTY instead of EMPTY which is no longer defined
+ * because it is too polluting.
+ *
+ * \code
+ * #include <Eigen/SuperLUSupport>
+ * \endcode
+ *
+ * In order to use this module, the superlu headers must be accessible from the include paths, and your binary must be
+ * linked to the superlu library and its dependencies. The dependencies depend on how superlu has been compiled. For a
+ * cmake based project, you can use our FindSuperLU.cmake module to help you in this task.
+ *
+ */
+
+// IWYU pragma: begin_exports
+#include "src/SuperLUSupport/SuperLUSupport.h"
+// IWYU pragma: end_exports
+
+#include "src/Core/util/ReenableStupidWarnings.h"
+
+#endif  // EIGEN_SUPERLUSUPPORT_MODULE_H

eigen-master/Eigen/ThreadPool

+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2016 Benoit Steiner <benoit.steiner.goog@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_THREADPOOL_MODULE_H
+#define EIGEN_THREADPOOL_MODULE_H
+
+#include "Core"
+
+#include "src/Core/util/DisableStupidWarnings.h"
+
+/** \defgroup ThreadPool_Module ThreadPool Module
+ *
+ * This module provides 2 threadpool implementations
+ *  - a simple reference implementation
+ *  - a faster non blocking implementation
+ *
+ * \code
+ * #include <Eigen/ThreadPool>
+ * \endcode
+ */
+
+#include <cstddef>
+#include <cstring>
+#include <time.h>
+
+#include <vector>
+#include <atomic>
+#include <condition_variable>
+#include <deque>
+#include <mutex>
+#include <thread>
+#include <functional>
+#include <memory>
+#include <utility>
+
+// There are non-parenthesized calls to "max" in the  <unordered_map> header,
+// which trigger a check in test/main.h causing compilation to fail.
+// We work around the check here by removing the check for max in
+// the case where we have to emulate thread_local.
+#ifdef max
+#undef max
+#endif
+#include <unordered_map>
+
+#include "src/Core/util/Meta.h"
+#include "src/Core/util/MaxSizeVector.h"
+
+#ifndef EIGEN_MUTEX
+#define EIGEN_MUTEX std::mutex
+#endif
+#ifndef EIGEN_MUTEX_LOCK
+#define EIGEN_MUTEX_LOCK std::unique_lock<std::mutex>
+#endif
+#ifndef EIGEN_CONDVAR
+#define EIGEN_CONDVAR std::condition_variable
+#endif
+
+// IWYU pragma: begin_exports
+#include "src/ThreadPool/ThreadLocal.h"
+#include "src/ThreadPool/ThreadYield.h"
+#include "src/ThreadPool/ThreadCancel.h"
+#include "src/ThreadPool/EventCount.h"
+#include "src/ThreadPool/RunQueue.h"
+#include "src/ThreadPool/ThreadPoolInterface.h"
+#include "src/ThreadPool/ThreadEnvironment.h"
+#include "src/ThreadPool/Barrier.h"
+#include "src/ThreadPool/NonBlockingThreadPool.h"
+#include "src/ThreadPool/CoreThreadPoolDevice.h"
+#include "src/ThreadPool/ForkJoin.h"
+// IWYU pragma: end_exports
+
+#include "src/Core/util/ReenableStupidWarnings.h"
+
+#endif  // EIGEN_CXX11_THREADPOOL_MODULE_H

eigen-master/Eigen/UmfPackSupport

+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_UMFPACKSUPPORT_MODULE_H
+#define EIGEN_UMFPACKSUPPORT_MODULE_H
+
+#include "SparseCore"
+
+#include "src/Core/util/DisableStupidWarnings.h"
+
+extern "C" {
+#include <umfpack.h>
+}
+
+/** \ingroup Support_modules
+ * \defgroup UmfPackSupport_Module UmfPackSupport module
+ *
+ * This module provides an interface to the UmfPack library which is part of the <a
+ * href="http://www.suitesparse.com">suitesparse</a> package. It provides the following factorization class:
+ * - class UmfPackLU: a multifrontal sequential LU factorization.
+ *
+ * \code
+ * #include <Eigen/UmfPackSupport>
+ * \endcode
+ *
+ * In order to use this module, the umfpack headers must be accessible from the include paths, and your binary must be
+ * linked to the umfpack library and its dependencies. The dependencies depend on how umfpack has been compiled. For a
+ * cmake based project, you can use our FindUmfPack.cmake module to help you in this task.
+ *
+ */
+
+// IWYU pragma: begin_exports
+#include "src/UmfPackSupport/UmfPackSupport.h"
+// IWYU pragma: endexports
+
+#include "src/Core/util/ReenableStupidWarnings.h"
+
+#endif  // EIGEN_UMFPACKSUPPORT_MODULE_H

eigen-master/Eigen/src/AccelerateSupport/AccelerateSupport.h

+#ifndef EIGEN_ACCELERATESUPPORT_H
+#define EIGEN_ACCELERATESUPPORT_H
+
+#include <Accelerate/Accelerate.h>
+
+#include <Eigen/Sparse>
+
+namespace Eigen {
+
+template <typename MatrixType_, int UpLo_, SparseFactorization_t Solver_, bool EnforceSquare_>
+class AccelerateImpl;
+
+/** \ingroup AccelerateSupport_Module
+ * \typedef AccelerateLLT
+ * \brief A direct Cholesky (LLT) factorization and solver based on Accelerate
+ *
+ * \warning Only single and double precision real scalar types are supported by Accelerate
+ *
+ * \tparam MatrixType_ the type of the sparse matrix A, it must be a SparseMatrix<>
+ * \tparam UpLo_ additional information about the matrix structure. Default is Lower.
+ *
+ * \sa \ref TutorialSparseSolverConcept, class AccelerateLLT
+ */
+template <typename MatrixType, int UpLo = Lower>
+using AccelerateLLT = AccelerateImpl<MatrixType, UpLo | Symmetric, SparseFactorizationCholesky, true>;
+
+/** \ingroup AccelerateSupport_Module
+ * \typedef AccelerateLDLT
+ * \brief The default Cholesky (LDLT) factorization and solver based on Accelerate
+ *
+ * \warning Only single and double precision real scalar types are supported by Accelerate
+ *
+ * \tparam MatrixType_ the type of the sparse matrix A, it must be a SparseMatrix<>
+ * \tparam UpLo_ additional information about the matrix structure. Default is Lower.
+ *
+ * \sa \ref TutorialSparseSolverConcept, class AccelerateLDLT
+ */
+template <typename MatrixType, int UpLo = Lower>
+using AccelerateLDLT = AccelerateImpl<MatrixType, UpLo | Symmetric, SparseFactorizationLDLT, true>;
+
+/** \ingroup AccelerateSupport_Module
+ * \typedef AccelerateLDLTUnpivoted
+ * \brief A direct Cholesky-like LDL^T factorization and solver based on Accelerate with only 1x1 pivots and no pivoting
+ *
+ * \warning Only single and double precision real scalar types are supported by Accelerate
+ *
+ * \tparam MatrixType_ the type of the sparse matrix A, it must be a SparseMatrix<>
+ * \tparam UpLo_ additional information about the matrix structure. Default is Lower.
+ *
+ * \sa \ref TutorialSparseSolverConcept, class AccelerateLDLTUnpivoted
+ */
+template <typename MatrixType, int UpLo = Lower>
+using AccelerateLDLTUnpivoted = AccelerateImpl<MatrixType, UpLo | Symmetric, SparseFactorizationLDLTUnpivoted, true>;
+
+/** \ingroup AccelerateSupport_Module
+ * \typedef AccelerateLDLTSBK
+ * \brief A direct Cholesky (LDLT) factorization and solver based on Accelerate with Supernode Bunch-Kaufman and static
+ * pivoting
+ *
+ * \warning Only single and double precision real scalar types are supported by Accelerate
+ *
+ * \tparam MatrixType_ the type of the sparse matrix A, it must be a SparseMatrix<>
+ * \tparam UpLo_ additional information about the matrix structure. Default is Lower.
+ *
+ * \sa \ref TutorialSparseSolverConcept, class AccelerateLDLTSBK
+ */
+template <typename MatrixType, int UpLo = Lower>
+using AccelerateLDLTSBK = AccelerateImpl<MatrixType, UpLo | Symmetric, SparseFactorizationLDLTSBK, true>;
+
+/** \ingroup AccelerateSupport_Module
+ * \typedef AccelerateLDLTTPP
+ * \brief A direct Cholesky (LDLT) factorization and solver based on Accelerate with full threshold partial pivoting
+ *
+ * \warning Only single and double precision real scalar types are supported by Accelerate
+ *
+ * \tparam MatrixType_ the type of the sparse matrix A, it must be a SparseMatrix<>
+ * \tparam UpLo_ additional information about the matrix structure. Default is Lower.
+ *
+ * \sa \ref TutorialSparseSolverConcept, class AccelerateLDLTTPP
+ */
+template <typename MatrixType, int UpLo = Lower>
+using AccelerateLDLTTPP = AccelerateImpl<MatrixType, UpLo | Symmetric, SparseFactorizationLDLTTPP, true>;
+
+/** \ingroup AccelerateSupport_Module
+ * \typedef AccelerateQR
+ * \brief A QR factorization and solver based on Accelerate
+ *
+ * \warning Only single and double precision real scalar types are supported by Accelerate
+ *
+ * \tparam MatrixType_ the type of the sparse matrix A, it must be a SparseMatrix<>
+ *
+ * \sa \ref TutorialSparseSolverConcept, class AccelerateQR
+ */
+template <typename MatrixType>
+using AccelerateQR = AccelerateImpl<MatrixType, 0, SparseFactorizationQR, false>;
+
+/** \ingroup AccelerateSupport_Module
+ * \typedef AccelerateCholeskyAtA
+ * \brief A QR factorization and solver based on Accelerate without storing Q (equivalent to A^TA = R^T R)
+ *
+ * \warning Only single and double precision real scalar types are supported by Accelerate
+ *
+ * \tparam MatrixType_ the type of the sparse matrix A, it must be a SparseMatrix<>
+ *
+ * \sa \ref TutorialSparseSolverConcept, class AccelerateCholeskyAtA
+ */
+template <typename MatrixType>
+using AccelerateCholeskyAtA = AccelerateImpl<MatrixType, 0, SparseFactorizationCholeskyAtA, false>;
+
+namespace internal {
+template <typename T>
+struct AccelFactorizationDeleter {
+  void operator()(T* sym) {
+    if (sym) {
+      SparseCleanup(*sym);
+      delete sym;
+      sym = nullptr;
+    }
+  }
+};
+
+template <typename DenseVecT, typename DenseMatT, typename SparseMatT, typename NumFactT>
+struct SparseTypesTraitBase {
+  typedef DenseVecT AccelDenseVector;
+  typedef DenseMatT AccelDenseMatrix;
+  typedef SparseMatT AccelSparseMatrix;
+
+  typedef SparseOpaqueSymbolicFactorization SymbolicFactorization;
+  typedef NumFactT NumericFactorization;
+
+  typedef AccelFactorizationDeleter<SymbolicFactorization> SymbolicFactorizationDeleter;
+  typedef AccelFactorizationDeleter<NumericFactorization> NumericFactorizationDeleter;
+};
+
+template <typename Scalar>
+struct SparseTypesTrait {};
+
+template <>
+struct SparseTypesTrait<double> : SparseTypesTraitBase<DenseVector_Double, DenseMatrix_Double, SparseMatrix_Double,
+                                                       SparseOpaqueFactorization_Double> {};
+
+template <>
+struct SparseTypesTrait<float>
+    : SparseTypesTraitBase<DenseVector_Float, DenseMatrix_Float, SparseMatrix_Float, SparseOpaqueFactorization_Float> {
+};
+
+}  // end namespace internal
+
+template <typename MatrixType_, int UpLo_, SparseFactorization_t Solver_, bool EnforceSquare_>
+class AccelerateImpl : public SparseSolverBase<AccelerateImpl<MatrixType_, UpLo_, Solver_, EnforceSquare_> > {
+ protected:
+  using Base = SparseSolverBase<AccelerateImpl>;
+  using Base::derived;
+  using Base::m_isInitialized;
+
+ public:
+  using Base::_solve_impl;
+
+  typedef MatrixType_ MatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::StorageIndex StorageIndex;
+  enum { ColsAtCompileTime = Dynamic, MaxColsAtCompileTime = Dynamic };
+  enum { UpLo = UpLo_ };
+
+  using AccelDenseVector = typename internal::SparseTypesTrait<Scalar>::AccelDenseVector;
+  using AccelDenseMatrix = typename internal::SparseTypesTrait<Scalar>::AccelDenseMatrix;
+  using AccelSparseMatrix = typename internal::SparseTypesTrait<Scalar>::AccelSparseMatrix;
+  using SymbolicFactorization = typename internal::SparseTypesTrait<Scalar>::SymbolicFactorization;
+  using NumericFactorization = typename internal::SparseTypesTrait<Scalar>::NumericFactorization;
+  using SymbolicFactorizationDeleter = typename internal::SparseTypesTrait<Scalar>::SymbolicFactorizationDeleter;
+  using NumericFactorizationDeleter = typename internal::SparseTypesTrait<Scalar>::NumericFactorizationDeleter;
+
+  AccelerateImpl() {
+    m_isInitialized = false;
+
+    auto check_flag_set = [](int value, int flag) { return ((value & flag) == flag); };
+
+    if (check_flag_set(UpLo_, Symmetric)) {
+      m_sparseKind = SparseSymmetric;
+      m_triType = (UpLo_ & Lower) ? SparseLowerTriangle : SparseUpperTriangle;
+    } else if (check_flag_set(UpLo_, UnitLower)) {
+      m_sparseKind = SparseUnitTriangular;
+      m_triType = SparseLowerTriangle;
+    } else if (check_flag_set(UpLo_, UnitUpper)) {
+      m_sparseKind = SparseUnitTriangular;
+      m_triType = SparseUpperTriangle;
+    } else if (check_flag_set(UpLo_, StrictlyLower)) {
+      m_sparseKind = SparseTriangular;
+      m_triType = SparseLowerTriangle;
+    } else if (check_flag_set(UpLo_, StrictlyUpper)) {
+      m_sparseKind = SparseTriangular;
+      m_triType = SparseUpperTriangle;
+    } else if (check_flag_set(UpLo_, Lower)) {
+      m_sparseKind = SparseTriangular;
+      m_triType = SparseLowerTriangle;
+    } else if (check_flag_set(UpLo_, Upper)) {
+      m_sparseKind = SparseTriangular;
+      m_triType = SparseUpperTriangle;
+    } else {
+      m_sparseKind = SparseOrdinary;
+      m_triType = (UpLo_ & Lower) ? SparseLowerTriangle : SparseUpperTriangle;
+    }
+
+    m_order = SparseOrderDefault;
+  }
+
+  explicit AccelerateImpl(const MatrixType& matrix) : AccelerateImpl() { compute(matrix); }
+
+  ~AccelerateImpl() {}
+
+  inline Index cols() const { return m_nCols; }
+  inline Index rows() const { return m_nRows; }
+
+  ComputationInfo info() const {
+    eigen_assert(m_isInitialized && "Decomposition is not initialized.");
+    return m_info;
+  }
+
+  void analyzePattern(const MatrixType& matrix);
+
+  void factorize(const MatrixType& matrix);
+
+  void compute(const MatrixType& matrix);
+
+  template <typename Rhs, typename Dest>
+  void _solve_impl(const MatrixBase<Rhs>& b, MatrixBase<Dest>& dest) const;
+
+  /** Sets the ordering algorithm to use. */
+  void setOrder(SparseOrder_t order) { m_order = order; }
+
+ private:
+  template <typename T>
+  void buildAccelSparseMatrix(const SparseMatrix<T>& a, AccelSparseMatrix& A, std::vector<long>& columnStarts) {
+    const Index nColumnsStarts = a.cols() + 1;
+
+    columnStarts.resize(nColumnsStarts);
+
+    for (Index i = 0; i < nColumnsStarts; i++) columnStarts[i] = a.outerIndexPtr()[i];
+
+    SparseAttributes_t attributes{};
+    attributes.transpose = false;
+    attributes.triangle = m_triType;
+    attributes.kind = m_sparseKind;
+
+    SparseMatrixStructure structure{};
+    structure.attributes = attributes;
+    structure.rowCount = static_cast<int>(a.rows());
+    structure.columnCount = static_cast<int>(a.cols());
+    structure.blockSize = 1;
+    structure.columnStarts = columnStarts.data();
+    structure.rowIndices = const_cast<int*>(a.innerIndexPtr());
+
+    A.structure = structure;
+    A.data = const_cast<T*>(a.valuePtr());
+  }
+
+  void doAnalysis(AccelSparseMatrix& A) {
+    m_numericFactorization.reset(nullptr);
+
+    SparseSymbolicFactorOptions opts{};
+    opts.control = SparseDefaultControl;
+    opts.orderMethod = m_order;
+    opts.order = nullptr;
+    opts.ignoreRowsAndColumns = nullptr;
+    opts.malloc = malloc;
+    opts.free = free;
+    opts.reportError = nullptr;
+
+    m_symbolicFactorization.reset(new SymbolicFactorization(SparseFactor(Solver_, A.structure, opts)));
+
+    SparseStatus_t status = m_symbolicFactorization->status;
+
+    updateInfoStatus(status);
+
+    if (status != SparseStatusOK) m_symbolicFactorization.reset(nullptr);
+  }
+
+  void doFactorization(AccelSparseMatrix& A) {
+    SparseStatus_t status = SparseStatusReleased;
+
+    if (m_symbolicFactorization) {
+      m_numericFactorization.reset(new NumericFactorization(SparseFactor(*m_symbolicFactorization, A)));
+
+      status = m_numericFactorization->status;
+
+      if (status != SparseStatusOK) m_numericFactorization.reset(nullptr);
+    }
+
+    updateInfoStatus(status);
+  }
+
+ protected:
+  void updateInfoStatus(SparseStatus_t status) const {
+    switch (status) {
+      case SparseStatusOK:
+        m_info = Success;
+        break;
+      case SparseFactorizationFailed:
+      case SparseMatrixIsSingular:
+        m_info = NumericalIssue;
+        break;
+      case SparseInternalError:
+      case SparseParameterError:
+      case SparseStatusReleased:
+      default:
+        m_info = InvalidInput;
+        break;
+    }
+  }
+
+  mutable ComputationInfo m_info;
+  Index m_nRows, m_nCols;
+  std::unique_ptr<SymbolicFactorization, SymbolicFactorizationDeleter> m_symbolicFactorization;
+  std::unique_ptr<NumericFactorization, NumericFactorizationDeleter> m_numericFactorization;
+  SparseKind_t m_sparseKind;
+  SparseTriangle_t m_triType;
+  SparseOrder_t m_order;
+};
+
+/** Computes the symbolic and numeric decomposition of matrix \a a */
+template <typename MatrixType_, int UpLo_, SparseFactorization_t Solver_, bool EnforceSquare_>
+void AccelerateImpl<MatrixType_, UpLo_, Solver_, EnforceSquare_>::compute(const MatrixType& a) {
+  if (EnforceSquare_) eigen_assert(a.rows() == a.cols());
+
+  m_nRows = a.rows();
+  m_nCols = a.cols();
+
+  AccelSparseMatrix A{};
+  std::vector<long> columnStarts;
+
+  buildAccelSparseMatrix(a, A, columnStarts);
+
+  doAnalysis(A);
+
+  if (m_symbolicFactorization) doFactorization(A);
+
+  m_isInitialized = true;
+}
+
+/** Performs a symbolic decomposition on the sparsity pattern of matrix \a a.
+ *
+ * This function is particularly useful when solving for several problems having the same structure.
+ *
+ * \sa factorize()
+ */
+template <typename MatrixType_, int UpLo_, SparseFactorization_t Solver_, bool EnforceSquare_>
+void AccelerateImpl<MatrixType_, UpLo_, Solver_, EnforceSquare_>::analyzePattern(const MatrixType& a) {
+  if (EnforceSquare_) eigen_assert(a.rows() == a.cols());
+
+  m_nRows = a.rows();
+  m_nCols = a.cols();
+
+  AccelSparseMatrix A{};
+  std::vector<long> columnStarts;
+
+  buildAccelSparseMatrix(a, A, columnStarts);
+
+  doAnalysis(A);
+
+  m_isInitialized = true;
+}
+
+/** Performs a numeric decomposition of matrix \a a.
+ *
+ * The given matrix must have the same sparsity pattern as the matrix on which the symbolic decomposition has been
+ * performed.
+ *
+ * \sa analyzePattern()
+ */
+template <typename MatrixType_, int UpLo_, SparseFactorization_t Solver_, bool EnforceSquare_>
+void AccelerateImpl<MatrixType_, UpLo_, Solver_, EnforceSquare_>::factorize(const MatrixType& a) {
+  eigen_assert(m_symbolicFactorization && "You must first call analyzePattern()");
+  eigen_assert(m_nRows == a.rows() && m_nCols == a.cols());
+
+  if (EnforceSquare_) eigen_assert(a.rows() == a.cols());
+
+  AccelSparseMatrix A{};
+  std::vector<long> columnStarts;
+
+  buildAccelSparseMatrix(a, A, columnStarts);
+
+  doFactorization(A);
+}
+
+template <typename MatrixType_, int UpLo_, SparseFactorization_t Solver_, bool EnforceSquare_>
+template <typename Rhs, typename Dest>
+void AccelerateImpl<MatrixType_, UpLo_, Solver_, EnforceSquare_>::_solve_impl(const MatrixBase<Rhs>& b,
+                                                                              MatrixBase<Dest>& x) const {
+  if (!m_numericFactorization) {
+    m_info = InvalidInput;
+    return;
+  }
+
+  eigen_assert(m_nRows == b.rows());
+  eigen_assert(((b.cols() == 1) || b.outerStride() == b.rows()));
+
+  SparseStatus_t status = SparseStatusOK;
+
+  Scalar* b_ptr = const_cast<Scalar*>(b.derived().data());
+  Scalar* x_ptr = const_cast<Scalar*>(x.derived().data());
+
+  AccelDenseMatrix xmat{};
+  xmat.attributes = SparseAttributes_t();
+  xmat.columnCount = static_cast<int>(x.cols());
+  xmat.rowCount = static_cast<int>(x.rows());
+  xmat.columnStride = xmat.rowCount;
+  xmat.data = x_ptr;
+
+  AccelDenseMatrix bmat{};
+  bmat.attributes = SparseAttributes_t();
+  bmat.columnCount = static_cast<int>(b.cols());
+  bmat.rowCount = static_cast<int>(b.rows());
+  bmat.columnStride = bmat.rowCount;
+  bmat.data = b_ptr;
+
+  SparseSolve(*m_numericFactorization, bmat, xmat);
+
+  updateInfoStatus(status);
+}
+
+}  // end namespace Eigen
+
+#endif  // EIGEN_ACCELERATESUPPORT_H

eigen-master/Eigen/src/AccelerateSupport/InternalHeaderCheck.h

+#ifndef EIGEN_ACCELERATESUPPORT_MODULE_H
+#error "Please include Eigen/AccelerateSupport instead of including headers inside the src directory directly."
+#endif

eigen-master/Eigen/src/Cholesky/InternalHeaderCheck.h

+#ifndef EIGEN_CHOLESKY_MODULE_H
+#error "Please include Eigen/Cholesky instead of including headers inside the src directory directly."
+#endif

eigen-master/Eigen/src/Cholesky/LDLT.h

+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2009 Keir Mierle <mierle@gmail.com>
+// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2011 Timothy E. Holy <tim.holy@gmail.com >
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_LDLT_H
+#define EIGEN_LDLT_H
+
+// IWYU pragma: private
+#include "./InternalHeaderCheck.h"
+
+namespace Eigen {
+
+namespace internal {
+template <typename MatrixType_, int UpLo_>
+struct traits<LDLT<MatrixType_, UpLo_> > : traits<MatrixType_> {
+  typedef MatrixXpr XprKind;
+  typedef SolverStorage StorageKind;
+  typedef int StorageIndex;
+  enum { Flags = 0 };
+};
+
+template <typename MatrixType, int UpLo>
+struct LDLT_Traits;
+
+// PositiveSemiDef means positive semi-definite and non-zero; same for NegativeSemiDef
+enum SignMatrix { PositiveSemiDef, NegativeSemiDef, ZeroSign, Indefinite };
+}  // namespace internal
+
+/** \ingroup Cholesky_Module
+ *
+ * \class LDLT
+ *
+ * \brief Robust Cholesky decomposition of a matrix with pivoting
+ *
+ * \tparam MatrixType_ the type of the matrix of which to compute the LDL^T Cholesky decomposition
+ * \tparam UpLo_ the triangular part that will be used for the decomposition: Lower (default) or Upper.
+ *             The other triangular part won't be read.
+ *
+ * Perform a robust Cholesky decomposition of a positive semidefinite or negative semidefinite
+ * matrix \f$ A \f$ such that \f$ A =  P^TLDL^*P \f$, where P is a permutation matrix, L
+ * is lower triangular with a unit diagonal and D is a diagonal matrix.
+ *
+ * The decomposition uses pivoting to ensure stability, so that D will have
+ * zeros in the bottom right rank(A) - n submatrix. Avoiding the square root
+ * on D also stabilizes the computation.
+ *
+ * Remember that Cholesky decompositions are not rank-revealing. Also, do not use a Cholesky
+ * decomposition to determine whether a system of equations has a solution.
+ *
+ * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
+ *
+ * \sa MatrixBase::ldlt(), SelfAdjointView::ldlt(), class LLT
+ */
+template <typename MatrixType_, int UpLo_>
+class LDLT : public SolverBase<LDLT<MatrixType_, UpLo_> > {
+ public:
+  typedef MatrixType_ MatrixType;
+  typedef SolverBase<LDLT> Base;
+  friend class SolverBase<LDLT>;
+
+  EIGEN_GENERIC_PUBLIC_INTERFACE(LDLT)
+  enum {
+    MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
+    MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
+    UpLo = UpLo_
+  };
+  typedef Matrix<Scalar, RowsAtCompileTime, 1, 0, MaxRowsAtCompileTime, 1> TmpMatrixType;
+
+  typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
+  typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
+
+  typedef internal::LDLT_Traits<MatrixType, UpLo> Traits;
+
+  /** \brief Default Constructor.
+   *
+   * The default constructor is useful in cases in which the user intends to
+   * perform decompositions via LDLT::compute(const MatrixType&).
+   */
+  LDLT() : m_matrix(), m_transpositions(), m_sign(internal::ZeroSign), m_isInitialized(false) {}
+
+  /** \brief Default Constructor with memory preallocation
+   *
+   * Like the default constructor but with preallocation of the internal data
+   * according to the specified problem \a size.
+   * \sa LDLT()
+   */
+  explicit LDLT(Index size)
+      : m_matrix(size, size),
+        m_transpositions(size),
+        m_temporary(size),
+        m_sign(internal::ZeroSign),
+        m_isInitialized(false) {}
+
+  /** \brief Constructor with decomposition
+   *
+   * This calculates the decomposition for the input \a matrix.
+   *
+   * \sa LDLT(Index size)
+   */
+  template <typename InputType>
+  explicit LDLT(const EigenBase<InputType>& matrix)
+      : m_matrix(matrix.rows(), matrix.cols()),
+        m_transpositions(matrix.rows()),
+        m_temporary(matrix.rows()),
+        m_sign(internal::ZeroSign),
+        m_isInitialized(false) {
+    compute(matrix.derived());
+  }
+
+  /** \brief Constructs a LDLT factorization from a given matrix
+   *
+   * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c
+   * MatrixType is a Eigen::Ref.
+   *
+   * \sa LDLT(const EigenBase&)
+   */
+  template <typename InputType>
+  explicit LDLT(EigenBase<InputType>& matrix)
+      : m_matrix(matrix.derived()),
+        m_transpositions(matrix.rows()),
+        m_temporary(matrix.rows()),
+        m_sign(internal::ZeroSign),
+        m_isInitialized(false) {
+    compute(matrix.derived());
+  }
+
+  /** Clear any existing decomposition
+   * \sa rankUpdate(w,sigma)
+   */
+  void setZero() { m_isInitialized = false; }
+
+  /** \returns a view of the upper triangular matrix U */
+  inline typename Traits::MatrixU matrixU() const {
+    eigen_assert(m_isInitialized && "LDLT is not initialized.");
+    return Traits::getU(m_matrix);
+  }
+
+  /** \returns a view of the lower triangular matrix L */
+  inline typename Traits::MatrixL matrixL() const {
+    eigen_assert(m_isInitialized && "LDLT is not initialized.");
+    return Traits::getL(m_matrix);
+  }
+
+  /** \returns the permutation matrix P as a transposition sequence.
+   */
+  inline const TranspositionType& transpositionsP() const {
+    eigen_assert(m_isInitialized && "LDLT is not initialized.");
+    return m_transpositions;
+  }
+
+  /** \returns the coefficients of the diagonal matrix D */
+  inline Diagonal<const MatrixType> vectorD() const {
+    eigen_assert(m_isInitialized && "LDLT is not initialized.");
+    return m_matrix.diagonal();
+  }
+
+  /** \returns true if the matrix is positive (semidefinite) */
+  inline bool isPositive() const {
+    eigen_assert(m_isInitialized && "LDLT is not initialized.");
+    return m_sign == internal::PositiveSemiDef || m_sign == internal::ZeroSign;
+  }
+
+  /** \returns true if the matrix is negative (semidefinite) */
+  inline bool isNegative(void) const {
+    eigen_assert(m_isInitialized && "LDLT is not initialized.");
+    return m_sign == internal::NegativeSemiDef || m_sign == internal::ZeroSign;
+  }
+
+#ifdef EIGEN_PARSED_BY_DOXYGEN
+  /** \returns a solution x of \f$ A x = b \f$ using the current decomposition of A.
+   *
+   * This function also supports in-place solves using the syntax <tt>x = decompositionObject.solve(x)</tt> .
+   *
+   * \note_about_checking_solutions
+   *
+   * More precisely, this method solves \f$ A x = b \f$ using the decomposition \f$ A = P^T L D L^* P \f$
+   * by solving the systems \f$ P^T y_1 = b \f$, \f$ L y_2 = y_1 \f$, \f$ D y_3 = y_2 \f$,
+   * \f$ L^* y_4 = y_3 \f$ and \f$ P x = y_4 \f$ in succession. If the matrix \f$ A \f$ is singular, then
+   * \f$ D \f$ will also be singular (all the other matrices are invertible). In that case, the
+   * least-square solution of \f$ D y_3 = y_2 \f$ is computed. This does not mean that this function
+   * computes the least-square solution of \f$ A x = b \f$ if \f$ A \f$ is singular.
+   *
+   * \sa MatrixBase::ldlt(), SelfAdjointView::ldlt()
+   */
+  template <typename Rhs>
+  inline const Solve<LDLT, Rhs> solve(const MatrixBase<Rhs>& b) const;
+#endif
+
+  template <typename Derived>
+  bool solveInPlace(MatrixBase<Derived>& bAndX) const;
+
+  template <typename InputType>
+  LDLT& compute(const EigenBase<InputType>& matrix);
+
+  /** \returns an estimate of the reciprocal condition number of the matrix of
+   *  which \c *this is the LDLT decomposition.
+   */
+  RealScalar rcond() const {
+    eigen_assert(m_isInitialized && "LDLT is not initialized.");
+    return internal::rcond_estimate_helper(m_l1_norm, *this);
+  }
+
+  template <typename Derived>
+  LDLT& rankUpdate(const MatrixBase<Derived>& w, const RealScalar& alpha = 1);
+
+  /** \returns the internal LDLT decomposition matrix
+   *
+   * TODO: document the storage layout
+   */
+  inline const MatrixType& matrixLDLT() const {
+    eigen_assert(m_isInitialized && "LDLT is not initialized.");
+    return m_matrix;
+  }
+
+  MatrixType reconstructedMatrix() const;
+
+  /** \returns the adjoint of \c *this, that is, a const reference to the decomposition itself as the underlying matrix
+   * is self-adjoint.
+   *
+   * This method is provided for compatibility with other matrix decompositions, thus enabling generic code such as:
+   * \code x = decomposition.adjoint().solve(b) \endcode
+   */
+  const LDLT& adjoint() const { return *this; }
+
+  EIGEN_DEVICE_FUNC constexpr Index rows() const noexcept { return m_matrix.rows(); }
+  EIGEN_DEVICE_FUNC constexpr Index cols() const noexcept { return m_matrix.cols(); }
+
+  /** \brief Reports whether previous computation was successful.
+   *
+   * \returns \c Success if computation was successful,
+   *          \c NumericalIssue if the factorization failed because of a zero pivot.
+   */
+  ComputationInfo info() const {
+    eigen_assert(m_isInitialized && "LDLT is not initialized.");
+    return m_info;
+  }
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+  template <typename RhsType, typename DstType>
+  void _solve_impl(const RhsType& rhs, DstType& dst) const;
+
+  template <bool Conjugate, typename RhsType, typename DstType>
+  void _solve_impl_transposed(const RhsType& rhs, DstType& dst) const;
+#endif
+
+ protected:
+  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
+
+  /** \internal
+   * Used to compute and store the Cholesky decomposition A = L D L^* = U^* D U.
+   * The strict upper part is used during the decomposition, the strict lower
+   * part correspond to the coefficients of L (its diagonal is equal to 1 and
+   * is not stored), and the diagonal entries correspond to D.
+   */
+  MatrixType m_matrix;
+  RealScalar m_l1_norm;
+  TranspositionType m_transpositions;
+  TmpMatrixType m_temporary;
+  internal::SignMatrix m_sign;
+  bool m_isInitialized;
+  ComputationInfo m_info;
+};
+
+namespace internal {
+
+template <int UpLo>
+struct ldlt_inplace;
+
+template <>
+struct ldlt_inplace<Lower> {
+  template <typename MatrixType, typename TranspositionType, typename Workspace>
+  static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign) {
+    using std::abs;
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename MatrixType::RealScalar RealScalar;
+    typedef typename TranspositionType::StorageIndex IndexType;
+    eigen_assert(mat.rows() == mat.cols());
+    const Index size = mat.rows();
+    bool found_zero_pivot = false;
+    bool ret = true;
+
+    if (size <= 1) {
+      transpositions.setIdentity();
+      if (size == 0)
+        sign = ZeroSign;
+      else if (numext::real(mat.coeff(0, 0)) > static_cast<RealScalar>(0))
+        sign = PositiveSemiDef;
+      else if (numext::real(mat.coeff(0, 0)) < static_cast<RealScalar>(0))
+        sign = NegativeSemiDef;
+      else
+        sign = ZeroSign;
+      return true;
+    }
+
+    for (Index k = 0; k < size; ++k) {
+      // Find largest diagonal element
+      Index index_of_biggest_in_corner;
+      mat.diagonal().tail(size - k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner);
+      index_of_biggest_in_corner += k;
+
+      transpositions.coeffRef(k) = IndexType(index_of_biggest_in_corner);
+      if (k != index_of_biggest_in_corner) {
+        // apply the transposition while taking care to consider only
+        // the lower triangular part
+        Index s = size - index_of_biggest_in_corner - 1;  // trailing size after the biggest element
+        mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
+        mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
+        std::swap(mat.coeffRef(k, k), mat.coeffRef(index_of_biggest_in_corner, index_of_biggest_in_corner));
+        for (Index i = k + 1; i < index_of_biggest_in_corner; ++i) {
+          Scalar tmp = mat.coeffRef(i, k);
+          mat.coeffRef(i, k) = numext::conj(mat.coeffRef(index_of_biggest_in_corner, i));
+          mat.coeffRef(index_of_biggest_in_corner, i) = numext::conj(tmp);
+        }
+        if (NumTraits<Scalar>::IsComplex)
+          mat.coeffRef(index_of_biggest_in_corner, k) = numext::conj(mat.coeff(index_of_biggest_in_corner, k));
+      }
+
+      // partition the matrix:
+      //       A00 |  -  |  -
+      // lu  = A10 | A11 |  -
+      //       A20 | A21 | A22
+      Index rs = size - k - 1;
+      Block<MatrixType, Dynamic, 1> A21(mat, k + 1, k, rs, 1);
+      Block<MatrixType, 1, Dynamic> A10(mat, k, 0, 1, k);
+      Block<MatrixType, Dynamic, Dynamic> A20(mat, k + 1, 0, rs, k);
+
+      if (k > 0) {
+        temp.head(k) = mat.diagonal().real().head(k).asDiagonal() * A10.adjoint();
+        mat.coeffRef(k, k) -= (A10 * temp.head(k)).value();
+        if (rs > 0) A21.noalias() -= A20 * temp.head(k);
+      }
+
+      // In some previous versions of Eigen (e.g., 3.2.1), the scaling was omitted if the pivot
+      // was smaller than the cutoff value. However, since LDLT is not rank-revealing
+      // we should only make sure that we do not introduce INF or NaN values.
+      // Remark that LAPACK also uses 0 as the cutoff value.
+      RealScalar realAkk = numext::real(mat.coeffRef(k, k));
+      bool pivot_is_valid = (abs(realAkk) > RealScalar(0));
+
+      if (k == 0 && !pivot_is_valid) {
+        // The entire diagonal is zero, there is nothing more to do
+        // except filling the transpositions, and checking whether the matrix is zero.
+        sign = ZeroSign;
+        for (Index j = 0; j < size; ++j) {
+          transpositions.coeffRef(j) = IndexType(j);
+          ret = ret && (mat.col(j).tail(size - j - 1).array() == Scalar(0)).all();
+        }
+        return ret;
+      }
+
+      if ((rs > 0) && pivot_is_valid)
+        A21 /= realAkk;
+      else if (rs > 0)
+        ret = ret && (A21.array() == Scalar(0)).all();
+
+      if (found_zero_pivot && pivot_is_valid)
+        ret = false;  // factorization failed
+      else if (!pivot_is_valid)
+        found_zero_pivot = true;
+
+      if (sign == PositiveSemiDef) {
+        if (realAkk < static_cast<RealScalar>(0)) sign = Indefinite;
+      } else if (sign == NegativeSemiDef) {
+        if (realAkk > static_cast<RealScalar>(0)) sign = Indefinite;
+      } else if (sign == ZeroSign) {
+        if (realAkk > static_cast<RealScalar>(0))
+          sign = PositiveSemiDef;
+        else if (realAkk < static_cast<RealScalar>(0))
+          sign = NegativeSemiDef;
+      }
+    }
+
+    return ret;
+  }
+
+  // Reference for the algorithm: Davis and Hager, "Multiple Rank
+  // Modifications of a Sparse Cholesky Factorization" (Algorithm 1)
+  // Trivial rearrangements of their computations (Timothy E. Holy)
+  // allow their algorithm to work for rank-1 updates even if the
+  // original matrix is not of full rank.
+  // Here only rank-1 updates are implemented, to reduce the
+  // requirement for intermediate storage and improve accuracy
+  template <typename MatrixType, typename WDerived>
+  static bool updateInPlace(MatrixType& mat, MatrixBase<WDerived>& w,
+                            const typename MatrixType::RealScalar& sigma = 1) {
+    using numext::isfinite;
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename MatrixType::RealScalar RealScalar;
+
+    const Index size = mat.rows();
+    eigen_assert(mat.cols() == size && w.size() == size);
+
+    RealScalar alpha = 1;
+
+    // Apply the update
+    for (Index j = 0; j < size; j++) {
+      // Check for termination due to an original decomposition of low-rank
+      if (!(isfinite)(alpha)) break;
+
+      // Update the diagonal terms
+      RealScalar dj = numext::real(mat.coeff(j, j));
+      Scalar wj = w.coeff(j);
+      RealScalar swj2 = sigma * numext::abs2(wj);
+      RealScalar gamma = dj * alpha + swj2;
+
+      mat.coeffRef(j, j) += swj2 / alpha;
+      alpha += swj2 / dj;
+
+      // Update the terms of L
+      Index rs = size - j - 1;
+      w.tail(rs) -= wj * mat.col(j).tail(rs);
+      if (!numext::is_exactly_zero(gamma)) mat.col(j).tail(rs) += (sigma * numext::conj(wj) / gamma) * w.tail(rs);
+    }
+    return true;
+  }
+
+  template <typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
+  static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w,
+                     const typename MatrixType::RealScalar& sigma = 1) {
+    // Apply the permutation to the input w
+    tmp = transpositions * w;
+
+    return ldlt_inplace<Lower>::updateInPlace(mat, tmp, sigma);
+  }
+};
+
+template <>
+struct ldlt_inplace<Upper> {
+  template <typename MatrixType, typename TranspositionType, typename Workspace>
+  static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp,
+                                            SignMatrix& sign) {
+    Transpose<MatrixType> matt(mat);
+    return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
+  }
+
+  template <typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
+  static EIGEN_STRONG_INLINE bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& tmp, WType& w,
+                                         const typename MatrixType::RealScalar& sigma = 1) {
+    Transpose<MatrixType> matt(mat);
+    return ldlt_inplace<Lower>::update(matt, transpositions, tmp, w.conjugate(), sigma);
+  }
+};
+
+template <typename MatrixType>
+struct LDLT_Traits<MatrixType, Lower> {
+  typedef const TriangularView<const MatrixType, UnitLower> MatrixL;
+  typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitUpper> MatrixU;
+  static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); }
+  static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); }
+};
+
+template <typename MatrixType>
+struct LDLT_Traits<MatrixType, Upper> {
+  typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitLower> MatrixL;
+  typedef const TriangularView<const MatrixType, UnitUpper> MatrixU;
+  static inline MatrixL getL(const MatrixType& m) { return MatrixL(m.adjoint()); }
+  static inline MatrixU getU(const MatrixType& m) { return MatrixU(m); }
+};
+
+}  // end namespace internal
+
+/** Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of \a matrix
+ */
+template <typename MatrixType, int UpLo_>
+template <typename InputType>
+LDLT<MatrixType, UpLo_>& LDLT<MatrixType, UpLo_>::compute(const EigenBase<InputType>& a) {
+  eigen_assert(a.rows() == a.cols());
+  const Index size = a.rows();
+
+  m_matrix = a.derived();
+
+  // Compute matrix L1 norm = max abs column sum.
+  m_l1_norm = RealScalar(0);
+  // TODO move this code to SelfAdjointView
+  for (Index col = 0; col < size; ++col) {
+    RealScalar abs_col_sum;
+    if (UpLo_ == Lower)
+      abs_col_sum =
+          m_matrix.col(col).tail(size - col).template lpNorm<1>() + m_matrix.row(col).head(col).template lpNorm<1>();
+    else
+      abs_col_sum =
+          m_matrix.col(col).head(col).template lpNorm<1>() + m_matrix.row(col).tail(size - col).template lpNorm<1>();
+    if (abs_col_sum > m_l1_norm) m_l1_norm = abs_col_sum;
+  }
+
+  m_transpositions.resize(size);
+  m_isInitialized = false;
+  m_temporary.resize(size);
+  m_sign = internal::ZeroSign;
+
+  m_info = internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign) ? Success
+                                                                                                    : NumericalIssue;
+
+  m_isInitialized = true;
+  return *this;
+}
+
+/** Update the LDLT decomposition:  given A = L D L^T, efficiently compute the decomposition of A + sigma w w^T.
+ * \param w a vector to be incorporated into the decomposition.
+ * \param sigma a scalar, +1 for updates and -1 for "downdates," which correspond to removing previously-added column
+ * vectors. Optional; default value is +1. \sa setZero()
+ */
+template <typename MatrixType, int UpLo_>
+template <typename Derived>
+LDLT<MatrixType, UpLo_>& LDLT<MatrixType, UpLo_>::rankUpdate(
+    const MatrixBase<Derived>& w, const typename LDLT<MatrixType, UpLo_>::RealScalar& sigma) {
+  typedef typename TranspositionType::StorageIndex IndexType;
+  const Index size = w.rows();
+  if (m_isInitialized) {
+    eigen_assert(m_matrix.rows() == size);
+  } else {
+    m_matrix.resize(size, size);
+    m_matrix.setZero();
+    m_transpositions.resize(size);
+    for (Index i = 0; i < size; i++) m_transpositions.coeffRef(i) = IndexType(i);
+    m_temporary.resize(size);
+    m_sign = sigma >= 0 ? internal::PositiveSemiDef : internal::NegativeSemiDef;
+    m_isInitialized = true;
+  }
+
+  internal::ldlt_inplace<UpLo>::update(m_matrix, m_transpositions, m_temporary, w, sigma);
+
+  return *this;
+}
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+template <typename MatrixType_, int UpLo_>
+template <typename RhsType, typename DstType>
+void LDLT<MatrixType_, UpLo_>::_solve_impl(const RhsType& rhs, DstType& dst) const {
+  _solve_impl_transposed<true>(rhs, dst);
+}
+
+template <typename MatrixType_, int UpLo_>
+template <bool Conjugate, typename RhsType, typename DstType>
+void LDLT<MatrixType_, UpLo_>::_solve_impl_transposed(const RhsType& rhs, DstType& dst) const {
+  // dst = P b
+  dst = m_transpositions * rhs;
+
+  // dst = L^-1 (P b)
+  // dst = L^-*T (P b)
+  matrixL().template conjugateIf<!Conjugate>().solveInPlace(dst);
+
+  // dst = D^-* (L^-1 P b)
+  // dst = D^-1 (L^-*T P b)
+  // more precisely, use pseudo-inverse of D (see bug 241)
+  using std::abs;
+  const typename Diagonal<const MatrixType>::RealReturnType vecD(vectorD());
+  // In some previous versions, tolerance was set to the max of 1/highest (or rather numeric_limits::min())
+  // and the maximal diagonal entry * epsilon as motivated by LAPACK's xGELSS:
+  // RealScalar tolerance = numext::maxi(vecD.array().abs().maxCoeff() * NumTraits<RealScalar>::epsilon(),RealScalar(1)
+  // / NumTraits<RealScalar>::highest()); However, LDLT is not rank revealing, and so adjusting the tolerance wrt to the
+  // highest diagonal element is not well justified and leads to numerical issues in some cases. Moreover, Lapack's
+  // xSYTRS routines use 0 for the tolerance. Using numeric_limits::min() gives us more robustness to denormals.
+  RealScalar tolerance = (std::numeric_limits<RealScalar>::min)();
+  for (Index i = 0; i < vecD.size(); ++i) {
+    if (abs(vecD(i)) > tolerance)
+      dst.row(i) /= vecD(i);
+    else
+      dst.row(i).setZero();
+  }
+
+  // dst = L^-* (D^-* L^-1 P b)
+  // dst = L^-T (D^-1 L^-*T P b)
+  matrixL().transpose().template conjugateIf<Conjugate>().solveInPlace(dst);
+
+  // dst = P^T (L^-* D^-* L^-1 P b) = A^-1 b
+  // dst = P^-T (L^-T D^-1 L^-*T P b) = A^-1 b
+  dst = m_transpositions.transpose() * dst;
+}
+#endif
+
+/** \internal use x = ldlt_object.solve(x);
+ *
+ * This is the \em in-place version of solve().
+ *
+ * \param bAndX represents both the right-hand side matrix b and result x.
+ *
+ * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
+ *
+ * This version avoids a copy when the right hand side matrix b is not
+ * needed anymore.
+ *
+ * \sa LDLT::solve(), MatrixBase::ldlt()
+ */
+template <typename MatrixType, int UpLo_>
+template <typename Derived>
+bool LDLT<MatrixType, UpLo_>::solveInPlace(MatrixBase<Derived>& bAndX) const {
+  eigen_assert(m_isInitialized && "LDLT is not initialized.");
+  eigen_assert(m_matrix.rows() == bAndX.rows());
+
+  bAndX = this->solve(bAndX);
+
+  return true;
+}
+
+/** \returns the matrix represented by the decomposition,
+ * i.e., it returns the product: P^T L D L^* P.
+ * This function is provided for debug purpose. */
+template <typename MatrixType, int UpLo_>
+MatrixType LDLT<MatrixType, UpLo_>::reconstructedMatrix() const {
+  eigen_assert(m_isInitialized && "LDLT is not initialized.");
+  const Index size = m_matrix.rows();
+  MatrixType res(size, size);
+
+  // P
+  res.setIdentity();
+  res = transpositionsP() * res;
+  // L^* P
+  res = matrixU() * res;
+  // D(L^*P)
+  res = vectorD().real().asDiagonal() * res;
+  // L(DL^*P)
+  res = matrixL() * res;
+  // P^T (LDL^*P)
+  res = transpositionsP().transpose() * res;
+
+  return res;
+}
+
+/** \cholesky_module
+ * \returns the Cholesky decomposition with full pivoting without square root of \c *this
+ * \sa MatrixBase::ldlt()
+ */
+template <typename MatrixType, unsigned int UpLo>
+inline const LDLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
+SelfAdjointView<MatrixType, UpLo>::ldlt() const {
+  return LDLT<PlainObject, UpLo>(m_matrix);
+}
+
+/** \cholesky_module
+ * \returns the Cholesky decomposition with full pivoting without square root of \c *this
+ * \sa SelfAdjointView::ldlt()
+ */
+template <typename Derived>
+inline const LDLT<typename MatrixBase<Derived>::PlainObject> MatrixBase<Derived>::ldlt() const {
+  return LDLT<PlainObject>(derived());
+}
+
+}  // end namespace Eigen
+
+#endif  // EIGEN_LDLT_H

eigen-master/Eigen/src/Cholesky/LLT.h

+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_LLT_H
+#define EIGEN_LLT_H
+
+// IWYU pragma: private
+#include "./InternalHeaderCheck.h"
+
+namespace Eigen {
+
+namespace internal {
+
+template <typename MatrixType_, int UpLo_>
+struct traits<LLT<MatrixType_, UpLo_> > : traits<MatrixType_> {
+  typedef MatrixXpr XprKind;
+  typedef SolverStorage StorageKind;
+  typedef int StorageIndex;
+  enum { Flags = 0 };
+};
+
+template <typename MatrixType, int UpLo>
+struct LLT_Traits;
+}  // namespace internal
+
+/** \ingroup Cholesky_Module
+ *
+ * \class LLT
+ *
+ * \brief Standard Cholesky decomposition (LL^T) of a matrix and associated features
+ *
+ * \tparam MatrixType_ the type of the matrix of which we are computing the LL^T Cholesky decomposition
+ * \tparam UpLo_ the triangular part that will be used for the decomposition: Lower (default) or Upper.
+ *               The other triangular part won't be read.
+ *
+ * This class performs a LL^T Cholesky decomposition of a symmetric, positive definite
+ * matrix A such that A = LL^* = U^*U, where L is lower triangular.
+ *
+ * While the Cholesky decomposition is particularly useful to solve selfadjoint problems like  D^*D x = b,
+ * for that purpose, we recommend the Cholesky decomposition without square root which is more stable
+ * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other
+ * situations like generalised eigen problems with hermitian matrices.
+ *
+ * Remember that Cholesky decompositions are not rank-revealing. This LLT decomposition is only stable on positive
+ * definite matrices, use LDLT instead for the semidefinite case. Also, do not use a Cholesky decomposition to determine
+ * whether a system of equations has a solution.
+ *
+ * Example: \include LLT_example.cpp
+ * Output: \verbinclude LLT_example.out
+ *
+ * \b Performance: for best performance, it is recommended to use a column-major storage format
+ * with the Lower triangular part (the default), or, equivalently, a row-major storage format
+ * with the Upper triangular part. Otherwise, you might get a 20% slowdown for the full factorization
+ * step, and rank-updates can be up to 3 times slower.
+ *
+ * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
+ *
+ * Note that during the decomposition, only the lower (or upper, as defined by UpLo_) triangular part of A is
+ * considered. Therefore, the strict lower part does not have to store correct values.
+ *
+ * \sa MatrixBase::llt(), SelfAdjointView::llt(), class LDLT
+ */
+template <typename MatrixType_, int UpLo_>
+class LLT : public SolverBase<LLT<MatrixType_, UpLo_> > {
+ public:
+  typedef MatrixType_ MatrixType;
+  typedef SolverBase<LLT> Base;
+  friend class SolverBase<LLT>;
+
+  EIGEN_GENERIC_PUBLIC_INTERFACE(LLT)
+  enum { MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime };
+
+  enum { PacketSize = internal::packet_traits<Scalar>::size, AlignmentMask = int(PacketSize) - 1, UpLo = UpLo_ };
+
+  typedef internal::LLT_Traits<MatrixType, UpLo> Traits;
+
+  /**
+   * \brief Default Constructor.
+   *
+   * The default constructor is useful in cases in which the user intends to
+   * perform decompositions via LLT::compute(const MatrixType&).
+   */
+  LLT() : m_matrix(), m_isInitialized(false) {}
+
+  /** \brief Default Constructor with memory preallocation
+   *
+   * Like the default constructor but with preallocation of the internal data
+   * according to the specified problem \a size.
+   * \sa LLT()
+   */
+  explicit LLT(Index size) : m_matrix(size, size), m_isInitialized(false) {}
+
+  template <typename InputType>
+  explicit LLT(const EigenBase<InputType>& matrix) : m_matrix(matrix.rows(), matrix.cols()), m_isInitialized(false) {
+    compute(matrix.derived());
+  }
+
+  /** \brief Constructs a LLT factorization from a given matrix
+   *
+   * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when
+   * \c MatrixType is a Eigen::Ref.
+   *
+   * \sa LLT(const EigenBase&)
+   */
+  template <typename InputType>
+  explicit LLT(EigenBase<InputType>& matrix) : m_matrix(matrix.derived()), m_isInitialized(false) {
+    compute(matrix.derived());
+  }
+
+  /** \returns a view of the upper triangular matrix U */
+  inline typename Traits::MatrixU matrixU() const {
+    eigen_assert(m_isInitialized && "LLT is not initialized.");
+    return Traits::getU(m_matrix);
+  }
+
+  /** \returns a view of the lower triangular matrix L */
+  inline typename Traits::MatrixL matrixL() const {
+    eigen_assert(m_isInitialized && "LLT is not initialized.");
+    return Traits::getL(m_matrix);
+  }
+
+#ifdef EIGEN_PARSED_BY_DOXYGEN
+  /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
+   *
+   * Since this LLT class assumes anyway that the matrix A is invertible, the solution
+   * theoretically exists and is unique regardless of b.
+   *
+   * Example: \include LLT_solve.cpp
+   * Output: \verbinclude LLT_solve.out
+   *
+   * \sa solveInPlace(), MatrixBase::llt(), SelfAdjointView::llt()
+   */
+  template <typename Rhs>
+  inline const Solve<LLT, Rhs> solve(const MatrixBase<Rhs>& b) const;
+#endif
+
+  template <typename Derived>
+  void solveInPlace(const MatrixBase<Derived>& bAndX) const;
+
+  template <typename InputType>
+  LLT& compute(const EigenBase<InputType>& matrix);
+
+  /** \returns an estimate of the reciprocal condition number of the matrix of
+   *  which \c *this is the Cholesky decomposition.
+   */
+  RealScalar rcond() const {
+    eigen_assert(m_isInitialized && "LLT is not initialized.");
+    eigen_assert(m_info == Success && "LLT failed because matrix appears to be negative");
+    return internal::rcond_estimate_helper(m_l1_norm, *this);
+  }
+
+  /** \returns the LLT decomposition matrix
+   *
+   * TODO: document the storage layout
+   */
+  inline const MatrixType& matrixLLT() const {
+    eigen_assert(m_isInitialized && "LLT is not initialized.");
+    return m_matrix;
+  }
+
+  MatrixType reconstructedMatrix() const;
+
+  /** \brief Reports whether previous computation was successful.
+   *
+   * \returns \c Success if computation was successful,
+   *          \c NumericalIssue if the matrix.appears not to be positive definite.
+   */
+  ComputationInfo info() const {
+    eigen_assert(m_isInitialized && "LLT is not initialized.");
+    return m_info;
+  }
+
+  /** \returns the adjoint of \c *this, that is, a const reference to the decomposition itself as the underlying matrix
+   * is self-adjoint.
+   *
+   * This method is provided for compatibility with other matrix decompositions, thus enabling generic code such as:
+   * \code x = decomposition.adjoint().solve(b) \endcode
+   */
+  const LLT& adjoint() const noexcept { return *this; }
+
+  constexpr Index rows() const noexcept { return m_matrix.rows(); }
+  constexpr Index cols() const noexcept { return m_matrix.cols(); }
+
+  template <typename VectorType>
+  LLT& rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+  template <typename RhsType, typename DstType>
+  void _solve_impl(const RhsType& rhs, DstType& dst) const;
+
+  template <bool Conjugate, typename RhsType, typename DstType>
+  void _solve_impl_transposed(const RhsType& rhs, DstType& dst) const;
+#endif
+
+ protected:
+  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
+
+  /** \internal
+   * Used to compute and store L
+   * The strict upper part is not used and even not initialized.
+   */
+  MatrixType m_matrix;
+  RealScalar m_l1_norm;
+  bool m_isInitialized;
+  ComputationInfo m_info;
+};
+
+namespace internal {
+
+template <typename Scalar, int UpLo>
+struct llt_inplace;
+
+template <typename MatrixType, typename VectorType>
+static Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec,
+                                   const typename MatrixType::RealScalar& sigma) {
+  using std::sqrt;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  typedef typename MatrixType::ColXpr ColXpr;
+  typedef internal::remove_all_t<ColXpr> ColXprCleaned;
+  typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
+  typedef Matrix<Scalar, Dynamic, 1> TempVectorType;
+  typedef typename TempVectorType::SegmentReturnType TempVecSegment;
+
+  Index n = mat.cols();
+  eigen_assert(mat.rows() == n && vec.size() == n);
+
+  TempVectorType temp;
+
+  if (sigma > 0) {
+    // This version is based on Givens rotations.
+    // It is faster than the other one below, but only works for updates,
+    // i.e., for sigma > 0
+    temp = sqrt(sigma) * vec;
+
+    for (Index i = 0; i < n; ++i) {
+      JacobiRotation<Scalar> g;
+      g.makeGivens(mat(i, i), -temp(i), &mat(i, i));
+
+      Index rs = n - i - 1;
+      if (rs > 0) {
+        ColXprSegment x(mat.col(i).tail(rs));
+        TempVecSegment y(temp.tail(rs));
+        apply_rotation_in_the_plane(x, y, g);
+      }
+    }
+  } else {
+    temp = vec;
+    RealScalar beta = 1;
+    for (Index j = 0; j < n; ++j) {
+      RealScalar Ljj = numext::real(mat.coeff(j, j));
+      RealScalar dj = numext::abs2(Ljj);
+      Scalar wj = temp.coeff(j);
+      RealScalar swj2 = sigma * numext::abs2(wj);
+      RealScalar gamma = dj * beta + swj2;
+
+      RealScalar x = dj + swj2 / beta;
+      if (x <= RealScalar(0)) return j;
+      RealScalar nLjj = sqrt(x);
+      mat.coeffRef(j, j) = nLjj;
+      beta += swj2 / dj;
+
+      // Update the terms of L
+      Index rs = n - j - 1;
+      if (rs) {
+        temp.tail(rs) -= (wj / Ljj) * mat.col(j).tail(rs);
+        if (!numext::is_exactly_zero(gamma))
+          mat.col(j).tail(rs) =
+              (nLjj / Ljj) * mat.col(j).tail(rs) + (nLjj * sigma * numext::conj(wj) / gamma) * temp.tail(rs);
+      }
+    }
+  }
+  return -1;
+}
+
+template <typename Scalar>
+struct llt_inplace<Scalar, Lower> {
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+  template <typename MatrixType>
+  static Index unblocked(MatrixType& mat) {
+    using std::sqrt;
+
+    eigen_assert(mat.rows() == mat.cols());
+    const Index size = mat.rows();
+    for (Index k = 0; k < size; ++k) {
+      Index rs = size - k - 1;  // remaining size
+
+      Block<MatrixType, Dynamic, 1> A21(mat, k + 1, k, rs, 1);
+      Block<MatrixType, 1, Dynamic> A10(mat, k, 0, 1, k);
+      Block<MatrixType, Dynamic, Dynamic> A20(mat, k + 1, 0, rs, k);
+
+      RealScalar x = numext::real(mat.coeff(k, k));
+      if (k > 0) x -= A10.squaredNorm();
+      if (x <= RealScalar(0)) return k;
+      mat.coeffRef(k, k) = x = sqrt(x);
+      if (k > 0 && rs > 0) A21.noalias() -= A20 * A10.adjoint();
+      if (rs > 0) A21 /= x;
+    }
+    return -1;
+  }
+
+  template <typename MatrixType>
+  static Index blocked(MatrixType& m) {
+    eigen_assert(m.rows() == m.cols());
+    Index size = m.rows();
+    if (size < 32) return unblocked(m);
+
+    Index blockSize = size / 8;
+    blockSize = (blockSize / 16) * 16;
+    blockSize = (std::min)((std::max)(blockSize, Index(8)), Index(128));
+
+    for (Index k = 0; k < size; k += blockSize) {
+      // partition the matrix:
+      //       A00 |  -  |  -
+      // lu  = A10 | A11 |  -
+      //       A20 | A21 | A22
+      Index bs = (std::min)(blockSize, size - k);
+      Index rs = size - k - bs;
+      Block<MatrixType, Dynamic, Dynamic> A11(m, k, k, bs, bs);
+      Block<MatrixType, Dynamic, Dynamic> A21(m, k + bs, k, rs, bs);
+      Block<MatrixType, Dynamic, Dynamic> A22(m, k + bs, k + bs, rs, rs);
+
+      Index ret;
+      if ((ret = unblocked(A11)) >= 0) return k + ret;
+      if (rs > 0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
+      if (rs > 0)
+        A22.template selfadjointView<Lower>().rankUpdate(A21,
+                                                         typename NumTraits<RealScalar>::Literal(-1));  // bottleneck
+    }
+    return -1;
+  }
+
+  template <typename MatrixType, typename VectorType>
+  static Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma) {
+    return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
+  }
+};
+
+template <typename Scalar>
+struct llt_inplace<Scalar, Upper> {
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+
+  template <typename MatrixType>
+  static EIGEN_STRONG_INLINE Index unblocked(MatrixType& mat) {
+    Transpose<MatrixType> matt(mat);
+    return llt_inplace<Scalar, Lower>::unblocked(matt);
+  }
+  template <typename MatrixType>
+  static EIGEN_STRONG_INLINE Index blocked(MatrixType& mat) {
+    Transpose<MatrixType> matt(mat);
+    return llt_inplace<Scalar, Lower>::blocked(matt);
+  }
+  template <typename MatrixType, typename VectorType>
+  static Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma) {
+    Transpose<MatrixType> matt(mat);
+    return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma);
+  }
+};
+
+template <typename MatrixType>
+struct LLT_Traits<MatrixType, Lower> {
+  typedef const TriangularView<const MatrixType, Lower> MatrixL;
+  typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU;
+  static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); }
+  static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); }
+  static bool inplace_decomposition(MatrixType& m) {
+    return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m) == -1;
+  }
+};
+
+template <typename MatrixType>
+struct LLT_Traits<MatrixType, Upper> {
+  typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL;
+  typedef const TriangularView<const MatrixType, Upper> MatrixU;
+  static inline MatrixL getL(const MatrixType& m) { return MatrixL(m.adjoint()); }
+  static inline MatrixU getU(const MatrixType& m) { return MatrixU(m); }
+  static bool inplace_decomposition(MatrixType& m) {
+    return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m) == -1;
+  }
+};
+
+}  // end namespace internal
+
+/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
+ *
+ * \returns a reference to *this
+ *
+ * Example: \include TutorialLinAlgComputeTwice.cpp
+ * Output: \verbinclude TutorialLinAlgComputeTwice.out
+ */
+template <typename MatrixType, int UpLo_>
+template <typename InputType>
+LLT<MatrixType, UpLo_>& LLT<MatrixType, UpLo_>::compute(const EigenBase<InputType>& a) {
+  eigen_assert(a.rows() == a.cols());
+  const Index size = a.rows();
+  m_matrix.resize(size, size);
+  if (!internal::is_same_dense(m_matrix, a.derived())) m_matrix = a.derived();
+
+  // Compute matrix L1 norm = max abs column sum.
+  m_l1_norm = RealScalar(0);
+  // TODO move this code to SelfAdjointView
+  for (Index col = 0; col < size; ++col) {
+    RealScalar abs_col_sum;
+    if (UpLo_ == Lower)
+      abs_col_sum =
+          m_matrix.col(col).tail(size - col).template lpNorm<1>() + m_matrix.row(col).head(col).template lpNorm<1>();
+    else
+      abs_col_sum =
+          m_matrix.col(col).head(col).template lpNorm<1>() + m_matrix.row(col).tail(size - col).template lpNorm<1>();
+    if (abs_col_sum > m_l1_norm) m_l1_norm = abs_col_sum;
+  }
+
+  m_isInitialized = true;
+  bool ok = Traits::inplace_decomposition(m_matrix);
+  m_info = ok ? Success : NumericalIssue;
+
+  return *this;
+}
+
+/** Performs a rank one update (or dowdate) of the current decomposition.
+ * If A = LL^* before the rank one update,
+ * then after it we have LL^* = A + sigma * v v^* where \a v must be a vector
+ * of same dimension.
+ */
+template <typename MatrixType_, int UpLo_>
+template <typename VectorType>
+LLT<MatrixType_, UpLo_>& LLT<MatrixType_, UpLo_>::rankUpdate(const VectorType& v, const RealScalar& sigma) {
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
+  eigen_assert(v.size() == m_matrix.cols());
+  eigen_assert(m_isInitialized);
+  if (internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix, v, sigma) >= 0)
+    m_info = NumericalIssue;
+  else
+    m_info = Success;
+
+  return *this;
+}
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+template <typename MatrixType_, int UpLo_>
+template <typename RhsType, typename DstType>
+void LLT<MatrixType_, UpLo_>::_solve_impl(const RhsType& rhs, DstType& dst) const {
+  _solve_impl_transposed<true>(rhs, dst);
+}
+
+template <typename MatrixType_, int UpLo_>
+template <bool Conjugate, typename RhsType, typename DstType>
+void LLT<MatrixType_, UpLo_>::_solve_impl_transposed(const RhsType& rhs, DstType& dst) const {
+  dst = rhs;
+
+  matrixL().template conjugateIf<!Conjugate>().solveInPlace(dst);
+  matrixU().template conjugateIf<!Conjugate>().solveInPlace(dst);
+}
+#endif
+
+/** \internal use x = llt_object.solve(x);
+ *
+ * This is the \em in-place version of solve().
+ *
+ * \param bAndX represents both the right-hand side matrix b and result x.
+ *
+ * This version avoids a copy when the right hand side matrix b is not needed anymore.
+ *
+ * \warning The parameter is only marked 'const' to make the C++ compiler accept a temporary expression here.
+ * This function will const_cast it, so constness isn't honored here.
+ *
+ * \sa LLT::solve(), MatrixBase::llt()
+ */
+template <typename MatrixType, int UpLo_>
+template <typename Derived>
+void LLT<MatrixType, UpLo_>::solveInPlace(const MatrixBase<Derived>& bAndX) const {
+  eigen_assert(m_isInitialized && "LLT is not initialized.");
+  eigen_assert(m_matrix.rows() == bAndX.rows());
+  matrixL().solveInPlace(bAndX);
+  matrixU().solveInPlace(bAndX);
+}
+
+/** \returns the matrix represented by the decomposition,
+ * i.e., it returns the product: L L^*.
+ * This function is provided for debug purpose. */
+template <typename MatrixType, int UpLo_>
+MatrixType LLT<MatrixType, UpLo_>::reconstructedMatrix() const {
+  eigen_assert(m_isInitialized && "LLT is not initialized.");
+  return matrixL() * matrixL().adjoint().toDenseMatrix();
+}
+
+/** \cholesky_module
+ * \returns the LLT decomposition of \c *this
+ * \sa SelfAdjointView::llt()
+ */
+template <typename Derived>
+inline const LLT<typename MatrixBase<Derived>::PlainObject> MatrixBase<Derived>::llt() const {
+  return LLT<PlainObject>(derived());
+}
+
+/** \cholesky_module
+ * \returns the LLT decomposition of \c *this
+ * \sa SelfAdjointView::llt()
+ */
+template <typename MatrixType, unsigned int UpLo>
+inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo> SelfAdjointView<MatrixType, UpLo>::llt()
+    const {
+  return LLT<PlainObject, UpLo>(m_matrix);
+}
+
+}  // end namespace Eigen
+
+#endif  // EIGEN_LLT_H

eigen-master/Eigen/src/Cholesky/LLT_LAPACKE.h

+/*
+ Copyright (c) 2011, Intel Corporation. All rights reserved.
+
+ Redistribution and use in source and binary forms, with or without modification,
+ are permitted provided that the following conditions are met:
+
+ * Redistributions of source code must retain the above copyright notice, this
+   list of conditions and the following disclaimer.
+ * Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+ * Neither the name of Intel Corporation nor the names of its contributors may
+   be used to endorse or promote products derived from this software without
+   specific prior written permission.
+
+ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+ ********************************************************************************
+ *   Content : Eigen bindings to LAPACKe
+ *     LLt decomposition based on LAPACKE_?potrf function.
+ ********************************************************************************
+*/
+
+#ifndef EIGEN_LLT_LAPACKE_H
+#define EIGEN_LLT_LAPACKE_H
+
+// IWYU pragma: private
+#include "./InternalHeaderCheck.h"
+
+namespace Eigen {
+
+namespace internal {
+
+namespace lapacke_helpers {
+// -------------------------------------------------------------------------------------------------------------------
+//        Dispatch for rank update handling upper and lower parts
+// -------------------------------------------------------------------------------------------------------------------
+
+template <UpLoType Mode>
+struct rank_update {};
+
+template <>
+struct rank_update<Lower> {
+  template <typename MatrixType, typename VectorType>
+  static Index run(MatrixType &mat, const VectorType &vec, const typename MatrixType::RealScalar &sigma) {
+    return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
+  }
+};
+
+template <>
+struct rank_update<Upper> {
+  template <typename MatrixType, typename VectorType>
+  static Index run(MatrixType &mat, const VectorType &vec, const typename MatrixType::RealScalar &sigma) {
+    Transpose<MatrixType> matt(mat);
+    return Eigen::internal::llt_rank_update_lower(matt, vec.conjugate(), sigma);
+  }
+};
+
+// -------------------------------------------------------------------------------------------------------------------
+//        Generic lapacke llt implementation that hands of to the dispatches
+// -------------------------------------------------------------------------------------------------------------------
+
+template <typename Scalar, UpLoType Mode>
+struct lapacke_llt {
+  EIGEN_STATIC_ASSERT(((Mode == Lower) || (Mode == Upper)), MODE_MUST_BE_UPPER_OR_LOWER)
+  template <typename MatrixType>
+  static Index blocked(MatrixType &m) {
+    eigen_assert(m.rows() == m.cols());
+    if (m.rows() == 0) {
+      return -1;
+    }
+    /* Set up parameters for ?potrf */
+    lapack_int size = to_lapack(m.rows());
+    lapack_int matrix_order = lapack_storage_of(m);
+    constexpr char uplo = Mode == Upper ? 'U' : 'L';
+    Scalar *a = &(m.coeffRef(0, 0));
+    lapack_int lda = to_lapack(m.outerStride());
+
+    lapack_int info = potrf(matrix_order, uplo, size, to_lapack(a), lda);
+    info = (info == 0) ? -1 : info > 0 ? info - 1 : size;
+    return info;
+  }
+
+  template <typename MatrixType, typename VectorType>
+  static Index rankUpdate(MatrixType &mat, const VectorType &vec, const typename MatrixType::RealScalar &sigma) {
+    return rank_update<Mode>::run(mat, vec, sigma);
+  }
+};
+}  // namespace lapacke_helpers
+// end namespace lapacke_helpers
+
+/*
+ * Here, we just put the generic implementation from lapacke_llt into a full specialization of the llt_inplace
+ * type. By being a full specialization, the versions defined here thus get precedence over the generic implementation
+ * in LLT.h for double, float and complex double, complex float types.
+ */
+
+#define EIGEN_LAPACKE_LLT(EIGTYPE)                                                             \
+  template <>                                                                                  \
+  struct llt_inplace<EIGTYPE, Lower> : public lapacke_helpers::lapacke_llt<EIGTYPE, Lower> {}; \
+  template <>                                                                                  \
+  struct llt_inplace<EIGTYPE, Upper> : public lapacke_helpers::lapacke_llt<EIGTYPE, Upper> {};
+
+EIGEN_LAPACKE_LLT(double)
+EIGEN_LAPACKE_LLT(float)
+EIGEN_LAPACKE_LLT(std::complex<double>)
+EIGEN_LAPACKE_LLT(std::complex<float>)
+
+#undef EIGEN_LAPACKE_LLT
+
+}  // end namespace internal
+
+}  // end namespace Eigen
+
+#endif  // EIGEN_LLT_LAPACKE_H

eigen-master/Eigen/src/CholmodSupport/CholmodSupport.h

+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_CHOLMODSUPPORT_H
+#define EIGEN_CHOLMODSUPPORT_H
+
+// IWYU pragma: private
+#include "./InternalHeaderCheck.h"
+
+namespace Eigen {
+
+namespace internal {
+
+template <typename Scalar>
+struct cholmod_configure_matrix;
+
+template <>
+struct cholmod_configure_matrix<double> {
+  template <typename CholmodType>
+  static void run(CholmodType& mat) {
+    mat.xtype = CHOLMOD_REAL;
+    mat.dtype = CHOLMOD_DOUBLE;
+  }
+};
+
+template <>
+struct cholmod_configure_matrix<std::complex<double> > {
+  template <typename CholmodType>
+  static void run(CholmodType& mat) {
+    mat.xtype = CHOLMOD_COMPLEX;
+    mat.dtype = CHOLMOD_DOUBLE;
+  }
+};
+
+// Other scalar types are not yet supported by Cholmod
+// template<> struct cholmod_configure_matrix<float> {
+//   template<typename CholmodType>
+//   static void run(CholmodType& mat) {
+//     mat.xtype = CHOLMOD_REAL;
+//     mat.dtype = CHOLMOD_SINGLE;
+//   }
+// };
+//
+// template<> struct cholmod_configure_matrix<std::complex<float> > {
+//   template<typename CholmodType>
+//   static void run(CholmodType& mat) {
+//     mat.xtype = CHOLMOD_COMPLEX;
+//     mat.dtype = CHOLMOD_SINGLE;
+//   }
+// };
+
+}  // namespace internal
+
+/** Wraps the Eigen sparse matrix \a mat into a Cholmod sparse matrix object.
+ * Note that the data are shared.
+ */
+template <typename Scalar_, int Options_, typename StorageIndex_>
+cholmod_sparse viewAsCholmod(Ref<SparseMatrix<Scalar_, Options_, StorageIndex_> > mat) {
+  cholmod_sparse res;
+  res.nzmax = mat.nonZeros();
+  res.nrow = mat.rows();
+  res.ncol = mat.cols();
+  res.p = mat.outerIndexPtr();
+  res.i = mat.innerIndexPtr();
+  res.x = mat.valuePtr();
+  res.z = 0;
+  res.sorted = 1;
+  if (mat.isCompressed()) {
+    res.packed = 1;
+    res.nz = 0;
+  } else {
+    res.packed = 0;
+    res.nz = mat.innerNonZeroPtr();
+  }
+
+  res.dtype = 0;
+  res.stype = -1;
+
+  if (internal::is_same<StorageIndex_, int>::value) {
+    res.itype = CHOLMOD_INT;
+  } else if (internal::is_same<StorageIndex_, SuiteSparse_long>::value) {
+    res.itype = CHOLMOD_LONG;
+  } else {
+    eigen_assert(false && "Index type not supported yet");
+  }
+
+  // setup res.xtype
+  internal::cholmod_configure_matrix<Scalar_>::run(res);
+
+  res.stype = 0;
+
+  return res;
+}
+
+template <typename Scalar_, int Options_, typename Index_>
+const cholmod_sparse viewAsCholmod(const SparseMatrix<Scalar_, Options_, Index_>& mat) {
+  cholmod_sparse res = viewAsCholmod(Ref<SparseMatrix<Scalar_, Options_, Index_> >(mat.const_cast_derived()));
+  return res;
+}
+
+template <typename Scalar_, int Options_, typename Index_>
+const cholmod_sparse viewAsCholmod(const SparseVector<Scalar_, Options_, Index_>& mat) {
+  cholmod_sparse res = viewAsCholmod(Ref<SparseMatrix<Scalar_, Options_, Index_> >(mat.const_cast_derived()));
+  return res;
+}
+
+/** Returns a view of the Eigen sparse matrix \a mat as Cholmod sparse matrix.
+ * The data are not copied but shared. */
+template <typename Scalar_, int Options_, typename Index_, unsigned int UpLo>
+cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<const SparseMatrix<Scalar_, Options_, Index_>, UpLo>& mat) {
+  cholmod_sparse res = viewAsCholmod(Ref<SparseMatrix<Scalar_, Options_, Index_> >(mat.matrix().const_cast_derived()));
+
+  if (UpLo == Upper) res.stype = 1;
+  if (UpLo == Lower) res.stype = -1;
+  // swap stype for rowmajor matrices (only works for real matrices)
+  EIGEN_STATIC_ASSERT((Options_ & RowMajorBit) == 0 || NumTraits<Scalar_>::IsComplex == 0,
+                      THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
+  if (Options_ & RowMajorBit) res.stype *= -1;
+
+  return res;
+}
+
+/** Returns a view of the Eigen \b dense matrix \a mat as Cholmod dense matrix.
+ * The data are not copied but shared. */
+template <typename Derived>
+cholmod_dense viewAsCholmod(MatrixBase<Derived>& mat) {
+  EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags & RowMajorBit) == 0,
+                      THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
+  typedef typename Derived::Scalar Scalar;
+
+  cholmod_dense res;
+  res.nrow = mat.rows();
+  res.ncol = mat.cols();
+  res.nzmax = res.nrow * res.ncol;
+  res.d = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride();
+  res.x = (void*)(mat.derived().data());
+  res.z = 0;
+
+  internal::cholmod_configure_matrix<Scalar>::run(res);
+
+  return res;
+}
+
+/** Returns a view of the Cholmod sparse matrix \a cm as an Eigen sparse matrix.
+ * The data are not copied but shared. */
+template <typename Scalar, typename StorageIndex>
+Map<const SparseMatrix<Scalar, ColMajor, StorageIndex> > viewAsEigen(cholmod_sparse& cm) {
+  return Map<const SparseMatrix<Scalar, ColMajor, StorageIndex> >(
+      cm.nrow, cm.ncol, static_cast<StorageIndex*>(cm.p)[cm.ncol], static_cast<StorageIndex*>(cm.p),
+      static_cast<StorageIndex*>(cm.i), static_cast<Scalar*>(cm.x));
+}
+
+/** Returns a view of the Cholmod sparse matrix factor \a cm as an Eigen sparse matrix.
+ * The data are not copied but shared. */
+template <typename Scalar, typename StorageIndex>
+Map<const SparseMatrix<Scalar, ColMajor, StorageIndex> > viewAsEigen(cholmod_factor& cm) {
+  return Map<const SparseMatrix<Scalar, ColMajor, StorageIndex> >(
+      cm.n, cm.n, static_cast<StorageIndex*>(cm.p)[cm.n], static_cast<StorageIndex*>(cm.p),
+      static_cast<StorageIndex*>(cm.i), static_cast<Scalar*>(cm.x));
+}
+
+namespace internal {
+
+// template specializations for int and long that call the correct cholmod method
+
+#define EIGEN_CHOLMOD_SPECIALIZE0(ret, name)                        \
+  template <typename StorageIndex_>                                 \
+  inline ret cm_##name(cholmod_common& Common) {                    \
+    return cholmod_##name(&Common);                                 \
+  }                                                                 \
+  template <>                                                       \
+  inline ret cm_##name<SuiteSparse_long>(cholmod_common & Common) { \
+    return cholmod_l_##name(&Common);                               \
+  }
+
+#define EIGEN_CHOLMOD_SPECIALIZE1(ret, name, t1, a1)                         \
+  template <typename StorageIndex_>                                          \
+  inline ret cm_##name(t1& a1, cholmod_common& Common) {                     \
+    return cholmod_##name(&a1, &Common);                                     \
+  }                                                                          \
+  template <>                                                                \
+  inline ret cm_##name<SuiteSparse_long>(t1 & a1, cholmod_common & Common) { \
+    return cholmod_l_##name(&a1, &Common);                                   \
+  }
+
+EIGEN_CHOLMOD_SPECIALIZE0(int, start)
+EIGEN_CHOLMOD_SPECIALIZE0(int, finish)
+
+EIGEN_CHOLMOD_SPECIALIZE1(int, free_factor, cholmod_factor*, L)
+EIGEN_CHOLMOD_SPECIALIZE1(int, free_dense, cholmod_dense*, X)
+EIGEN_CHOLMOD_SPECIALIZE1(int, free_sparse, cholmod_sparse*, A)
+
+EIGEN_CHOLMOD_SPECIALIZE1(cholmod_factor*, analyze, cholmod_sparse, A)
+EIGEN_CHOLMOD_SPECIALIZE1(cholmod_sparse*, factor_to_sparse, cholmod_factor, L)
+
+template <typename StorageIndex_>
+inline cholmod_dense* cm_solve(int sys, cholmod_factor& L, cholmod_dense& B, cholmod_common& Common) {
+  return cholmod_solve(sys, &L, &B, &Common);
+}
+template <>
+inline cholmod_dense* cm_solve<SuiteSparse_long>(int sys, cholmod_factor& L, cholmod_dense& B, cholmod_common& Common) {
+  return cholmod_l_solve(sys, &L, &B, &Common);
+}
+
+template <typename StorageIndex_>
+inline cholmod_sparse* cm_spsolve(int sys, cholmod_factor& L, cholmod_sparse& B, cholmod_common& Common) {
+  return cholmod_spsolve(sys, &L, &B, &Common);
+}
+template <>
+inline cholmod_sparse* cm_spsolve<SuiteSparse_long>(int sys, cholmod_factor& L, cholmod_sparse& B,
+                                                    cholmod_common& Common) {
+  return cholmod_l_spsolve(sys, &L, &B, &Common);
+}
+
+template <typename StorageIndex_>
+inline int cm_factorize_p(cholmod_sparse* A, double beta[2], StorageIndex_* fset, std::size_t fsize, cholmod_factor* L,
+                          cholmod_common& Common) {
+  return cholmod_factorize_p(A, beta, fset, fsize, L, &Common);
+}
+template <>
+inline int cm_factorize_p<SuiteSparse_long>(cholmod_sparse* A, double beta[2], SuiteSparse_long* fset,
+                                            std::size_t fsize, cholmod_factor* L, cholmod_common& Common) {
+  return cholmod_l_factorize_p(A, beta, fset, fsize, L, &Common);
+}
+
+#undef EIGEN_CHOLMOD_SPECIALIZE0
+#undef EIGEN_CHOLMOD_SPECIALIZE1
+
+}  // namespace internal
+
+enum CholmodMode { CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt };
+
+/** \ingroup CholmodSupport_Module
+ * \class CholmodBase
+ * \brief The base class for the direct Cholesky factorization of Cholmod
+ * \sa class CholmodSupernodalLLT, class CholmodSimplicialLDLT, class CholmodSimplicialLLT
+ */
+template <typename MatrixType_, int UpLo_, typename Derived>
+class CholmodBase : public SparseSolverBase<Derived> {
+ protected:
+  typedef SparseSolverBase<Derived> Base;
+  using Base::derived;
+  using Base::m_isInitialized;
+
+ public:
+  typedef MatrixType_ MatrixType;
+  enum { UpLo = UpLo_ };
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  typedef MatrixType CholMatrixType;
+  typedef typename MatrixType::StorageIndex StorageIndex;
+  enum { ColsAtCompileTime = MatrixType::ColsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime };
+
+ public:
+  CholmodBase() : m_cholmodFactor(0), m_info(Success), m_factorizationIsOk(false), m_analysisIsOk(false) {
+    EIGEN_STATIC_ASSERT((internal::is_same<double, RealScalar>::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY);
+    m_shiftOffset[0] = m_shiftOffset[1] = 0.0;
+    internal::cm_start<StorageIndex>(m_cholmod);
+  }
+
+  explicit CholmodBase(const MatrixType& matrix)
+      : m_cholmodFactor(0), m_info(Success), m_factorizationIsOk(false), m_analysisIsOk(false) {
+    EIGEN_STATIC_ASSERT((internal::is_same<double, RealScalar>::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY);
+    m_shiftOffset[0] = m_shiftOffset[1] = 0.0;
+    internal::cm_start<StorageIndex>(m_cholmod);
+    compute(matrix);
+  }
+
+  ~CholmodBase() {
+    if (m_cholmodFactor) internal::cm_free_factor<StorageIndex>(m_cholmodFactor, m_cholmod);
+    internal::cm_finish<StorageIndex>(m_cholmod);
+  }
+
+  inline StorageIndex cols() const { return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }
+  inline StorageIndex rows() const { return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }
+
+  /** \brief Reports whether previous computation was successful.
+   *
+   * \returns \c Success if computation was successful,
+   *          \c NumericalIssue if the matrix.appears to be negative.
+   */
+  ComputationInfo info() const {
+    eigen_assert(m_isInitialized && "Decomposition is not initialized.");
+    return m_info;
+  }
+
+  /** Computes the sparse Cholesky decomposition of \a matrix */
+  Derived& compute(const MatrixType& matrix) {
+    analyzePattern(matrix);
+    factorize(matrix);
+    return derived();
+  }
+
+  /** Performs a symbolic decomposition on the sparsity pattern of \a matrix.
+   *
+   * This function is particularly useful when solving for several problems having the same structure.
+   *
+   * \sa factorize()
+   */
+  void analyzePattern(const MatrixType& matrix) {
+    if (m_cholmodFactor) {
+      internal::cm_free_factor<StorageIndex>(m_cholmodFactor, m_cholmod);
+      m_cholmodFactor = 0;
+    }
+    cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
+    m_cholmodFactor = internal::cm_analyze<StorageIndex>(A, m_cholmod);
+
+    this->m_isInitialized = true;
+    this->m_info = Success;
+    m_analysisIsOk = true;
+    m_factorizationIsOk = false;
+  }
+
+  /** Performs a numeric decomposition of \a matrix
+   *
+   * The given matrix must have the same sparsity pattern as the matrix on which the symbolic decomposition has been
+   * performed.
+   *
+   * \sa analyzePattern()
+   */
+  void factorize(const MatrixType& matrix) {
+    eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
+    cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
+    internal::cm_factorize_p<StorageIndex>(&A, m_shiftOffset, 0, 0, m_cholmodFactor, m_cholmod);
+
+    // If the factorization failed, either the input matrix was zero (so m_cholmodFactor == nullptr), or minor is the
+    // column at which it failed. On success minor == n.
+    this->m_info =
+        (m_cholmodFactor != nullptr && m_cholmodFactor->minor == m_cholmodFactor->n ? Success : NumericalIssue);
+    m_factorizationIsOk = true;
+  }
+
+  /** Returns a reference to the Cholmod's configuration structure to get a full control over the performed operations.
+   *  See the Cholmod user guide for details. */
+  cholmod_common& cholmod() { return m_cholmod; }
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+  /** \internal */
+  template <typename Rhs, typename Dest>
+  void _solve_impl(const MatrixBase<Rhs>& b, MatrixBase<Dest>& dest) const {
+    eigen_assert(m_factorizationIsOk &&
+                 "The decomposition is not in a valid state for solving, you must first call either compute() or "
+                 "symbolic()/numeric()");
+    const Index size = m_cholmodFactor->n;
+    EIGEN_UNUSED_VARIABLE(size);
+    eigen_assert(size == b.rows());
+
+    // Cholmod needs column-major storage without inner-stride, which corresponds to the default behavior of Ref.
+    Ref<const Matrix<typename Rhs::Scalar, Dynamic, Dynamic, ColMajor> > b_ref(b.derived());
+
+    cholmod_dense b_cd = viewAsCholmod(b_ref);
+    cholmod_dense* x_cd = internal::cm_solve<StorageIndex>(CHOLMOD_A, *m_cholmodFactor, b_cd, m_cholmod);
+    if (!x_cd) {
+      this->m_info = NumericalIssue;
+      return;
+    }
+    // TODO optimize this copy by swapping when possible (be careful with alignment, etc.)
+    // NOTE Actually, the copy can be avoided by calling cholmod_solve2 instead of cholmod_solve
+    dest = Matrix<Scalar, Dest::RowsAtCompileTime, Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),
+                                                                                 b.rows(), b.cols());
+    internal::cm_free_dense<StorageIndex>(x_cd, m_cholmod);
+  }
+
+  /** \internal */
+  template <typename RhsDerived, typename DestDerived>
+  void _solve_impl(const SparseMatrixBase<RhsDerived>& b, SparseMatrixBase<DestDerived>& dest) const {
+    eigen_assert(m_factorizationIsOk &&
+                 "The decomposition is not in a valid state for solving, you must first call either compute() or "
+                 "symbolic()/numeric()");
+    const Index size = m_cholmodFactor->n;
+    EIGEN_UNUSED_VARIABLE(size);
+    eigen_assert(size == b.rows());
+
+    // note: cs stands for Cholmod Sparse
+    Ref<SparseMatrix<typename RhsDerived::Scalar, ColMajor, typename RhsDerived::StorageIndex> > b_ref(
+        b.const_cast_derived());
+    cholmod_sparse b_cs = viewAsCholmod(b_ref);
+    cholmod_sparse* x_cs = internal::cm_spsolve<StorageIndex>(CHOLMOD_A, *m_cholmodFactor, b_cs, m_cholmod);
+    if (!x_cs) {
+      this->m_info = NumericalIssue;
+      return;
+    }
+    // TODO optimize this copy by swapping when possible (be careful with alignment, etc.)
+    // NOTE cholmod_spsolve in fact just calls the dense solver for blocks of 4 columns at a time (similar to Eigen's
+    // sparse solver)
+    dest.derived() = viewAsEigen<typename DestDerived::Scalar, typename DestDerived::StorageIndex>(*x_cs);
+    internal::cm_free_sparse<StorageIndex>(x_cs, m_cholmod);
+  }
+#endif  // EIGEN_PARSED_BY_DOXYGEN
+
+  /** Sets the shift parameter that will be used to adjust the diagonal coefficients during the numerical factorization.
+   *
+   * During the numerical factorization, an offset term is added to the diagonal coefficients:\n
+   * \c d_ii = \a offset + \c d_ii
+   *
+   * The default is \a offset=0.
+   *
+   * \returns a reference to \c *this.
+   */
+  Derived& setShift(const RealScalar& offset) {
+    m_shiftOffset[0] = double(offset);
+    return derived();
+  }
+
+  /** \returns the determinant of the underlying matrix from the current factorization */
+  Scalar determinant() const {
+    using std::exp;
+    return exp(logDeterminant());
+  }
+
+  /** \returns the log determinant of the underlying matrix from the current factorization */
+  Scalar logDeterminant() const {
+    using numext::real;
+    using std::log;
+    eigen_assert(m_factorizationIsOk &&
+                 "The decomposition is not in a valid state for solving, you must first call either compute() or "
+                 "symbolic()/numeric()");
+
+    RealScalar logDet = 0;
+    Scalar* x = static_cast<Scalar*>(m_cholmodFactor->x);
+    if (m_cholmodFactor->is_super) {
+      // Supernodal factorization stored as a packed list of dense column-major blocks,
+      // as described by the following structure:
+
+      // super[k] == index of the first column of the j-th super node
+      StorageIndex* super = static_cast<StorageIndex*>(m_cholmodFactor->super);
+      // pi[k] == offset to the description of row indices
+      StorageIndex* pi = static_cast<StorageIndex*>(m_cholmodFactor->pi);
+      // px[k] == offset to the respective dense block
+      StorageIndex* px = static_cast<StorageIndex*>(m_cholmodFactor->px);
+
+      Index nb_super_nodes = m_cholmodFactor->nsuper;
+      for (Index k = 0; k < nb_super_nodes; ++k) {
+        StorageIndex ncols = super[k + 1] - super[k];
+        StorageIndex nrows = pi[k + 1] - pi[k];
+
+        Map<const Array<Scalar, 1, Dynamic>, 0, InnerStride<> > sk(x + px[k], ncols, InnerStride<>(nrows + 1));
+        logDet += sk.real().log().sum();
+      }
+    } else {
+      // Simplicial factorization stored as standard CSC matrix.
+      StorageIndex* p = static_cast<StorageIndex*>(m_cholmodFactor->p);
+      Index size = m_cholmodFactor->n;
+      for (Index k = 0; k < size; ++k) logDet += log(real(x[p[k]]));
+    }
+    if (m_cholmodFactor->is_ll) logDet *= 2.0;
+    return logDet;
+  }
+
+  template <typename Stream>
+  void dumpMemory(Stream& /*s*/) {}
+
+ protected:
+  mutable cholmod_common m_cholmod;
+  cholmod_factor* m_cholmodFactor;
+  double m_shiftOffset[2];
+  mutable ComputationInfo m_info;
+  int m_factorizationIsOk;
+  int m_analysisIsOk;
+};
+
+/** \ingroup CholmodSupport_Module
+ * \class CholmodSimplicialLLT
+ * \brief A simplicial direct Cholesky (LLT) factorization and solver based on Cholmod
+ *
+ * This class allows to solve for A.X = B sparse linear problems via a simplicial LL^T Cholesky factorization
+ * using the Cholmod library.
+ * This simplicial variant is equivalent to Eigen's built-in SimplicialLLT class. Therefore, it has little practical
+ * interest. The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices X and B can be
+ * either dense or sparse.
+ *
+ * \tparam MatrixType_ the type of the sparse matrix A, it must be a SparseMatrix<>
+ * \tparam UpLo_ the triangular part that will be used for the computations. It can be Lower
+ *               or Upper. Default is Lower.
+ *
+ * \implsparsesolverconcept
+ *
+ * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non
+ * compressed.
+ *
+ * \warning Only double precision real and complex scalar types are supported by Cholmod.
+ *
+ * \sa \ref TutorialSparseSolverConcept, class CholmodSupernodalLLT, class SimplicialLLT
+ */
+template <typename MatrixType_, int UpLo_ = Lower>
+class CholmodSimplicialLLT : public CholmodBase<MatrixType_, UpLo_, CholmodSimplicialLLT<MatrixType_, UpLo_> > {
+  typedef CholmodBase<MatrixType_, UpLo_, CholmodSimplicialLLT> Base;
+  using Base::m_cholmod;
+
+ public:
+  typedef MatrixType_ MatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  typedef typename MatrixType::StorageIndex StorageIndex;
+  typedef TriangularView<const MatrixType, Eigen::Lower> MatrixL;
+  typedef TriangularView<const typename MatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
+
+  CholmodSimplicialLLT() : Base() { init(); }
+
+  CholmodSimplicialLLT(const MatrixType& matrix) : Base() {
+    init();
+    this->compute(matrix);
+  }
+
+  ~CholmodSimplicialLLT() {}
+
+  /** \returns an expression of the factor L */
+  inline MatrixL matrixL() const { return viewAsEigen<Scalar, StorageIndex>(*Base::m_cholmodFactor); }
+
+  /** \returns an expression of the factor U (= L^*) */
+  inline MatrixU matrixU() const { return matrixL().adjoint(); }
+
+ protected:
+  void init() {
+    m_cholmod.final_asis = 0;
+    m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
+    m_cholmod.final_ll = 1;
+  }
+};
+
+/** \ingroup CholmodSupport_Module
+ * \class CholmodSimplicialLDLT
+ * \brief A simplicial direct Cholesky (LDLT) factorization and solver based on Cholmod
+ *
+ * This class allows to solve for A.X = B sparse linear problems via a simplicial LDL^T Cholesky factorization
+ * using the Cholmod library.
+ * This simplicial variant is equivalent to Eigen's built-in SimplicialLDLT class. Therefore, it has little practical
+ * interest. The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices X and B can be
+ * either dense or sparse.
+ *
+ * \tparam MatrixType_ the type of the sparse matrix A, it must be a SparseMatrix<>
+ * \tparam UpLo_ the triangular part that will be used for the computations. It can be Lower
+ *               or Upper. Default is Lower.
+ *
+ * \implsparsesolverconcept
+ *
+ * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non
+ * compressed.
+ *
+ * \warning Only double precision real and complex scalar types are supported by Cholmod.
+ *
+ * \sa \ref TutorialSparseSolverConcept, class CholmodSupernodalLLT, class SimplicialLDLT
+ */
+template <typename MatrixType_, int UpLo_ = Lower>
+class CholmodSimplicialLDLT : public CholmodBase<MatrixType_, UpLo_, CholmodSimplicialLDLT<MatrixType_, UpLo_> > {
+  typedef CholmodBase<MatrixType_, UpLo_, CholmodSimplicialLDLT> Base;
+  using Base::m_cholmod;
+
+ public:
+  typedef MatrixType_ MatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  typedef typename MatrixType::StorageIndex StorageIndex;
+  typedef Matrix<Scalar, Dynamic, 1> VectorType;
+  typedef TriangularView<const MatrixType, Eigen::UnitLower> MatrixL;
+  typedef TriangularView<const typename MatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
+
+  CholmodSimplicialLDLT() : Base() { init(); }
+
+  CholmodSimplicialLDLT(const MatrixType& matrix) : Base() {
+    init();
+    this->compute(matrix);
+  }
+
+  ~CholmodSimplicialLDLT() {}
+
+  /** \returns a vector expression of the diagonal D */
+  inline VectorType vectorD() const {
+    auto cholmodL = viewAsEigen<Scalar, StorageIndex>(*Base::m_cholmodFactor);
+
+    VectorType D{cholmodL.rows()};
+
+    for (Index k = 0; k < cholmodL.outerSize(); ++k) {
+      typename decltype(cholmodL)::InnerIterator it{cholmodL, k};
+      D(k) = it.value();
+    }
+
+    return D;
+  }
+
+  /** \returns an expression of the factor L */
+  inline MatrixL matrixL() const { return viewAsEigen<Scalar, StorageIndex>(*Base::m_cholmodFactor); }
+
+  /** \returns an expression of the factor U (= L^*) */
+  inline MatrixU matrixU() const { return matrixL().adjoint(); }
+
+ protected:
+  void init() {
+    m_cholmod.final_asis = 1;
+    m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
+  }
+};
+
+/** \ingroup CholmodSupport_Module
+ * \class CholmodSupernodalLLT
+ * \brief A supernodal Cholesky (LLT) factorization and solver based on Cholmod
+ *
+ * This class allows to solve for A.X = B sparse linear problems via a supernodal LL^T Cholesky factorization
+ * using the Cholmod library.
+ * This supernodal variant performs best on dense enough problems, e.g., 3D FEM, or very high order 2D FEM.
+ * The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices
+ * X and B can be either dense or sparse.
+ *
+ * \tparam MatrixType_ the type of the sparse matrix A, it must be a SparseMatrix<>
+ * \tparam UpLo_ the triangular part that will be used for the computations. It can be Lower
+ *               or Upper. Default is Lower.
+ *
+ * \implsparsesolverconcept
+ *
+ * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non
+ * compressed.
+ *
+ * \warning Only double precision real and complex scalar types are supported by Cholmod.
+ *
+ * \sa \ref TutorialSparseSolverConcept
+ */
+template <typename MatrixType_, int UpLo_ = Lower>
+class CholmodSupernodalLLT : public CholmodBase<MatrixType_, UpLo_, CholmodSupernodalLLT<MatrixType_, UpLo_> > {
+  typedef CholmodBase<MatrixType_, UpLo_, CholmodSupernodalLLT> Base;
+  using Base::m_cholmod;
+
+ public:
+  typedef MatrixType_ MatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  typedef typename MatrixType::StorageIndex StorageIndex;
+
+  CholmodSupernodalLLT() : Base() { init(); }
+
+  CholmodSupernodalLLT(const MatrixType& matrix) : Base() {
+    init();
+    this->compute(matrix);
+  }
+
+  ~CholmodSupernodalLLT() {}
+
+  /** \returns an expression of the factor L */
+  inline MatrixType matrixL() const {
+    // Convert Cholmod factor's supernodal storage format to Eigen's CSC storage format
+    cholmod_sparse* cholmodL = internal::cm_factor_to_sparse(*Base::m_cholmodFactor, m_cholmod);
+    MatrixType L = viewAsEigen<Scalar, StorageIndex>(*cholmodL);
+    internal::cm_free_sparse<StorageIndex>(cholmodL, m_cholmod);
+
+    return L;
+  }
+
+  /** \returns an expression of the factor U (= L^*) */
+  inline MatrixType matrixU() const { return matrixL().adjoint(); }
+
+ protected:
+  void init() {
+    m_cholmod.final_asis = 1;
+    m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
+  }
+};
+
+/** \ingroup CholmodSupport_Module
+ * \class CholmodDecomposition
+ * \brief A general Cholesky factorization and solver based on Cholmod
+ *
+ * This class allows to solve for A.X = B sparse linear problems via a LL^T or LDL^T Cholesky factorization
+ * using the Cholmod library. The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices
+ * X and B can be either dense or sparse.
+ *
+ * This variant permits to change the underlying Cholesky method at runtime.
+ * On the other hand, it does not provide access to the result of the factorization.
+ * The default is to let Cholmod automatically choose between a simplicial and supernodal factorization.
+ *
+ * \tparam MatrixType_ the type of the sparse matrix A, it must be a SparseMatrix<>
+ * \tparam UpLo_ the triangular part that will be used for the computations. It can be Lower
+ *               or Upper. Default is Lower.
+ *
+ * \implsparsesolverconcept
+ *
+ * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non
+ * compressed.
+ *
+ * \warning Only double precision real and complex scalar types are supported by Cholmod.
+ *
+ * \sa \ref TutorialSparseSolverConcept
+ */
+template <typename MatrixType_, int UpLo_ = Lower>
+class CholmodDecomposition : public CholmodBase<MatrixType_, UpLo_, CholmodDecomposition<MatrixType_, UpLo_> > {
+  typedef CholmodBase<MatrixType_, UpLo_, CholmodDecomposition> Base;
+  using Base::m_cholmod;
+
+ public:
+  typedef MatrixType_ MatrixType;
+
+  CholmodDecomposition() : Base() { init(); }
+
+  CholmodDecomposition(const MatrixType& matrix) : Base() {
+    init();
+    this->compute(matrix);
+  }
+
+  ~CholmodDecomposition() {}
+
+  void setMode(CholmodMode mode) {
+    switch (mode) {
+      case CholmodAuto:
+        m_cholmod.final_asis = 1;
+        m_cholmod.supernodal = CHOLMOD_AUTO;
+        break;
+      case CholmodSimplicialLLt:
+        m_cholmod.final_asis = 0;
+        m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
+        m_cholmod.final_ll = 1;
+        break;
+      case CholmodSupernodalLLt:
+        m_cholmod.final_asis = 1;
+        m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
+        break;
+      case CholmodLDLt:
+        m_cholmod.final_asis = 1;
+        m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
+        break;
+      default:
+        break;
+    }
+  }
+
+ protected:
+  void init() {
+    m_cholmod.final_asis = 1;
+    m_cholmod.supernodal = CHOLMOD_AUTO;
+  }
+};
+
+}  // end namespace Eigen
+
+#endif  // EIGEN_CHOLMODSUPPORT_H

eigen-master/Eigen/src/CholmodSupport/InternalHeaderCheck.h

+#ifndef EIGEN_CHOLMODSUPPORT_MODULE_H
+#error "Please include Eigen/CholmodSupport instead of including headers inside the src directory directly."
+#endif

eigen-master/Eigen/src/Core/ArithmeticSequence.h

+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2017 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_ARITHMETIC_SEQUENCE_H
+#define EIGEN_ARITHMETIC_SEQUENCE_H
+
+// IWYU pragma: private
+#include "./InternalHeaderCheck.h"
+
+namespace Eigen {
+
+namespace internal {
+
+// Helper to cleanup the type of the increment:
+template <typename T>
+struct cleanup_seq_incr {
+  typedef typename cleanup_index_type<T, DynamicIndex>::type type;
+};
+
+}  // namespace internal
+
+//--------------------------------------------------------------------------------
+// seq(first,last,incr) and seqN(first,size,incr)
+//--------------------------------------------------------------------------------
+
+template <typename FirstType = Index, typename SizeType = Index, typename IncrType = internal::FixedInt<1> >
+class ArithmeticSequence;
+
+template <typename FirstType, typename SizeType, typename IncrType>
+ArithmeticSequence<typename internal::cleanup_index_type<FirstType>::type,
+                   typename internal::cleanup_index_type<SizeType>::type,
+                   typename internal::cleanup_seq_incr<IncrType>::type>
+seqN(FirstType first, SizeType size, IncrType incr);
+
+/** \class ArithmeticSequence
+ * \ingroup Core_Module
+ *
+ * This class represents an arithmetic progression \f$ a_0, a_1, a_2, ..., a_{n-1}\f$ defined by
+ * its \em first value \f$ a_0 \f$, its \em size (aka length) \em n, and the \em increment (aka stride)
+ * that is equal to \f$ a_{i+1}-a_{i}\f$ for any \em i.
+ *
+ * It is internally used as the return type of the Eigen::seq and Eigen::seqN functions, and as the input arguments
+ * of DenseBase::operator()(const RowIndices&, const ColIndices&), and most of the time this is the
+ * only way it is used.
+ *
+ * \tparam FirstType type of the first element, usually an Index,
+ *                   but internally it can be a symbolic expression
+ * \tparam SizeType type representing the size of the sequence, usually an Index
+ *                  or a compile time integral constant. Internally, it can also be a symbolic expression
+ * \tparam IncrType type of the increment, can be a runtime Index, or a compile time integral constant (default is
+ * compile-time 1)
+ *
+ * \sa Eigen::seq, Eigen::seqN, DenseBase::operator()(const RowIndices&, const ColIndices&), class IndexedView
+ */
+template <typename FirstType, typename SizeType, typename IncrType>
+class ArithmeticSequence {
+ public:
+  constexpr ArithmeticSequence() = default;
+  constexpr ArithmeticSequence(FirstType first, SizeType size) : m_first(first), m_size(size) {}
+  constexpr ArithmeticSequence(FirstType first, SizeType size, IncrType incr)
+      : m_first(first), m_size(size), m_incr(incr) {}
+
+  enum {
+    // SizeAtCompileTime = internal::get_fixed_value<SizeType>::value,
+    IncrAtCompileTime = internal::get_fixed_value<IncrType, DynamicIndex>::value
+  };
+
+  /** \returns the size, i.e., number of elements, of the sequence */
+  constexpr Index size() const { return m_size; }
+
+  /** \returns the first element \f$ a_0 \f$ in the sequence */
+  constexpr Index first() const { return m_first; }
+
+  /** \returns the value \f$ a_i \f$ at index \a i in the sequence. */
+  constexpr Index operator[](Index i) const { return m_first + i * m_incr; }
+
+  constexpr const FirstType& firstObject() const { return m_first; }
+  constexpr const SizeType& sizeObject() const { return m_size; }
+  constexpr const IncrType& incrObject() const { return m_incr; }
+
+ protected:
+  FirstType m_first;
+  SizeType m_size;
+  IncrType m_incr;
+
+ public:
+  constexpr auto reverse() const -> decltype(Eigen::seqN(m_first + (m_size + fix<-1>()) * m_incr, m_size, -m_incr)) {
+    return seqN(m_first + (m_size + fix<-1>()) * m_incr, m_size, -m_incr);
+  }
+};
+
+/** \returns an ArithmeticSequence starting at \a first, of length \a size, and increment \a incr
+ *
+ * \sa seqN(FirstType,SizeType), seq(FirstType,LastType,IncrType) */
+template <typename FirstType, typename SizeType, typename IncrType>
+ArithmeticSequence<typename internal::cleanup_index_type<FirstType>::type,
+                   typename internal::cleanup_index_type<SizeType>::type,
+                   typename internal::cleanup_seq_incr<IncrType>::type>
+seqN(FirstType first, SizeType size, IncrType incr) {
+  return ArithmeticSequence<typename internal::cleanup_index_type<FirstType>::type,
+                            typename internal::cleanup_index_type<SizeType>::type,
+                            typename internal::cleanup_seq_incr<IncrType>::type>(first, size, incr);
+}
+
+/** \returns an ArithmeticSequence starting at \a first, of length \a size, and unit increment
+ *
+ * \sa seqN(FirstType,SizeType,IncrType), seq(FirstType,LastType) */
+template <typename FirstType, typename SizeType>
+ArithmeticSequence<typename internal::cleanup_index_type<FirstType>::type,
+                   typename internal::cleanup_index_type<SizeType>::type>
+seqN(FirstType first, SizeType size) {
+  return ArithmeticSequence<typename internal::cleanup_index_type<FirstType>::type,
+                            typename internal::cleanup_index_type<SizeType>::type>(first, size);
+}
+
+#ifdef EIGEN_PARSED_BY_DOXYGEN
+
+/** \returns an ArithmeticSequence starting at \a f, up (or down) to \a l, and with positive (or negative) increment \a
+ * incr
+ *
+ * It is essentially an alias to:
+ * \code
+ * seqN(f, (l-f+incr)/incr, incr);
+ * \endcode
+ *
+ * \sa seqN(FirstType,SizeType,IncrType), seq(FirstType,LastType)
+ */
+template <typename FirstType, typename LastType, typename IncrType>
+auto seq(FirstType f, LastType l, IncrType incr);
+
+/** \returns an ArithmeticSequence starting at \a f, up (or down) to \a l, and unit increment
+ *
+ * It is essentially an alias to:
+ * \code
+ * seqN(f,l-f+1);
+ * \endcode
+ *
+ * \sa seqN(FirstType,SizeType), seq(FirstType,LastType,IncrType)
+ */
+template <typename FirstType, typename LastType>
+auto seq(FirstType f, LastType l);
+
+#else  // EIGEN_PARSED_BY_DOXYGEN
+
+template <typename FirstType, typename LastType>
+auto seq(FirstType f, LastType l)
+    -> decltype(seqN(typename internal::cleanup_index_type<FirstType>::type(f),
+                     (typename internal::cleanup_index_type<LastType>::type(l) -
+                      typename internal::cleanup_index_type<FirstType>::type(f) + fix<1>()))) {
+  return seqN(typename internal::cleanup_index_type<FirstType>::type(f),
+              (typename internal::cleanup_index_type<LastType>::type(l) -
+               typename internal::cleanup_index_type<FirstType>::type(f) + fix<1>()));
+}
+
+template <typename FirstType, typename LastType, typename IncrType>
+auto seq(FirstType f, LastType l, IncrType incr)
+    -> decltype(seqN(typename internal::cleanup_index_type<FirstType>::type(f),
+                     (typename internal::cleanup_index_type<LastType>::type(l) -
+                      typename internal::cleanup_index_type<FirstType>::type(f) +
+                      typename internal::cleanup_seq_incr<IncrType>::type(incr)) /
+                         typename internal::cleanup_seq_incr<IncrType>::type(incr),
+                     typename internal::cleanup_seq_incr<IncrType>::type(incr))) {
+  typedef typename internal::cleanup_seq_incr<IncrType>::type CleanedIncrType;
+  return seqN(typename internal::cleanup_index_type<FirstType>::type(f),
+              (typename internal::cleanup_index_type<LastType>::type(l) -
+               typename internal::cleanup_index_type<FirstType>::type(f) + CleanedIncrType(incr)) /
+                  CleanedIncrType(incr),
+              CleanedIncrType(incr));
+}
+
+#endif  // EIGEN_PARSED_BY_DOXYGEN
+
+namespace placeholders {
+
+/** \cpp11
+ * \returns a symbolic ArithmeticSequence representing the last \a size elements with increment \a incr.
+ *
+ * It is a shortcut for: \code seqN(last-(size-fix<1>)*incr, size, incr) \endcode
+ *
+ * \sa lastN(SizeType), seqN(FirstType,SizeType), seq(FirstType,LastType,IncrType) */
+template <typename SizeType, typename IncrType>
+auto lastN(SizeType size, IncrType incr)
+    -> decltype(seqN(Eigen::placeholders::last - (size - fix<1>()) * incr, size, incr)) {
+  return seqN(Eigen::placeholders::last - (size - fix<1>()) * incr, size, incr);
+}
+
+/** \cpp11
+ * \returns a symbolic ArithmeticSequence representing the last \a size elements with a unit increment.
+ *
+ *  It is a shortcut for: \code seq(last+fix<1>-size, last) \endcode
+ *
+ * \sa lastN(SizeType,IncrType, seqN(FirstType,SizeType), seq(FirstType,LastType) */
+template <typename SizeType>
+auto lastN(SizeType size) -> decltype(seqN(Eigen::placeholders::last + fix<1>() - size, size)) {
+  return seqN(Eigen::placeholders::last + fix<1>() - size, size);
+}
+
+}  // namespace placeholders
+
+/** \namespace Eigen::indexing
+  * \ingroup Core_Module
+  *
+  * The sole purpose of this namespace is to be able to import all functions
+  * and symbols that are expected to be used within operator() for indexing
+  * and slicing. If you already imported the whole Eigen namespace:
+  * \code using namespace Eigen; \endcode
+  * then you are already all set. Otherwise, if you don't want/cannot import
+  * the whole Eigen namespace, the following line:
+  * \code using namespace Eigen::indexing; \endcode
+  * is equivalent to:
+  * \code
+  using Eigen::fix;
+  using Eigen::seq;
+  using Eigen::seqN;
+  using Eigen::placeholders::all;
+  using Eigen::placeholders::last;
+  using Eigen::placeholders::lastN;  // c++11 only
+  using Eigen::placeholders::lastp1;
+  \endcode
+  */
+namespace indexing {
+using Eigen::fix;
+using Eigen::seq;
+using Eigen::seqN;
+using Eigen::placeholders::all;
+using Eigen::placeholders::last;
+using Eigen::placeholders::lastN;
+using Eigen::placeholders::lastp1;
+}  // namespace indexing
+
+}  // end namespace Eigen
+
+#endif  // EIGEN_ARITHMETIC_SEQUENCE_H

eigen-master/Eigen/src/Core/Array.h

+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_ARRAY_H
+#define EIGEN_ARRAY_H
+
+// IWYU pragma: private
+#include "./InternalHeaderCheck.h"
+
+namespace Eigen {
+
+namespace internal {
+template <typename Scalar_, int Rows_, int Cols_, int Options_, int MaxRows_, int MaxCols_>
+struct traits<Array<Scalar_, Rows_, Cols_, Options_, MaxRows_, MaxCols_>>
+    : traits<Matrix<Scalar_, Rows_, Cols_, Options_, MaxRows_, MaxCols_>> {
+  typedef ArrayXpr XprKind;
+  typedef ArrayBase<Array<Scalar_, Rows_, Cols_, Options_, MaxRows_, MaxCols_>> XprBase;
+};
+}  // namespace internal
+
+/** \class Array
+ * \ingroup Core_Module
+ *
+ * \brief General-purpose arrays with easy API for coefficient-wise operations
+ *
+ * The %Array class is very similar to the Matrix class. It provides
+ * general-purpose one- and two-dimensional arrays. The difference between the
+ * %Array and the %Matrix class is primarily in the API: the API for the
+ * %Array class provides easy access to coefficient-wise operations, while the
+ * API for the %Matrix class provides easy access to linear-algebra
+ * operations.
+ *
+ * See documentation of class Matrix for detailed information on the template parameters
+ * storage layout.
+ *
+ * This class can be extended with the help of the plugin mechanism described on the page
+ * \ref TopicCustomizing_Plugins by defining the preprocessor symbol \c EIGEN_ARRAY_PLUGIN.
+ *
+ * \sa \blank \ref TutorialArrayClass, \ref TopicClassHierarchy
+ */
+template <typename Scalar_, int Rows_, int Cols_, int Options_, int MaxRows_, int MaxCols_>
+class Array : public PlainObjectBase<Array<Scalar_, Rows_, Cols_, Options_, MaxRows_, MaxCols_>> {
+ public:
+  typedef PlainObjectBase<Array> Base;
+  EIGEN_DENSE_PUBLIC_INTERFACE(Array)
+
+  enum { Options = Options_ };
+  typedef typename Base::PlainObject PlainObject;
+
+ protected:
+  template <typename Derived, typename OtherDerived, bool IsVector>
+  friend struct internal::conservative_resize_like_impl;
+
+  using Base::m_storage;
+
+ public:
+  using Base::base;
+  using Base::coeff;
+  using Base::coeffRef;
+
+  /**
+   * The usage of
+   *   using Base::operator=;
+   * fails on MSVC. Since the code below is working with GCC and MSVC, we skipped
+   * the usage of 'using'. This should be done only for operator=.
+   */
+  template <typename OtherDerived>
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Array& operator=(const EigenBase<OtherDerived>& other) {
+    return Base::operator=(other);
+  }
+
+  /** Set all the entries to \a value.
+   * \sa DenseBase::setConstant(), DenseBase::fill()
+   */
+  /* This overload is needed because the usage of
+   *   using Base::operator=;
+   * fails on MSVC. Since the code below is working with GCC and MSVC, we skipped
+   * the usage of 'using'. This should be done only for operator=.
+   */
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Array& operator=(const Scalar& value) {
+    Base::setConstant(value);
+    return *this;
+  }
+
+  /** Copies the value of the expression \a other into \c *this with automatic resizing.
+   *
+   * *this might be resized to match the dimensions of \a other. If *this was a null matrix (not already initialized),
+   * it will be initialized.
+   *
+   * Note that copying a row-vector into a vector (and conversely) is allowed.
+   * The resizing, if any, is then done in the appropriate way so that row-vectors
+   * remain row-vectors and vectors remain vectors.
+   */
+  template <typename OtherDerived>
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Array& operator=(const DenseBase<OtherDerived>& other) {
+    return Base::_set(other);
+  }
+
+  /**
+   * \brief Assigns arrays to each other.
+   *
+   * \note This is a special case of the templated operator=. Its purpose is
+   * to prevent a default operator= from hiding the templated operator=.
+   *
+   * \callgraph
+   */
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Array& operator=(const Array& other) { return Base::_set(other); }
+
+  /** Default constructor.
+   *
+   * For fixed-size matrices, does nothing.
+   *
+   * For dynamic-size matrices, creates an empty matrix of size 0. Does not allocate any array. Such a matrix
+   * is called a null matrix. This constructor is the unique way to create null matrices: resizing
+   * a matrix to 0 is not supported.
+   *
+   * \sa resize(Index,Index)
+   */
+#ifdef EIGEN_INITIALIZE_COEFFS
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Array() : Base() { EIGEN_INITIALIZE_COEFFS_IF_THAT_OPTION_IS_ENABLED }
+#else
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Array() = default;
+#endif
+  /** \brief Move constructor */
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Array(Array&&) = default;
+  EIGEN_DEVICE_FUNC Array& operator=(Array&& other) noexcept(std::is_nothrow_move_assignable<Scalar>::value) {
+    Base::operator=(std::move(other));
+    return *this;
+  }
+
+  /** \brief Construct a row of column vector with fixed size from an arbitrary number of coefficients.
+   *
+   * \only_for_vectors
+   *
+   * This constructor is for 1D array or vectors with more than 4 coefficients.
+   *
+   * \warning To construct a column (resp. row) vector of fixed length, the number of values passed to this
+   * constructor must match the the fixed number of rows (resp. columns) of \c *this.
+   *
+   *
+   * Example: \include Array_variadic_ctor_cxx11.cpp
+   * Output: \verbinclude Array_variadic_ctor_cxx11.out
+   *
+   * \sa Array(const std::initializer_list<std::initializer_list<Scalar>>&)
+   * \sa Array(const Scalar&), Array(const Scalar&,const Scalar&)
+   */
+  template <typename... ArgTypes>
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Array(const Scalar& a0, const Scalar& a1, const Scalar& a2, const Scalar& a3,
+                                              const ArgTypes&... args)
+      : Base(a0, a1, a2, a3, args...) {}
+
+  /** \brief Constructs an array and initializes it from the coefficients given as initializer-lists grouped by row.
+   * \cpp11
+   *
+   * In the general case, the constructor takes a list of rows, each row being represented as a list of coefficients:
+   *
+   * Example: \include Array_initializer_list_23_cxx11.cpp
+   * Output: \verbinclude Array_initializer_list_23_cxx11.out
+   *
+   * Each of the inner initializer lists must contain the exact same number of elements, otherwise an assertion is
+   * triggered.
+   *
+   * In the case of a compile-time column 1D array, implicit transposition from a single row is allowed.
+   * Therefore <code> Array<int,Dynamic,1>{{1,2,3,4,5}}</code> is legal and the more verbose syntax
+   * <code>Array<int,Dynamic,1>{{1},{2},{3},{4},{5}}</code> can be avoided:
+   *
+   * Example: \include Array_initializer_list_vector_cxx11.cpp
+   * Output: \verbinclude Array_initializer_list_vector_cxx11.out
+   *
+   * In the case of fixed-sized arrays, the initializer list sizes must exactly match the array sizes,
+   * and implicit transposition is allowed for compile-time 1D arrays only.
+   *
+   * \sa  Array(const Scalar& a0, const Scalar& a1, const Scalar& a2, const Scalar& a3, const ArgTypes&... args)
+   */
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Array(
+      const std::initializer_list<std::initializer_list<Scalar>>& list)
+      : Base(list) {}
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+  template <typename T>
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE explicit Array(const T& x) {
+    Base::template _init1<T>(x);
+  }
+
+  template <typename T0, typename T1>
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Array(const T0& val0, const T1& val1) {
+    this->template _init2<T0, T1>(val0, val1);
+  }
+
+#else
+  /** \brief Constructs a fixed-sized array initialized with coefficients starting at \a data */
+  EIGEN_DEVICE_FUNC explicit Array(const Scalar* data);
+  /** Constructs a vector or row-vector with given dimension. \only_for_vectors
+   *
+   * Note that this is only useful for dynamic-size vectors. For fixed-size vectors,
+   * it is redundant to pass the dimension here, so it makes more sense to use the default
+   * constructor Array() instead.
+   */
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE explicit Array(Index dim);
+  /** constructs an initialized 1x1 Array with the given coefficient
+   * \sa const Scalar& a0, const Scalar& a1, const Scalar& a2, const Scalar& a3, const ArgTypes&... args */
+  Array(const Scalar& value);
+  /** constructs an uninitialized array with \a rows rows and \a cols columns.
+   *
+   * This is useful for dynamic-size arrays. For fixed-size arrays,
+   * it is redundant to pass these parameters, so one should use the default constructor
+   * Array() instead. */
+  Array(Index rows, Index cols);
+  /** constructs an initialized 2D vector with given coefficients
+   * \sa Array(const Scalar& a0, const Scalar& a1, const Scalar& a2, const Scalar& a3, const ArgTypes&... args) */
+  Array(const Scalar& val0, const Scalar& val1);
+#endif  // end EIGEN_PARSED_BY_DOXYGEN
+
+  /** constructs an initialized 3D vector with given coefficients
+   * \sa Array(const Scalar& a0, const Scalar& a1, const Scalar& a2, const Scalar& a3, const ArgTypes&... args)
+   */
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Array(const Scalar& val0, const Scalar& val1, const Scalar& val2) {
+    EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Array, 3)
+    m_storage.data()[0] = val0;
+    m_storage.data()[1] = val1;
+    m_storage.data()[2] = val2;
+  }
+  /** constructs an initialized 4D vector with given coefficients
+   * \sa Array(const Scalar& a0, const Scalar& a1, const Scalar& a2, const Scalar& a3, const ArgTypes&... args)
+   */
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Array(const Scalar& val0, const Scalar& val1, const Scalar& val2,
+                                              const Scalar& val3) {
+    EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Array, 4)
+    m_storage.data()[0] = val0;
+    m_storage.data()[1] = val1;
+    m_storage.data()[2] = val2;
+    m_storage.data()[3] = val3;
+  }
+
+  /** Copy constructor */
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE constexpr Array(const Array&) = default;
+
+ private:
+  struct PrivateType {};
+
+ public:
+  /** \sa MatrixBase::operator=(const EigenBase<OtherDerived>&) */
+  template <typename OtherDerived>
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Array(
+      const EigenBase<OtherDerived>& other,
+      std::enable_if_t<internal::is_convertible<typename OtherDerived::Scalar, Scalar>::value, PrivateType> =
+          PrivateType())
+      : Base(other.derived()) {}
+
+  EIGEN_DEVICE_FUNC constexpr Index innerStride() const noexcept { return 1; }
+  EIGEN_DEVICE_FUNC constexpr Index outerStride() const noexcept { return this->innerSize(); }
+
+#ifdef EIGEN_ARRAY_PLUGIN
+#include EIGEN_ARRAY_PLUGIN
+#endif
+
+ private:
+  template <typename MatrixType, typename OtherDerived, bool SwapPointers>
+  friend struct internal::matrix_swap_impl;
+};
+
+/** \defgroup arraytypedefs Global array typedefs
+ * \ingroup Core_Module
+ *
+ * %Eigen defines several typedef shortcuts for most common 1D and 2D array types.
+ *
+ * The general patterns are the following:
+ *
+ * \c ArrayRowsColsType where \c Rows and \c Cols can be \c 2,\c 3,\c 4 for fixed size square matrices or \c X for
+ * dynamic size, and where \c Type can be \c i for integer, \c f for float, \c d for double, \c cf for complex float, \c
+ * cd for complex double.
+ *
+ * For example, \c Array33d is a fixed-size 3x3 array type of doubles, and \c ArrayXXf is a dynamic-size matrix of
+ * floats.
+ *
+ * There are also \c ArraySizeType which are self-explanatory. For example, \c Array4cf is
+ * a fixed-size 1D array of 4 complex floats.
+ *
+ * With \cpp11, template alias are also defined for common sizes.
+ * They follow the same pattern as above except that the scalar type suffix is replaced by a
+ * template parameter, i.e.:
+ *   - `ArrayRowsCols<Type>` where `Rows` and `Cols` can be \c 2,\c 3,\c 4, or \c X for fixed or dynamic size.
+ *   - `ArraySize<Type>` where `Size` can be \c 2,\c 3,\c 4 or \c X for fixed or dynamic size 1D arrays.
+ *
+ * \sa class Array
+ */
+
+#define EIGEN_MAKE_ARRAY_TYPEDEFS(Type, TypeSuffix, Size, SizeSuffix)        \
+  /** \ingroup arraytypedefs */                                              \
+  typedef Array<Type, Size, Size> Array##SizeSuffix##SizeSuffix##TypeSuffix; \
+  /** \ingroup arraytypedefs */                                              \
+  typedef Array<Type, Size, 1> Array##SizeSuffix##TypeSuffix;
+
+#define EIGEN_MAKE_ARRAY_FIXED_TYPEDEFS(Type, TypeSuffix, Size)  \
+  /** \ingroup arraytypedefs */                                  \
+  typedef Array<Type, Size, Dynamic> Array##Size##X##TypeSuffix; \
+  /** \ingroup arraytypedefs */                                  \
+  typedef Array<Type, Dynamic, Size> Array##X##Size##TypeSuffix;
+
+#define EIGEN_MAKE_ARRAY_TYPEDEFS_ALL_SIZES(Type, TypeSuffix) \
+  EIGEN_MAKE_ARRAY_TYPEDEFS(Type, TypeSuffix, 2, 2)           \
+  EIGEN_MAKE_ARRAY_TYPEDEFS(Type, TypeSuffix, 3, 3)           \
+  EIGEN_MAKE_ARRAY_TYPEDEFS(Type, TypeSuffix, 4, 4)           \
+  EIGEN_MAKE_ARRAY_TYPEDEFS(Type, TypeSuffix, Dynamic, X)     \
+  EIGEN_MAKE_ARRAY_FIXED_TYPEDEFS(Type, TypeSuffix, 2)        \
+  EIGEN_MAKE_ARRAY_FIXED_TYPEDEFS(Type, TypeSuffix, 3)        \
+  EIGEN_MAKE_ARRAY_FIXED_TYPEDEFS(Type, TypeSuffix, 4)
+
+EIGEN_MAKE_ARRAY_TYPEDEFS_ALL_SIZES(int, i)
+EIGEN_MAKE_ARRAY_TYPEDEFS_ALL_SIZES(float, f)
+EIGEN_MAKE_ARRAY_TYPEDEFS_ALL_SIZES(double, d)
+EIGEN_MAKE_ARRAY_TYPEDEFS_ALL_SIZES(std::complex<float>, cf)
+EIGEN_MAKE_ARRAY_TYPEDEFS_ALL_SIZES(std::complex<double>, cd)
+
+#undef EIGEN_MAKE_ARRAY_TYPEDEFS_ALL_SIZES
+#undef EIGEN_MAKE_ARRAY_TYPEDEFS
+#undef EIGEN_MAKE_ARRAY_FIXED_TYPEDEFS
+
+#define EIGEN_MAKE_ARRAY_TYPEDEFS(Size, SizeSuffix)              \
+  /** \ingroup arraytypedefs */                                  \
+  /** \brief \cpp11 */                                           \
+  template <typename Type>                                       \
+  using Array##SizeSuffix##SizeSuffix = Array<Type, Size, Size>; \
+  /** \ingroup arraytypedefs */                                  \
+  /** \brief \cpp11 */                                           \
+  template <typename Type>                                       \
+  using Array##SizeSuffix = Array<Type, Size, 1>;
+
+#define EIGEN_MAKE_ARRAY_FIXED_TYPEDEFS(Size)        \
+  /** \ingroup arraytypedefs */                      \
+  /** \brief \cpp11 */                               \
+  template <typename Type>                           \
+  using Array##Size##X = Array<Type, Size, Dynamic>; \
+  /** \ingroup arraytypedefs */                      \
+  /** \brief \cpp11 */                               \
+  template <typename Type>                           \
+  using Array##X##Size = Array<Type, Dynamic, Size>;
+
+EIGEN_MAKE_ARRAY_TYPEDEFS(2, 2)
+EIGEN_MAKE_ARRAY_TYPEDEFS(3, 3)
+EIGEN_MAKE_ARRAY_TYPEDEFS(4, 4)
+EIGEN_MAKE_ARRAY_TYPEDEFS(Dynamic, X)
+EIGEN_MAKE_ARRAY_FIXED_TYPEDEFS(2)
+EIGEN_MAKE_ARRAY_FIXED_TYPEDEFS(3)
+EIGEN_MAKE_ARRAY_FIXED_TYPEDEFS(4)
+
+#undef EIGEN_MAKE_ARRAY_TYPEDEFS
+#undef EIGEN_MAKE_ARRAY_FIXED_TYPEDEFS
+
+#define EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE_AND_SIZE(TypeSuffix, SizeSuffix) \
+  using Eigen::Matrix##SizeSuffix##TypeSuffix;                               \
+  using Eigen::Vector##SizeSuffix##TypeSuffix;                               \
+  using Eigen::RowVector##SizeSuffix##TypeSuffix;
+
+#define EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE(TypeSuffix)       \
+  EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE_AND_SIZE(TypeSuffix, 2) \
+  EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE_AND_SIZE(TypeSuffix, 3) \
+  EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE_AND_SIZE(TypeSuffix, 4) \
+  EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE_AND_SIZE(TypeSuffix, X)
+
+#define EIGEN_USING_ARRAY_TYPEDEFS        \
+  EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE(i)  \
+  EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE(f)  \
+  EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE(d)  \
+  EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE(cf) \
+  EIGEN_USING_ARRAY_TYPEDEFS_FOR_TYPE(cd)
+
+}  // end namespace Eigen
+
+#endif  // EIGEN_ARRAY_H

eigen-master/Eigen/src/Core/ArrayBase.h

+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_ARRAYBASE_H
+#define EIGEN_ARRAYBASE_H
+
+// IWYU pragma: private
+#include "./InternalHeaderCheck.h"
+
+namespace Eigen {
+
+template <typename ExpressionType>
+class MatrixWrapper;
+
+/** \class ArrayBase
+ * \ingroup Core_Module
+ *
+ * \brief Base class for all 1D and 2D array, and related expressions
+ *
+ * An array is similar to a dense vector or matrix. While matrices are mathematical
+ * objects with well defined linear algebra operators, an array is just a collection
+ * of scalar values arranged in a one or two dimensional fashion. As the main consequence,
+ * all operations applied to an array are performed coefficient wise. Furthermore,
+ * arrays support scalar math functions of the c++ standard library (e.g., std::sin(x)), and convenient
+ * constructors allowing to easily write generic code working for both scalar values
+ * and arrays.
+ *
+ * This class is the base that is inherited by all array expression types.
+ *
+ * \tparam Derived is the derived type, e.g., an array or an expression type.
+ *
+ * This class can be extended with the help of the plugin mechanism described on the page
+ * \ref TopicCustomizing_Plugins by defining the preprocessor symbol \c EIGEN_ARRAYBASE_PLUGIN.
+ *
+ * \sa class MatrixBase, \ref TopicClassHierarchy
+ */
+template <typename Derived>
+class ArrayBase : public DenseBase<Derived> {
+ public:
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+  /** The base class for a given storage type. */
+  typedef ArrayBase StorageBaseType;
+
+  typedef ArrayBase Eigen_BaseClassForSpecializationOfGlobalMathFuncImpl;
+
+  typedef typename internal::traits<Derived>::StorageKind StorageKind;
+  typedef typename internal::traits<Derived>::Scalar Scalar;
+  typedef typename internal::packet_traits<Scalar>::type PacketScalar;
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+
+  typedef DenseBase<Derived> Base;
+  using Base::ColsAtCompileTime;
+  using Base::Flags;
+  using Base::IsVectorAtCompileTime;
+  using Base::MaxColsAtCompileTime;
+  using Base::MaxRowsAtCompileTime;
+  using Base::MaxSizeAtCompileTime;
+  using Base::RowsAtCompileTime;
+  using Base::SizeAtCompileTime;
+
+  using Base::coeff;
+  using Base::coeffRef;
+  using Base::cols;
+  using Base::const_cast_derived;
+  using Base::derived;
+  using Base::lazyAssign;
+  using Base::rows;
+  using Base::size;
+  using Base::operator-;
+  using Base::operator=;
+  using Base::operator+=;
+  using Base::operator-=;
+  using Base::operator*=;
+  using Base::operator/=;
+
+  typedef typename Base::CoeffReturnType CoeffReturnType;
+
+  typedef typename Base::PlainObject PlainObject;
+
+  /** \internal Represents a matrix with all coefficients equal to one another*/
+  typedef CwiseNullaryOp<internal::scalar_constant_op<Scalar>, PlainObject> ConstantReturnType;
+#endif  // not EIGEN_PARSED_BY_DOXYGEN
+
+#define EIGEN_CURRENT_STORAGE_BASE_CLASS Eigen::ArrayBase
+#define EIGEN_DOC_UNARY_ADDONS(X, Y)
+#include "../plugins/MatrixCwiseUnaryOps.inc"
+#include "../plugins/ArrayCwiseUnaryOps.inc"
+#include "../plugins/CommonCwiseBinaryOps.inc"
+#include "../plugins/MatrixCwiseBinaryOps.inc"
+#include "../plugins/ArrayCwiseBinaryOps.inc"
+#ifdef EIGEN_ARRAYBASE_PLUGIN
+#include EIGEN_ARRAYBASE_PLUGIN
+#endif
+#undef EIGEN_CURRENT_STORAGE_BASE_CLASS
+#undef EIGEN_DOC_UNARY_ADDONS
+
+  /** Special case of the template operator=, in order to prevent the compiler
+   * from generating a default operator= (issue hit with g++ 4.1)
+   */
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Derived& operator=(const ArrayBase& other) {
+    internal::call_assignment(derived(), other.derived());
+    return derived();
+  }
+
+  /** Set all the entries to \a value.
+   * \sa DenseBase::setConstant(), DenseBase::fill() */
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Derived& operator=(const Scalar& value) {
+    Base::setConstant(value);
+    return derived();
+  }
+
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Derived& operator+=(const Scalar& other) {
+    internal::call_assignment(this->derived(), PlainObject::Constant(rows(), cols(), other),
+                              internal::add_assign_op<Scalar, Scalar>());
+    return derived();
+  }
+
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Derived& operator-=(const Scalar& other) {
+    internal::call_assignment(this->derived(), PlainObject::Constant(rows(), cols(), other),
+                              internal::sub_assign_op<Scalar, Scalar>());
+    return derived();
+  }
+
+  /** replaces \c *this by \c *this + \a other.
+   *
+   * \returns a reference to \c *this
+   */
+  template <typename OtherDerived>
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Derived& operator+=(const ArrayBase<OtherDerived>& other) {
+    call_assignment(derived(), other.derived(), internal::add_assign_op<Scalar, typename OtherDerived::Scalar>());
+    return derived();
+  }
+
+  /** replaces \c *this by \c *this - \a other.
+   *
+   * \returns a reference to \c *this
+   */
+  template <typename OtherDerived>
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Derived& operator-=(const ArrayBase<OtherDerived>& other) {
+    call_assignment(derived(), other.derived(), internal::sub_assign_op<Scalar, typename OtherDerived::Scalar>());
+    return derived();
+  }
+
+  /** replaces \c *this by \c *this * \a other coefficient wise.
+   *
+   * \returns a reference to \c *this
+   */
+  template <typename OtherDerived>
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Derived& operator*=(const ArrayBase<OtherDerived>& other) {
+    call_assignment(derived(), other.derived(), internal::mul_assign_op<Scalar, typename OtherDerived::Scalar>());
+    return derived();
+  }
+
+  /** replaces \c *this by \c *this / \a other coefficient wise.
+   *
+   * \returns a reference to \c *this
+   */
+  template <typename OtherDerived>
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Derived& operator/=(const ArrayBase<OtherDerived>& other) {
+    call_assignment(derived(), other.derived(), internal::div_assign_op<Scalar, typename OtherDerived::Scalar>());
+    return derived();
+  }
+
+ public:
+  EIGEN_DEVICE_FUNC ArrayBase<Derived>& array() { return *this; }
+  EIGEN_DEVICE_FUNC const ArrayBase<Derived>& array() const { return *this; }
+
+  /** \returns an \link Eigen::MatrixBase Matrix \endlink expression of this array
+   * \sa MatrixBase::array() */
+  EIGEN_DEVICE_FUNC MatrixWrapper<Derived> matrix() { return MatrixWrapper<Derived>(derived()); }
+  EIGEN_DEVICE_FUNC const MatrixWrapper<const Derived> matrix() const {
+    return MatrixWrapper<const Derived>(derived());
+  }
+
+  //     template<typename Dest>
+  //     inline void evalTo(Dest& dst) const { dst = matrix(); }
+
+ protected:
+  EIGEN_DEFAULT_COPY_CONSTRUCTOR(ArrayBase)
+  EIGEN_DEFAULT_EMPTY_CONSTRUCTOR_AND_DESTRUCTOR(ArrayBase)
+
+ private:
+  explicit ArrayBase(Index);
+  ArrayBase(Index, Index);
+  template <typename OtherDerived>
+  explicit ArrayBase(const ArrayBase<OtherDerived>&);
+
+ protected:
+  // mixing arrays and matrices is not legal
+  template <typename OtherDerived>
+  Derived& operator+=(const MatrixBase<OtherDerived>&) {
+    EIGEN_STATIC_ASSERT(std::ptrdiff_t(sizeof(typename OtherDerived::Scalar)) == -1,
+                        YOU_CANNOT_MIX_ARRAYS_AND_MATRICES);
+    return *this;
+  }
+  // mixing arrays and matrices is not legal
+  template <typename OtherDerived>
+  Derived& operator-=(const MatrixBase<OtherDerived>&) {
+    EIGEN_STATIC_ASSERT(std::ptrdiff_t(sizeof(typename OtherDerived::Scalar)) == -1,
+                        YOU_CANNOT_MIX_ARRAYS_AND_MATRICES);
+    return *this;
+  }
+};
+
+}  // end namespace Eigen
+
+#endif  // EIGEN_ARRAYBASE_H