diff options
| -rw-r--r-- | Eigen/OrderingMethods | 3 | ||||
| -rw-r--r-- | Eigen/SparseCholesky | 8 | ||||
| -rw-r--r-- | Eigen/src/OrderingMethods/Amd.h | 4 | ||||
| -rw-r--r-- | Eigen/src/OrderingMethods/Ordering.h | 26 | ||||
| -rw-r--r-- | Eigen/src/SparseCholesky/SimplicialCholesky.h | 189 | ||||
| -rw-r--r-- | Eigen/src/SparseCholesky/SimplicialCholesky_impl.h | 199 | ||||
| -rw-r--r-- | unsupported/Eigen/IterativeSolvers | 3 | ||||
| -rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h | 4 | ||||
| -rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/IterationController.h | 4 |
9 files changed, 221 insertions, 219 deletions
diff --git a/Eigen/OrderingMethods b/Eigen/OrderingMethods index 423cfb0cd..7c0f1ffff 100644 --- a/Eigen/OrderingMethods +++ b/Eigen/OrderingMethods @@ -56,7 +56,10 @@ * \endcode */ +#ifndef EIGEN_MPL2_ONLY #include "src/OrderingMethods/Amd.h" +#endif + #include "src/OrderingMethods/Ordering.h" #include "src/Core/util/ReenableStupidWarnings.h" diff --git a/Eigen/SparseCholesky b/Eigen/SparseCholesky index f94e2e765..800f17bc4 100644 --- a/Eigen/SparseCholesky +++ b/Eigen/SparseCholesky @@ -20,11 +20,19 @@ * \endcode */ +#ifdef EIGEN_MPL2_ONLY +#error The SparseCholesky module has nothing to offer in MPL2 only mode +#endif + #include "src/misc/Solve.h" #include "src/misc/SparseSolve.h" #include "src/SparseCholesky/SimplicialCholesky.h" +#ifndef EIGEN_MPL2_ONLY +#include "src/SparseCholesky/SimplicialCholesky_impl.h" +#endif + #include "src/Core/util/ReenableStupidWarnings.h" #endif // EIGEN_SPARSECHOLESKY_MODULE_H diff --git a/Eigen/src/OrderingMethods/Amd.h b/Eigen/src/OrderingMethods/Amd.h index 8878ef863..2ef6aa64c 100644 --- a/Eigen/src/OrderingMethods/Amd.h +++ b/Eigen/src/OrderingMethods/Amd.h @@ -2,10 +2,6 @@ // for linear algebra. // // Copyright (C) 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/. /* diff --git a/Eigen/src/OrderingMethods/Ordering.h b/Eigen/src/OrderingMethods/Ordering.h index 2471316b9..fb1da1fe2 100644 --- a/Eigen/src/OrderingMethods/Ordering.h +++ b/Eigen/src/OrderingMethods/Ordering.h @@ -4,29 +4,13 @@ // // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> // -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see <http://www.gnu.org/licenses/>. +// 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_ORDERING_H #define EIGEN_ORDERING_H -#include "Amd.h" namespace Eigen { #include "Eigen_Colamd.h" @@ -53,6 +37,8 @@ void ordering_helper_at_plus_a(const MatrixType& mat, MatrixType& symmat) } +#ifndef EIGEN_MPL2_ONLY + /** \ingroup OrderingMethods_Module * \class AMDOrdering * @@ -94,6 +80,8 @@ class AMDOrdering } }; +#endif // EIGEN_MPL2_ONLY + /** \ingroup OrderingMethods_Module * \class NaturalOrdering * diff --git a/Eigen/src/SparseCholesky/SimplicialCholesky.h b/Eigen/src/SparseCholesky/SimplicialCholesky.h index 59bddb1e4..62747279d 100644 --- a/Eigen/src/SparseCholesky/SimplicialCholesky.h +++ b/Eigen/src/SparseCholesky/SimplicialCholesky.h @@ -1,52 +1,12 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2008-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/. -/* - -NOTE: the _symbolic, and _numeric functions has been adapted from - the LDL library: - -LDL Copyright (c) 2005 by Timothy A. Davis. All Rights Reserved. - -LDL License: - - Your use or distribution of LDL or any modified version of - LDL implies that you agree to this License. - - This library is free software; you can redistribute it and/or - modify it under the terms of the GNU Lesser General Public - License as published by the Free Software Foundation; either - version 2.1 of the License, or (at your option) any later version. - - This library is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public - License along with this library; if not, write to the Free Software - Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 - USA - - Permission is hereby granted to use or copy this program under the - terms of the GNU LGPL, provided that the Copyright, this License, - and the Availability of the original version is retained on all copies. - User documentation of any code that uses this code or any modified - version of this code must cite the Copyright, this License, the - Availability note, and "Used by permission." Permission to modify - the code and to distribute modified code is granted, provided the - Copyright, this License, and the Availability note are retained, - and a notice that the code was modified is included. - */ - -#include "../Core/util/NonMPL2.h" - #ifndef EIGEN_SIMPLICIAL_CHOLESKY_H #define EIGEN_SIMPLICIAL_CHOLESKY_H @@ -672,153 +632,6 @@ void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, CholMatrixTy ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); } -template<typename Derived> -void SimplicialCholeskyBase<Derived>::analyzePattern_preordered(const CholMatrixType& ap, bool doLDLT) -{ - const Index size = ap.rows(); - m_matrix.resize(size, size); - m_parent.resize(size); - m_nonZerosPerCol.resize(size); - - ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0); - - for(Index k = 0; k < size; ++k) - { - /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */ - m_parent[k] = -1; /* parent of k is not yet known */ - tags[k] = k; /* mark node k as visited */ - m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */ - for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it) - { - Index i = it.index(); - if(i < k) - { - /* follow path from i to root of etree, stop at flagged node */ - for(; tags[i] != k; i = m_parent[i]) - { - /* find parent of i if not yet determined */ - if (m_parent[i] == -1) - m_parent[i] = k; - m_nonZerosPerCol[i]++; /* L (k,i) is nonzero */ - tags[i] = k; /* mark i as visited */ - } - } - } - } - - /* construct Lp index array from m_nonZerosPerCol column counts */ - Index* Lp = m_matrix.outerIndexPtr(); - Lp[0] = 0; - for(Index k = 0; k < size; ++k) - Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1); - - m_matrix.resizeNonZeros(Lp[size]); - - m_isInitialized = true; - m_info = Success; - m_analysisIsOk = true; - m_factorizationIsOk = false; -} - - -template<typename Derived> -template<bool DoLDLT> -void SimplicialCholeskyBase<Derived>::factorize_preordered(const CholMatrixType& ap) -{ - using std::sqrt; - - eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); - eigen_assert(ap.rows()==ap.cols()); - const Index size = ap.rows(); - eigen_assert(m_parent.size()==size); - eigen_assert(m_nonZerosPerCol.size()==size); - - const Index* Lp = m_matrix.outerIndexPtr(); - Index* Li = m_matrix.innerIndexPtr(); - Scalar* Lx = m_matrix.valuePtr(); - - ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0); - ei_declare_aligned_stack_constructed_variable(Index, pattern, size, 0); - ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0); - - bool ok = true; - m_diag.resize(DoLDLT ? size : 0); - - for(Index k = 0; k < size; ++k) - { - // compute nonzero pattern of kth row of L, in topological order - y[k] = 0.0; // Y(0:k) is now all zero - Index top = size; // stack for pattern is empty - tags[k] = k; // mark node k as visited - m_nonZerosPerCol[k] = 0; // count of nonzeros in column k of L - for(typename MatrixType::InnerIterator it(ap,k); it; ++it) - { - Index i = it.index(); - if(i <= k) - { - y[i] += internal::conj(it.value()); /* scatter A(i,k) into Y (sum duplicates) */ - Index len; - for(len = 0; tags[i] != k; i = m_parent[i]) - { - pattern[len++] = i; /* L(k,i) is nonzero */ - tags[i] = k; /* mark i as visited */ - } - while(len > 0) - pattern[--top] = pattern[--len]; - } - } - - /* compute numerical values kth row of L (a sparse triangular solve) */ - - RealScalar d = internal::real(y[k]) * m_shiftScale + m_shiftOffset; // get D(k,k), apply the shift function, and clear Y(k) - y[k] = 0.0; - for(; top < size; ++top) - { - Index i = pattern[top]; /* pattern[top:n-1] is pattern of L(:,k) */ - Scalar yi = y[i]; /* get and clear Y(i) */ - y[i] = 0.0; - - /* the nonzero entry L(k,i) */ - Scalar l_ki; - if(DoLDLT) - l_ki = yi / m_diag[i]; - else - yi = l_ki = yi / Lx[Lp[i]]; - - Index p2 = Lp[i] + m_nonZerosPerCol[i]; - Index p; - for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p) - y[Li[p]] -= internal::conj(Lx[p]) * yi; - d -= internal::real(l_ki * internal::conj(yi)); - Li[p] = k; /* store L(k,i) in column form of L */ - Lx[p] = l_ki; - ++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */ - } - if(DoLDLT) - { - m_diag[k] = d; - if(d == RealScalar(0)) - { - ok = false; /* failure, D(k,k) is zero */ - break; - } - } - else - { - Index p = Lp[k] + m_nonZerosPerCol[k]++; - Li[p] = k ; /* store L(k,k) = sqrt (d) in column k */ - if(d <= RealScalar(0)) { - ok = false; /* failure, matrix is not positive definite */ - break; - } - Lx[p] = sqrt(d) ; - } - } - - m_info = ok ? Success : NumericalIssue; - m_factorizationIsOk = true; -} - namespace internal { template<typename Derived, typename Rhs> diff --git a/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h b/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h new file mode 100644 index 000000000..4b249868f --- /dev/null +++ b/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h @@ -0,0 +1,199 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2012 Gael Guennebaud <gael.guennebaud@inria.fr> + +/* + +NOTE: thes functions vave been adapted from the LDL library: + +LDL Copyright (c) 2005 by Timothy A. Davis. All Rights Reserved. + +LDL License: + + Your use or distribution of LDL or any modified version of + LDL implies that you agree to this License. + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 + USA + + Permission is hereby granted to use or copy this program under the + terms of the GNU LGPL, provided that the Copyright, this License, + and the Availability of the original version is retained on all copies. + User documentation of any code that uses this code or any modified + version of this code must cite the Copyright, this License, the + Availability note, and "Used by permission." Permission to modify + the code and to distribute modified code is granted, provided the + Copyright, this License, and the Availability note are retained, + and a notice that the code was modified is included. + */ + +#include "../Core/util/NonMPL2.h" + +#ifndef EIGEN_SIMPLICIAL_CHOLESKY_IMPL_H +#define EIGEN_SIMPLICIAL_CHOLESKY_IMPL_H + +namespace Eigen { + +template<typename Derived> +void SimplicialCholeskyBase<Derived>::analyzePattern_preordered(const CholMatrixType& ap, bool doLDLT) +{ + const Index size = ap.rows(); + m_matrix.resize(size, size); + m_parent.resize(size); + m_nonZerosPerCol.resize(size); + + ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0); + + for(Index k = 0; k < size; ++k) + { + /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */ + m_parent[k] = -1; /* parent of k is not yet known */ + tags[k] = k; /* mark node k as visited */ + m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */ + for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it) + { + Index i = it.index(); + if(i < k) + { + /* follow path from i to root of etree, stop at flagged node */ + for(; tags[i] != k; i = m_parent[i]) + { + /* find parent of i if not yet determined */ + if (m_parent[i] == -1) + m_parent[i] = k; + m_nonZerosPerCol[i]++; /* L (k,i) is nonzero */ + tags[i] = k; /* mark i as visited */ + } + } + } + } + + /* construct Lp index array from m_nonZerosPerCol column counts */ + Index* Lp = m_matrix.outerIndexPtr(); + Lp[0] = 0; + for(Index k = 0; k < size; ++k) + Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1); + + m_matrix.resizeNonZeros(Lp[size]); + + m_isInitialized = true; + m_info = Success; + m_analysisIsOk = true; + m_factorizationIsOk = false; +} + + +template<typename Derived> +template<bool DoLDLT> +void SimplicialCholeskyBase<Derived>::factorize_preordered(const CholMatrixType& ap) +{ + using std::sqrt; + + eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); + eigen_assert(ap.rows()==ap.cols()); + const Index size = ap.rows(); + eigen_assert(m_parent.size()==size); + eigen_assert(m_nonZerosPerCol.size()==size); + + const Index* Lp = m_matrix.outerIndexPtr(); + Index* Li = m_matrix.innerIndexPtr(); + Scalar* Lx = m_matrix.valuePtr(); + + ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0); + ei_declare_aligned_stack_constructed_variable(Index, pattern, size, 0); + ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0); + + bool ok = true; + m_diag.resize(DoLDLT ? size : 0); + + for(Index k = 0; k < size; ++k) + { + // compute nonzero pattern of kth row of L, in topological order + y[k] = 0.0; // Y(0:k) is now all zero + Index top = size; // stack for pattern is empty + tags[k] = k; // mark node k as visited + m_nonZerosPerCol[k] = 0; // count of nonzeros in column k of L + for(typename MatrixType::InnerIterator it(ap,k); it; ++it) + { + Index i = it.index(); + if(i <= k) + { + y[i] += internal::conj(it.value()); /* scatter A(i,k) into Y (sum duplicates) */ + Index len; + for(len = 0; tags[i] != k; i = m_parent[i]) + { + pattern[len++] = i; /* L(k,i) is nonzero */ + tags[i] = k; /* mark i as visited */ + } + while(len > 0) + pattern[--top] = pattern[--len]; + } + } + + /* compute numerical values kth row of L (a sparse triangular solve) */ + + RealScalar d = internal::real(y[k]) * m_shiftScale + m_shiftOffset; // get D(k,k), apply the shift function, and clear Y(k) + y[k] = 0.0; + for(; top < size; ++top) + { + Index i = pattern[top]; /* pattern[top:n-1] is pattern of L(:,k) */ + Scalar yi = y[i]; /* get and clear Y(i) */ + y[i] = 0.0; + + /* the nonzero entry L(k,i) */ + Scalar l_ki; + if(DoLDLT) + l_ki = yi / m_diag[i]; + else + yi = l_ki = yi / Lx[Lp[i]]; + + Index p2 = Lp[i] + m_nonZerosPerCol[i]; + Index p; + for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p) + y[Li[p]] -= internal::conj(Lx[p]) * yi; + d -= internal::real(l_ki * internal::conj(yi)); + Li[p] = k; /* store L(k,i) in column form of L */ + Lx[p] = l_ki; + ++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */ + } + if(DoLDLT) + { + m_diag[k] = d; + if(d == RealScalar(0)) + { + ok = false; /* failure, D(k,k) is zero */ + break; + } + } + else + { + Index p = Lp[k] + m_nonZerosPerCol[k]++; + Li[p] = k ; /* store L(k,k) = sqrt (d) in column k */ + if(d <= RealScalar(0)) { + ok = false; /* failure, matrix is not positive definite */ + break; + } + Lx[p] = sqrt(d) ; + } + } + + m_info = ok ? Success : NumericalIssue; + m_factorizationIsOk = true; +} + +} // end namespace Eigen + +#endif // EIGEN_SIMPLICIAL_CHOLESKY_IMPL_H diff --git a/unsupported/Eigen/IterativeSolvers b/unsupported/Eigen/IterativeSolvers index 04341b22e..aa15403db 100644 --- a/unsupported/Eigen/IterativeSolvers +++ b/unsupported/Eigen/IterativeSolvers @@ -27,8 +27,11 @@ #include "../../Eigen/src/misc/Solve.h" #include "../../Eigen/src/misc/SparseSolve.h" +#ifndef EIGEN_MPL2_ONLY #include "src/IterativeSolvers/IterationController.h" #include "src/IterativeSolvers/ConstrainedConjGrad.h" +#endif + #include "src/IterativeSolvers/IncompleteLU.h" #include "../../Eigen/Jacobi" #include "../../Eigen/Householder" diff --git a/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h b/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h index b83bf7aef..3f18beeeb 100644 --- a/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h +++ b/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h @@ -2,10 +2,6 @@ // 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/. /* NOTE The functions of this file have been adapted from the GMM++ library */ diff --git a/unsupported/Eigen/src/IterativeSolvers/IterationController.h b/unsupported/Eigen/src/IterativeSolvers/IterationController.h index ea86dfef4..c9c1a4be2 100644 --- a/unsupported/Eigen/src/IterativeSolvers/IterationController.h +++ b/unsupported/Eigen/src/IterativeSolvers/IterationController.h @@ -2,10 +2,6 @@ // for linear algebra. // // Copyright (C) 2008-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/. /* NOTE The class IterationController has been adapted from the iteration * class of the GMM++ and ITL libraries. |