| // 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 StorageIndex size = StorageIndex(ap.rows()); |
| m_matrix.resize(size, size); |
| m_parent.resize(size); |
| m_nonZerosPerCol.resize(size); |
| |
| ei_declare_aligned_stack_constructed_variable(StorageIndex, tags, size, 0); |
| |
| for(StorageIndex 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) |
| { |
| StorageIndex 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 */ |
| StorageIndex* Lp = m_matrix.outerIndexPtr(); |
| Lp[0] = 0; |
| for(StorageIndex 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()); |
| eigen_assert(m_parent.size()==ap.rows()); |
| eigen_assert(m_nonZerosPerCol.size()==ap.rows()); |
| |
| const StorageIndex size = StorageIndex(ap.rows()); |
| const StorageIndex* Lp = m_matrix.outerIndexPtr(); |
| StorageIndex* 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(StorageIndex, pattern, size, 0); |
| ei_declare_aligned_stack_constructed_variable(StorageIndex, tags, size, 0); |
| |
| bool ok = true; |
| m_diag.resize(DoLDLT ? size : 0); |
| |
| for(StorageIndex 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 |
| StorageIndex 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 CholMatrixType::InnerIterator it(ap,k); it; ++it) |
| { |
| StorageIndex i = it.index(); |
| if(i <= k) |
| { |
| y[i] += numext::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 = numext::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]] -= numext::conj(Lx[p]) * yi; |
| d -= numext::real(l_ki * numext::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 |