[136] | 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@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
|
---|