Home | History | Annotate | Download | only in SparseCholesky
      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 StorageIndex size = StorageIndex(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(StorageIndex, tags, size, 0);
     59 
     60   for(StorageIndex 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       StorageIndex 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   StorageIndex* Lp = m_matrix.outerIndexPtr();
     86   Lp[0] = 0;
     87   for(StorageIndex 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   eigen_assert(m_parent.size()==ap.rows());
    108   eigen_assert(m_nonZerosPerCol.size()==ap.rows());
    109 
    110   const StorageIndex size = StorageIndex(ap.rows());
    111   const StorageIndex* Lp = m_matrix.outerIndexPtr();
    112   StorageIndex* 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(StorageIndex,  pattern, size, 0);
    117   ei_declare_aligned_stack_constructed_variable(StorageIndex,  tags, size, 0);
    118 
    119   bool ok = true;
    120   m_diag.resize(DoLDLT ? size : 0);
    121 
    122   for(StorageIndex 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     StorageIndex 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 CholMatrixType::InnerIterator it(ap,k); it; ++it)
    130     {
    131       StorageIndex 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