1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2012 Gael Guennebaud <gael.guennebaud (at) inria.fr> 5 6 /* 7 8 NOTE: thes functions vave been adapted from the LDL library: 9 10 LDL Copyright (c) 2005 by Timothy A. Davis. All Rights Reserved. 11 12 LDL License: 13 14 Your use or distribution of LDL or any modified version of 15 LDL implies that you agree to this License. 16 17 This library is free software; you can redistribute it and/or 18 modify it under the terms of the GNU Lesser General Public 19 License as published by the Free Software Foundation; either 20 version 2.1 of the License, or (at your option) any later version. 21 22 This library is distributed in the hope that it will be useful, 23 but WITHOUT ANY WARRANTY; without even the implied warranty of 24 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 25 Lesser General Public License for more details. 26 27 You should have received a copy of the GNU Lesser General Public 28 License along with this library; if not, write to the Free Software 29 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 30 USA 31 32 Permission is hereby granted to use or copy this program under the 33 terms of the GNU LGPL, provided that the Copyright, this License, 34 and the Availability of the original version is retained on all copies. 35 User documentation of any code that uses this code or any modified 36 version of this code must cite the Copyright, this License, the 37 Availability note, and "Used by permission." Permission to modify 38 the code and to distribute modified code is granted, provided the 39 Copyright, this License, and the Availability note are retained, 40 and a notice that the code was modified is included. 41 */ 42 43 #include "../Core/util/NonMPL2.h" 44 45 #ifndef EIGEN_SIMPLICIAL_CHOLESKY_IMPL_H 46 #define EIGEN_SIMPLICIAL_CHOLESKY_IMPL_H 47 48 namespace Eigen { 49 50 template<typename Derived> 51 void SimplicialCholeskyBase<Derived>::analyzePattern_preordered(const CholMatrixType& ap, bool doLDLT) 52 { 53 const Index size = ap.rows(); 54 m_matrix.resize(size, size); 55 m_parent.resize(size); 56 m_nonZerosPerCol.resize(size); 57 58 ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0); 59 60 for(Index k = 0; k < size; ++k) 61 { 62 /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */ 63 m_parent[k] = -1; /* parent of k is not yet known */ 64 tags[k] = k; /* mark node k as visited */ 65 m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */ 66 for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it) 67 { 68 Index i = it.index(); 69 if(i < k) 70 { 71 /* follow path from i to root of etree, stop at flagged node */ 72 for(; tags[i] != k; i = m_parent[i]) 73 { 74 /* find parent of i if not yet determined */ 75 if (m_parent[i] == -1) 76 m_parent[i] = k; 77 m_nonZerosPerCol[i]++; /* L (k,i) is nonzero */ 78 tags[i] = k; /* mark i as visited */ 79 } 80 } 81 } 82 } 83 84 /* construct Lp index array from m_nonZerosPerCol column counts */ 85 Index* Lp = m_matrix.outerIndexPtr(); 86 Lp[0] = 0; 87 for(Index k = 0; k < size; ++k) 88 Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1); 89 90 m_matrix.resizeNonZeros(Lp[size]); 91 92 m_isInitialized = true; 93 m_info = Success; 94 m_analysisIsOk = true; 95 m_factorizationIsOk = false; 96 } 97 98 99 template<typename Derived> 100 template<bool DoLDLT> 101 void SimplicialCholeskyBase<Derived>::factorize_preordered(const CholMatrixType& ap) 102 { 103 using std::sqrt; 104 105 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); 106 eigen_assert(ap.rows()==ap.cols()); 107 const Index size = ap.rows(); 108 eigen_assert(m_parent.size()==size); 109 eigen_assert(m_nonZerosPerCol.size()==size); 110 111 const Index* Lp = m_matrix.outerIndexPtr(); 112 Index* Li = m_matrix.innerIndexPtr(); 113 Scalar* Lx = m_matrix.valuePtr(); 114 115 ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0); 116 ei_declare_aligned_stack_constructed_variable(Index, pattern, size, 0); 117 ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0); 118 119 bool ok = true; 120 m_diag.resize(DoLDLT ? size : 0); 121 122 for(Index k = 0; k < size; ++k) 123 { 124 // compute nonzero pattern of kth row of L, in topological order 125 y[k] = 0.0; // Y(0:k) is now all zero 126 Index top = size; // stack for pattern is empty 127 tags[k] = k; // mark node k as visited 128 m_nonZerosPerCol[k] = 0; // count of nonzeros in column k of L 129 for(typename MatrixType::InnerIterator it(ap,k); it; ++it) 130 { 131 Index i = it.index(); 132 if(i <= k) 133 { 134 y[i] += numext::conj(it.value()); /* scatter A(i,k) into Y (sum duplicates) */ 135 Index len; 136 for(len = 0; tags[i] != k; i = m_parent[i]) 137 { 138 pattern[len++] = i; /* L(k,i) is nonzero */ 139 tags[i] = k; /* mark i as visited */ 140 } 141 while(len > 0) 142 pattern[--top] = pattern[--len]; 143 } 144 } 145 146 /* compute numerical values kth row of L (a sparse triangular solve) */ 147 148 RealScalar d = numext::real(y[k]) * m_shiftScale + m_shiftOffset; // get D(k,k), apply the shift function, and clear Y(k) 149 y[k] = 0.0; 150 for(; top < size; ++top) 151 { 152 Index i = pattern[top]; /* pattern[top:n-1] is pattern of L(:,k) */ 153 Scalar yi = y[i]; /* get and clear Y(i) */ 154 y[i] = 0.0; 155 156 /* the nonzero entry L(k,i) */ 157 Scalar l_ki; 158 if(DoLDLT) 159 l_ki = yi / m_diag[i]; 160 else 161 yi = l_ki = yi / Lx[Lp[i]]; 162 163 Index p2 = Lp[i] + m_nonZerosPerCol[i]; 164 Index p; 165 for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p) 166 y[Li[p]] -= numext::conj(Lx[p]) * yi; 167 d -= numext::real(l_ki * numext::conj(yi)); 168 Li[p] = k; /* store L(k,i) in column form of L */ 169 Lx[p] = l_ki; 170 ++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */ 171 } 172 if(DoLDLT) 173 { 174 m_diag[k] = d; 175 if(d == RealScalar(0)) 176 { 177 ok = false; /* failure, D(k,k) is zero */ 178 break; 179 } 180 } 181 else 182 { 183 Index p = Lp[k] + m_nonZerosPerCol[k]++; 184 Li[p] = k ; /* store L(k,k) = sqrt (d) in column k */ 185 if(d <= RealScalar(0)) { 186 ok = false; /* failure, matrix is not positive definite */ 187 break; 188 } 189 Lx[p] = sqrt(d) ; 190 } 191 } 192 193 m_info = ok ? Success : NumericalIssue; 194 m_factorizationIsOk = true; 195 } 196 197 } // end namespace Eigen 198 199 #endif // EIGEN_SIMPLICIAL_CHOLESKY_IMPL_H 200
