SimplicialCholesky.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 /*
26 
27 NOTE: the _symbolic, and _numeric functions has been adapted from
28  the LDL library:
29 
30 LDL Copyright (c) 2005 by Timothy A. Davis. All Rights Reserved.
31 
32 LDL License:
33 
34  Your use or distribution of LDL or any modified version of
35  LDL implies that you agree to this License.
36 
37  This library is free software; you can redistribute it and/or
38  modify it under the terms of the GNU Lesser General Public
39  License as published by the Free Software Foundation; either
40  version 2.1 of the License, or (at your option) any later version.
41 
42  This library is distributed in the hope that it will be useful,
43  but WITHOUT ANY WARRANTY; without even the implied warranty of
44  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
45  Lesser General Public License for more details.
46 
47  You should have received a copy of the GNU Lesser General Public
48  License along with this library; if not, write to the Free Software
49  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
50  USA
51 
52  Permission is hereby granted to use or copy this program under the
53  terms of the GNU LGPL, provided that the Copyright, this License,
54  and the Availability of the original version is retained on all copies.
55  User documentation of any code that uses this code or any modified
56  version of this code must cite the Copyright, this License, the
57  Availability note, and "Used by permission." Permission to modify
58  the code and to distribute modified code is granted, provided the
59  Copyright, this License, and the Availability note are retained,
60  and a notice that the code was modified is included.
61  */
62 
63 #ifndef EIGEN_SIMPLICIAL_CHOLESKY_H
64 #define EIGEN_SIMPLICIAL_CHOLESKY_H
65 
66 namespace Eigen {
67 
71 };
72 
85 template<typename Derived>
86 class SimplicialCholeskyBase : internal::noncopyable
87 {
88  public:
89  typedef typename internal::traits<Derived>::MatrixType MatrixType;
90  enum { UpLo = internal::traits<Derived>::UpLo };
91  typedef typename MatrixType::Scalar Scalar;
92  typedef typename MatrixType::RealScalar RealScalar;
93  typedef typename MatrixType::Index Index;
96 
97  public:
98 
102  {}
103 
106  {
107  derived().compute(matrix);
108  }
109 
111  {
112  }
113 
114  Derived& derived() { return *static_cast<Derived*>(this); }
115  const Derived& derived() const { return *static_cast<const Derived*>(this); }
116 
117  inline Index cols() const { return m_matrix.cols(); }
118  inline Index rows() const { return m_matrix.rows(); }
119 
126  {
127  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
128  return m_info;
129  }
130 
135  template<typename Rhs>
136  inline const internal::solve_retval<SimplicialCholeskyBase, Rhs>
137  solve(const MatrixBase<Rhs>& b) const
138  {
139  eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
140  eigen_assert(rows()==b.rows()
141  && "SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b");
142  return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
143  }
144 
149  template<typename Rhs>
150  inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>
151  solve(const SparseMatrixBase<Rhs>& b) const
152  {
153  eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
154  eigen_assert(rows()==b.rows()
155  && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
156  return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
157  }
158 
162  { return m_P; }
163 
167  { return m_Pinv; }
168 
178  Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1)
179  {
180  m_shiftOffset = offset;
181  m_shiftScale = scale;
182  return derived();
183  }
184 
185 #ifndef EIGEN_PARSED_BY_DOXYGEN
186 
187  template<typename Stream>
188  void dumpMemory(Stream& s)
189  {
190  int total = 0;
191  s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
192  s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
193  s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
194  s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
195  s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
196  s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
197  s << " TOTAL: " << (total>> 20) << "Mb" << "\n";
198  }
199 
201  template<typename Rhs,typename Dest>
202  void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
203  {
204  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
205  eigen_assert(m_matrix.rows()==b.rows());
206 
207  if(m_info!=Success)
208  return;
209 
210  if(m_P.size()>0)
211  dest = m_Pinv * b;
212  else
213  dest = b;
214 
215  if(m_matrix.nonZeros()>0) // otherwise L==I
216  derived().matrixL().solveInPlace(dest);
217 
218  if(m_diag.size()>0)
219  dest = m_diag.asDiagonal().inverse() * dest;
220 
221  if (m_matrix.nonZeros()>0) // otherwise I==I
222  derived().matrixU().solveInPlace(dest);
223 
224  if(m_P.size()>0)
225  dest = m_P * dest;
226  }
227 
229  template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
230  void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
231  {
232  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
233  eigen_assert(m_matrix.rows()==b.rows());
234 
235  // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
236  static const int NbColsAtOnce = 4;
237  int rhsCols = b.cols();
238  int size = b.rows();
240  for(int k=0; k<rhsCols; k+=NbColsAtOnce)
241  {
242  int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
243  tmp.leftCols(actualCols) = b.middleCols(k,actualCols);
244  tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols));
245  dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView();
246  }
247  }
248 
249 #endif // EIGEN_PARSED_BY_DOXYGEN
250 
251  protected:
252 
254  template<bool DoLDLT>
255  void compute(const MatrixType& matrix)
256  {
257  eigen_assert(matrix.rows()==matrix.cols());
258  Index size = matrix.cols();
259  CholMatrixType ap(size,size);
260  ordering(matrix, ap);
261  analyzePattern_preordered(ap, DoLDLT);
262  factorize_preordered<DoLDLT>(ap);
263  }
264 
265  template<bool DoLDLT>
266  void factorize(const MatrixType& a)
267  {
268  eigen_assert(a.rows()==a.cols());
269  int size = a.cols();
270  CholMatrixType ap(size,size);
271  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv);
272  factorize_preordered<DoLDLT>(ap);
273  }
274 
275  template<bool DoLDLT>
276  void factorize_preordered(const CholMatrixType& a);
277 
278  void analyzePattern(const MatrixType& a, bool doLDLT)
279  {
280  eigen_assert(a.rows()==a.cols());
281  int size = a.cols();
282  CholMatrixType ap(size,size);
283  ordering(a, ap);
284  analyzePattern_preordered(ap,doLDLT);
285  }
286  void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);
287 
288  void ordering(const MatrixType& a, CholMatrixType& ap);
289 
291  struct keep_diag {
292  inline bool operator() (const Index& row, const Index& col, const Scalar&) const
293  {
294  return row!=col;
295  }
296  };
297 
302 
304  VectorType m_diag; // the diagonal coefficients (LDLT mode)
305  VectorXi m_parent; // elimination tree
309 
312 };
313 
314 template<typename _MatrixType, int _UpLo = Lower> class SimplicialLLT;
315 template<typename _MatrixType, int _UpLo = Lower> class SimplicialLDLT;
316 template<typename _MatrixType, int _UpLo = Lower> class SimplicialCholesky;
317 
318 namespace internal {
319 
320 template<typename _MatrixType, int _UpLo> struct traits<SimplicialLLT<_MatrixType,_UpLo> >
321 {
322  typedef _MatrixType MatrixType;
323  enum { UpLo = _UpLo };
324  typedef typename MatrixType::Scalar Scalar;
325  typedef typename MatrixType::Index Index;
326  typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
327  typedef SparseTriangularView<CholMatrixType, Eigen::Lower> MatrixL;
328  typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
329  static inline MatrixL getL(const MatrixType& m) { return m; }
330  static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
331 };
332 
333 template<typename _MatrixType,int _UpLo> struct traits<SimplicialLDLT<_MatrixType,_UpLo> >
334 {
335  typedef _MatrixType MatrixType;
336  enum { UpLo = _UpLo };
337  typedef typename MatrixType::Scalar Scalar;
338  typedef typename MatrixType::Index Index;
339  typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
340  typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower> MatrixL;
341  typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
342  static inline MatrixL getL(const MatrixType& m) { return m; }
343  static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
344 };
345 
346 template<typename _MatrixType, int _UpLo> struct traits<SimplicialCholesky<_MatrixType,_UpLo> >
347 {
348  typedef _MatrixType MatrixType;
349  enum { UpLo = _UpLo };
350 };
351 
352 }
353 
368 template<typename _MatrixType, int _UpLo>
369  class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo> >
370 {
371 public:
372  typedef _MatrixType MatrixType;
373  enum { UpLo = _UpLo };
375  typedef typename MatrixType::Scalar Scalar;
376  typedef typename MatrixType::RealScalar RealScalar;
377  typedef typename MatrixType::Index Index;
380  typedef internal::traits<SimplicialLLT> Traits;
381  typedef typename Traits::MatrixL MatrixL;
382  typedef typename Traits::MatrixU MatrixU;
383 public:
387  SimplicialLLT(const MatrixType& matrix)
388  : Base(matrix) {}
389 
391  inline const MatrixL matrixL() const {
392  eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
393  return Traits::getL(Base::m_matrix);
394  }
395 
397  inline const MatrixU matrixU() const {
398  eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
399  return Traits::getU(Base::m_matrix);
400  }
401 
404  {
405  Base::template compute<false>(matrix);
406  return *this;
407  }
408 
415  void analyzePattern(const MatrixType& a)
416  {
417  Base::analyzePattern(a, false);
418  }
419 
426  void factorize(const MatrixType& a)
427  {
428  Base::template factorize<false>(a);
429  }
430 
433  {
434  Scalar detL = Base::m_matrix.diagonal().prod();
435  return internal::abs2(detL);
436  }
437 };
438 
453 template<typename _MatrixType, int _UpLo>
454  class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo> >
455 {
456 public:
457  typedef _MatrixType MatrixType;
458  enum { UpLo = _UpLo };
460  typedef typename MatrixType::Scalar Scalar;
461  typedef typename MatrixType::RealScalar RealScalar;
462  typedef typename MatrixType::Index Index;
465  typedef internal::traits<SimplicialLDLT> Traits;
466  typedef typename Traits::MatrixL MatrixL;
467  typedef typename Traits::MatrixU MatrixU;
468 public:
471 
473  SimplicialLDLT(const MatrixType& matrix)
474  : Base(matrix) {}
475 
477  inline const VectorType vectorD() const {
478  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
479  return Base::m_diag;
480  }
482  inline const MatrixL matrixL() const {
483  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
484  return Traits::getL(Base::m_matrix);
485  }
486 
488  inline const MatrixU matrixU() const {
489  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
490  return Traits::getU(Base::m_matrix);
491  }
492 
495  {
496  Base::template compute<true>(matrix);
497  return *this;
498  }
499 
506  void analyzePattern(const MatrixType& a)
507  {
508  Base::analyzePattern(a, true);
509  }
510 
517  void factorize(const MatrixType& a)
518  {
519  Base::template factorize<true>(a);
520  }
521 
524  {
525  return Base::m_diag.prod();
526  }
527 };
528 
535 template<typename _MatrixType, int _UpLo>
536  class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo> >
537 {
538 public:
539  typedef _MatrixType MatrixType;
540  enum { UpLo = _UpLo };
542  typedef typename MatrixType::Scalar Scalar;
543  typedef typename MatrixType::RealScalar RealScalar;
544  typedef typename MatrixType::Index Index;
547  typedef internal::traits<SimplicialCholesky> Traits;
548  typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
549  typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
550  public:
552 
554  : Base(), m_LDLT(true)
555  {
556  compute(matrix);
557  }
558 
560  {
561  switch(mode)
562  {
564  m_LDLT = false;
565  break;
567  m_LDLT = true;
568  break;
569  default:
570  break;
571  }
572 
573  return *this;
574  }
575 
576  inline const VectorType vectorD() const {
577  eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
578  return Base::m_diag;
579  }
580  inline const CholMatrixType rawMatrix() const {
581  eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
582  return Base::m_matrix;
583  }
584 
587  {
588  if(m_LDLT)
589  Base::template compute<true>(matrix);
590  else
591  Base::template compute<false>(matrix);
592  return *this;
593  }
594 
601  void analyzePattern(const MatrixType& a)
602  {
604  }
605 
612  void factorize(const MatrixType& a)
613  {
614  if(m_LDLT)
615  Base::template factorize<true>(a);
616  else
617  Base::template factorize<false>(a);
618  }
619 
621  template<typename Rhs,typename Dest>
622  void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
623  {
624  eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
625  eigen_assert(Base::m_matrix.rows()==b.rows());
626 
627  if(Base::m_info!=Success)
628  return;
629 
630  if(Base::m_P.size()>0)
631  dest = Base::m_Pinv * b;
632  else
633  dest = b;
634 
635  if(Base::m_matrix.nonZeros()>0) // otherwise L==I
636  {
637  if(m_LDLT)
638  LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
639  else
640  LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
641  }
642 
643  if(Base::m_diag.size()>0)
644  dest = Base::m_diag.asDiagonal().inverse() * dest;
645 
646  if (Base::m_matrix.nonZeros()>0) // otherwise I==I
647  {
648  if(m_LDLT)
649  LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
650  else
651  LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
652  }
653 
654  if(Base::m_P.size()>0)
655  dest = Base::m_P * dest;
656  }
657 
659  {
660  if(m_LDLT)
661  {
662  return Base::m_diag.prod();
663  }
664  else
665  {
667  return internal::abs2(detL);
668  }
669  }
670 
671  protected:
672  bool m_LDLT;
673 };
674 
675 template<typename Derived>
677 {
678  eigen_assert(a.rows()==a.cols());
679  const Index size = a.rows();
680  // TODO allows to configure the permutation
681  {
682  CholMatrixType C;
683  C = a.template selfadjointView<UpLo>();
684  // remove diagonal entries:
685  // seems not to be needed
686  // C.prune(keep_diag());
688  }
689 
690  if(m_P.size()>0)
691  m_Pinv = m_P.inverse();
692  else
693  m_Pinv.resize(0);
694 
695  ap.resize(size,size);
696  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv);
697 }
698 
699 template<typename Derived>
701 {
702  const Index size = ap.rows();
703  m_matrix.resize(size, size);
704  m_parent.resize(size);
705  m_nonZerosPerCol.resize(size);
706 
708 
709  for(Index k = 0; k < size; ++k)
710  {
711  /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
712  m_parent[k] = -1; /* parent of k is not yet known */
713  tags[k] = k; /* mark node k as visited */
714  m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */
715  for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
716  {
717  Index i = it.index();
718  if(i < k)
719  {
720  /* follow path from i to root of etree, stop at flagged node */
721  for(; tags[i] != k; i = m_parent[i])
722  {
723  /* find parent of i if not yet determined */
724  if (m_parent[i] == -1)
725  m_parent[i] = k;
726  m_nonZerosPerCol[i]++; /* L (k,i) is nonzero */
727  tags[i] = k; /* mark i as visited */
728  }
729  }
730  }
731  }
732 
733  /* construct Lp index array from m_nonZerosPerCol column counts */
734  Index* Lp = m_matrix.outerIndexPtr();
735  Lp[0] = 0;
736  for(Index k = 0; k < size; ++k)
737  Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1);
738 
739  m_matrix.resizeNonZeros(Lp[size]);
740 
741  m_isInitialized = true;
742  m_info = Success;
743  m_analysisIsOk = true;
744  m_factorizationIsOk = false;
745 }
746 
747 
748 template<typename Derived>
749 template<bool DoLDLT>
751 {
752  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
753  eigen_assert(ap.rows()==ap.cols());
754  const Index size = ap.rows();
755  eigen_assert(m_parent.size()==size);
756  eigen_assert(m_nonZerosPerCol.size()==size);
757 
758  const Index* Lp = m_matrix.outerIndexPtr();
759  Index* Li = m_matrix.innerIndexPtr();
760  Scalar* Lx = m_matrix.valuePtr();
761 
765 
766  bool ok = true;
767  m_diag.resize(DoLDLT ? size : 0);
768 
769  for(Index k = 0; k < size; ++k)
770  {
771  // compute nonzero pattern of kth row of L, in topological order
772  y[k] = 0.0; // Y(0:k) is now all zero
773  Index top = size; // stack for pattern is empty
774  tags[k] = k; // mark node k as visited
775  m_nonZerosPerCol[k] = 0; // count of nonzeros in column k of L
776  for(typename MatrixType::InnerIterator it(ap,k); it; ++it)
777  {
778  Index i = it.index();
779  if(i <= k)
780  {
781  y[i] += internal::conj(it.value()); /* scatter A(i,k) into Y (sum duplicates) */
782  Index len;
783  for(len = 0; tags[i] != k; i = m_parent[i])
784  {
785  pattern[len++] = i; /* L(k,i) is nonzero */
786  tags[i] = k; /* mark i as visited */
787  }
788  while(len > 0)
789  pattern[--top] = pattern[--len];
790  }
791  }
792 
793  /* compute numerical values kth row of L (a sparse triangular solve) */
794 
795  RealScalar d = internal::real(y[k]) * m_shiftScale + m_shiftOffset; // get D(k,k), apply the shift function, and clear Y(k)
796  y[k] = 0.0;
797  for(; top < size; ++top)
798  {
799  Index i = pattern[top]; /* pattern[top:n-1] is pattern of L(:,k) */
800  Scalar yi = y[i]; /* get and clear Y(i) */
801  y[i] = 0.0;
802 
803  /* the nonzero entry L(k,i) */
804  Scalar l_ki;
805  if(DoLDLT)
806  l_ki = yi / m_diag[i];
807  else
808  yi = l_ki = yi / Lx[Lp[i]];
809 
810  Index p2 = Lp[i] + m_nonZerosPerCol[i];
811  Index p;
812  for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p)
813  y[Li[p]] -= internal::conj(Lx[p]) * yi;
814  d -= internal::real(l_ki * internal::conj(yi));
815  Li[p] = k; /* store L(k,i) in column form of L */
816  Lx[p] = l_ki;
817  ++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */
818  }
819  if(DoLDLT)
820  {
821  m_diag[k] = d;
822  if(d == RealScalar(0))
823  {
824  ok = false; /* failure, D(k,k) is zero */
825  break;
826  }
827  }
828  else
829  {
830  Index p = Lp[k] + m_nonZerosPerCol[k]++;
831  Li[p] = k ; /* store L(k,k) = sqrt (d) in column k */
832  if(d <= RealScalar(0)) {
833  ok = false; /* failure, matrix is not positive definite */
834  break;
835  }
836  Lx[p] = internal::sqrt(d) ;
837  }
838  }
839 
840  m_info = ok ? Success : NumericalIssue;
841  m_factorizationIsOk = true;
842 }
843 
844 namespace internal {
845 
846 template<typename Derived, typename Rhs>
847 struct solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
848  : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
849 {
851  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
852 
853  template<typename Dest> void evalTo(Dest& dst) const
854  {
855  dec().derived()._solve(rhs(),dst);
856  }
857 };
858 
859 template<typename Derived, typename Rhs>
860 struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
861  : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
862 {
863  typedef SimplicialCholeskyBase<Derived> Dec;
865 
866  template<typename Dest> void evalTo(Dest& dst) const
867  {
868  dec().derived()._solve_sparse(rhs(),dst);
869  }
870 };
871 
872 } // end namespace internal
873 
874 } // end namespace Eigen
875 
876 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H