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 
88 template<typename Derived>
89 class SimplicialCholeskyBase : internal::noncopyable
90 {
91  public:
92  typedef typename internal::traits<Derived>::MatrixType MatrixType;
93  enum { UpLo = internal::traits<Derived>::UpLo };
94  typedef typename MatrixType::Scalar Scalar;
95  typedef typename MatrixType::RealScalar RealScalar;
96  typedef typename MatrixType::Index Index;
99 
100  public:
101 
105  {}
106 
109  {
110  derived().compute(matrix);
111  }
112 
114  {
115  }
116 
117  Derived& derived() { return *static_cast<Derived*>(this); }
118  const Derived& derived() const { return *static_cast<const Derived*>(this); }
119 
120  inline Index cols() const { return m_matrix.cols(); }
121  inline Index rows() const { return m_matrix.rows(); }
122 
129  {
130  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
131  return m_info;
132  }
133 
138  template<typename Rhs>
139  inline const internal::solve_retval<SimplicialCholeskyBase, Rhs>
140  solve(const MatrixBase<Rhs>& b) const
141  {
142  eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
143  eigen_assert(rows()==b.rows()
144  && "SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b");
145  return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
146  }
147 
152  template<typename Rhs>
153  inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>
154  solve(const SparseMatrixBase<Rhs>& b) const
155  {
156  eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
157  eigen_assert(rows()==b.rows()
158  && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
159  return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
160  }
161 
165  { return m_P; }
166 
170  { return m_Pinv; }
171 
181  Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1)
182  {
183  m_shiftOffset = offset;
184  m_shiftScale = scale;
185  return derived();
186  }
187 
188 #ifndef EIGEN_PARSED_BY_DOXYGEN
189 
190  template<typename Stream>
191  void dumpMemory(Stream& s)
192  {
193  int total = 0;
194  s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
195  s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
196  s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
197  s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
198  s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
199  s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
200  s << " TOTAL: " << (total>> 20) << "Mb" << "\n";
201  }
202 
204  template<typename Rhs,typename Dest>
205  void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
206  {
207  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
208  eigen_assert(m_matrix.rows()==b.rows());
209 
210  if(m_info!=Success)
211  return;
212 
213  if(m_P.size()>0)
214  dest = m_P * b;
215  else
216  dest = b;
217 
218  if(m_matrix.nonZeros()>0) // otherwise L==I
219  derived().matrixL().solveInPlace(dest);
220 
221  if(m_diag.size()>0)
222  dest = m_diag.asDiagonal().inverse() * dest;
223 
224  if (m_matrix.nonZeros()>0) // otherwise U==I
225  derived().matrixU().solveInPlace(dest);
226 
227  if(m_P.size()>0)
228  dest = m_Pinv * dest;
229  }
230 
232  template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
233  void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
234  {
235  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
236  eigen_assert(m_matrix.rows()==b.rows());
237 
238  // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
239  static const int NbColsAtOnce = 4;
240  int rhsCols = b.cols();
241  int size = b.rows();
243  for(int k=0; k<rhsCols; k+=NbColsAtOnce)
244  {
245  int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
246  tmp.leftCols(actualCols) = b.middleCols(k,actualCols);
247  tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols));
248  dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView();
249  }
250  }
251 
252 #endif // EIGEN_PARSED_BY_DOXYGEN
253 
254  protected:
255 
257  template<bool DoLDLT>
258  void compute(const MatrixType& matrix)
259  {
260  eigen_assert(matrix.rows()==matrix.cols());
261  Index size = matrix.cols();
262  CholMatrixType ap(size,size);
263  ordering(matrix, ap);
264  analyzePattern_preordered(ap, DoLDLT);
265  factorize_preordered<DoLDLT>(ap);
266  }
267 
268  template<bool DoLDLT>
269  void factorize(const MatrixType& a)
270  {
271  eigen_assert(a.rows()==a.cols());
272  int size = a.cols();
273  CholMatrixType ap(size,size);
274  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
275  factorize_preordered<DoLDLT>(ap);
276  }
277 
278  template<bool DoLDLT>
279  void factorize_preordered(const CholMatrixType& a);
280 
281  void analyzePattern(const MatrixType& a, bool doLDLT)
282  {
283  eigen_assert(a.rows()==a.cols());
284  int size = a.cols();
285  CholMatrixType ap(size,size);
286  ordering(a, ap);
287  analyzePattern_preordered(ap,doLDLT);
288  }
289  void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);
290 
291  void ordering(const MatrixType& a, CholMatrixType& ap);
292 
294  struct keep_diag {
295  inline bool operator() (const Index& row, const Index& col, const Scalar&) const
296  {
297  return row!=col;
298  }
299  };
300 
305 
307  VectorType m_diag; // the diagonal coefficients (LDLT mode)
308  VectorXi m_parent; // elimination tree
312 
315 };
316 
317 template<typename _MatrixType, int _UpLo = Lower> class SimplicialLLT;
318 template<typename _MatrixType, int _UpLo = Lower> class SimplicialLDLT;
319 template<typename _MatrixType, int _UpLo = Lower> class SimplicialCholesky;
320 
321 namespace internal {
322 
323 template<typename _MatrixType, int _UpLo> struct traits<SimplicialLLT<_MatrixType,_UpLo> >
324 {
325  typedef _MatrixType MatrixType;
326  enum { UpLo = _UpLo };
327  typedef typename MatrixType::Scalar Scalar;
328  typedef typename MatrixType::Index Index;
329  typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
330  typedef SparseTriangularView<CholMatrixType, Eigen::Lower> MatrixL;
331  typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
332  static inline MatrixL getL(const MatrixType& m) { return m; }
333  static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
334 };
335 
336 template<typename _MatrixType,int _UpLo> struct traits<SimplicialLDLT<_MatrixType,_UpLo> >
337 {
338  typedef _MatrixType MatrixType;
339  enum { UpLo = _UpLo };
340  typedef typename MatrixType::Scalar Scalar;
341  typedef typename MatrixType::Index Index;
342  typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
343  typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower> MatrixL;
344  typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
345  static inline MatrixL getL(const MatrixType& m) { return m; }
346  static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
347 };
348 
349 template<typename _MatrixType, int _UpLo> struct traits<SimplicialCholesky<_MatrixType,_UpLo> >
350 {
351  typedef _MatrixType MatrixType;
352  enum { UpLo = _UpLo };
353 };
354 
355 }
356 
374 template<typename _MatrixType, int _UpLo>
375  class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo> >
376 {
377 public:
378  typedef _MatrixType MatrixType;
379  enum { UpLo = _UpLo };
381  typedef typename MatrixType::Scalar Scalar;
382  typedef typename MatrixType::RealScalar RealScalar;
383  typedef typename MatrixType::Index Index;
386  typedef internal::traits<SimplicialLLT> Traits;
387  typedef typename Traits::MatrixL MatrixL;
388  typedef typename Traits::MatrixU MatrixU;
389 public:
393  SimplicialLLT(const MatrixType& matrix)
394  : Base(matrix) {}
395 
397  inline const MatrixL matrixL() const {
398  eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
399  return Traits::getL(Base::m_matrix);
400  }
401 
403  inline const MatrixU matrixU() const {
404  eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
405  return Traits::getU(Base::m_matrix);
406  }
407 
410  {
411  Base::template compute<false>(matrix);
412  return *this;
413  }
414 
421  void analyzePattern(const MatrixType& a)
422  {
423  Base::analyzePattern(a, false);
424  }
425 
432  void factorize(const MatrixType& a)
433  {
434  Base::template factorize<false>(a);
435  }
436 
439  {
440  Scalar detL = Base::m_matrix.diagonal().prod();
441  return internal::abs2(detL);
442  }
443 };
444 
462 template<typename _MatrixType, int _UpLo>
463  class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo> >
464 {
465 public:
466  typedef _MatrixType MatrixType;
467  enum { UpLo = _UpLo };
469  typedef typename MatrixType::Scalar Scalar;
470  typedef typename MatrixType::RealScalar RealScalar;
471  typedef typename MatrixType::Index Index;
474  typedef internal::traits<SimplicialLDLT> Traits;
475  typedef typename Traits::MatrixL MatrixL;
476  typedef typename Traits::MatrixU MatrixU;
477 public:
480 
482  SimplicialLDLT(const MatrixType& matrix)
483  : Base(matrix) {}
484 
486  inline const VectorType vectorD() const {
487  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
488  return Base::m_diag;
489  }
491  inline const MatrixL matrixL() const {
492  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
493  return Traits::getL(Base::m_matrix);
494  }
495 
497  inline const MatrixU matrixU() const {
498  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
499  return Traits::getU(Base::m_matrix);
500  }
501 
504  {
505  Base::template compute<true>(matrix);
506  return *this;
507  }
508 
515  void analyzePattern(const MatrixType& a)
516  {
517  Base::analyzePattern(a, true);
518  }
519 
526  void factorize(const MatrixType& a)
527  {
528  Base::template factorize<true>(a);
529  }
530 
533  {
534  return Base::m_diag.prod();
535  }
536 };
537 
544 template<typename _MatrixType, int _UpLo>
545  class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo> >
546 {
547 public:
548  typedef _MatrixType MatrixType;
549  enum { UpLo = _UpLo };
551  typedef typename MatrixType::Scalar Scalar;
552  typedef typename MatrixType::RealScalar RealScalar;
553  typedef typename MatrixType::Index Index;
556  typedef internal::traits<SimplicialCholesky> Traits;
557  typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
558  typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
559  public:
561 
563  : Base(), m_LDLT(true)
564  {
565  compute(matrix);
566  }
567 
569  {
570  switch(mode)
571  {
573  m_LDLT = false;
574  break;
576  m_LDLT = true;
577  break;
578  default:
579  break;
580  }
581 
582  return *this;
583  }
584 
585  inline const VectorType vectorD() const {
586  eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
587  return Base::m_diag;
588  }
589  inline const CholMatrixType rawMatrix() const {
590  eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
591  return Base::m_matrix;
592  }
593 
596  {
597  if(m_LDLT)
598  Base::template compute<true>(matrix);
599  else
600  Base::template compute<false>(matrix);
601  return *this;
602  }
603 
610  void analyzePattern(const MatrixType& a)
611  {
613  }
614 
621  void factorize(const MatrixType& a)
622  {
623  if(m_LDLT)
624  Base::template factorize<true>(a);
625  else
626  Base::template factorize<false>(a);
627  }
628 
630  template<typename Rhs,typename Dest>
631  void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
632  {
633  eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
634  eigen_assert(Base::m_matrix.rows()==b.rows());
635 
636  if(Base::m_info!=Success)
637  return;
638 
639  if(Base::m_P.size()>0)
640  dest = Base::m_P * b;
641  else
642  dest = b;
643 
644  if(Base::m_matrix.nonZeros()>0) // otherwise L==I
645  {
646  if(m_LDLT)
647  LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
648  else
649  LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
650  }
651 
652  if(Base::m_diag.size()>0)
653  dest = Base::m_diag.asDiagonal().inverse() * dest;
654 
655  if (Base::m_matrix.nonZeros()>0) // otherwise I==I
656  {
657  if(m_LDLT)
658  LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
659  else
660  LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
661  }
662 
663  if(Base::m_P.size()>0)
664  dest = Base::m_Pinv * dest;
665  }
666 
668  {
669  if(m_LDLT)
670  {
671  return Base::m_diag.prod();
672  }
673  else
674  {
676  return internal::abs2(detL);
677  }
678  }
679 
680  protected:
681  bool m_LDLT;
682 };
683 
684 template<typename Derived>
686 {
687  eigen_assert(a.rows()==a.cols());
688  const Index size = a.rows();
689  // TODO allows to configure the permutation
690  // Note that amd compute the inverse permutation
691  {
692  CholMatrixType C;
693  C = a.template selfadjointView<UpLo>();
694  // remove diagonal entries:
695  // seems not to be needed
696  // C.prune(keep_diag());
698  }
699 
700  if(m_Pinv.size()>0)
701  m_P = m_Pinv.inverse();
702  else
703  m_P.resize(0);
704 
705  ap.resize(size,size);
706  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
707 }
708 
709 template<typename Derived>
711 {
712  const Index size = ap.rows();
713  m_matrix.resize(size, size);
714  m_parent.resize(size);
715  m_nonZerosPerCol.resize(size);
716 
718 
719  for(Index k = 0; k < size; ++k)
720  {
721  /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
722  m_parent[k] = -1; /* parent of k is not yet known */
723  tags[k] = k; /* mark node k as visited */
724  m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */
725  for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
726  {
727  Index i = it.index();
728  if(i < k)
729  {
730  /* follow path from i to root of etree, stop at flagged node */
731  for(; tags[i] != k; i = m_parent[i])
732  {
733  /* find parent of i if not yet determined */
734  if (m_parent[i] == -1)
735  m_parent[i] = k;
736  m_nonZerosPerCol[i]++; /* L (k,i) is nonzero */
737  tags[i] = k; /* mark i as visited */
738  }
739  }
740  }
741  }
742 
743  /* construct Lp index array from m_nonZerosPerCol column counts */
744  Index* Lp = m_matrix.outerIndexPtr();
745  Lp[0] = 0;
746  for(Index k = 0; k < size; ++k)
747  Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1);
748 
749  m_matrix.resizeNonZeros(Lp[size]);
750 
751  m_isInitialized = true;
752  m_info = Success;
753  m_analysisIsOk = true;
754  m_factorizationIsOk = false;
755 }
756 
757 
758 template<typename Derived>
759 template<bool DoLDLT>
761 {
762  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
763  eigen_assert(ap.rows()==ap.cols());
764  const Index size = ap.rows();
765  eigen_assert(m_parent.size()==size);
766  eigen_assert(m_nonZerosPerCol.size()==size);
767 
768  const Index* Lp = m_matrix.outerIndexPtr();
769  Index* Li = m_matrix.innerIndexPtr();
770  Scalar* Lx = m_matrix.valuePtr();
771 
775 
776  bool ok = true;
777  m_diag.resize(DoLDLT ? size : 0);
778 
779  for(Index k = 0; k < size; ++k)
780  {
781  // compute nonzero pattern of kth row of L, in topological order
782  y[k] = 0.0; // Y(0:k) is now all zero
783  Index top = size; // stack for pattern is empty
784  tags[k] = k; // mark node k as visited
785  m_nonZerosPerCol[k] = 0; // count of nonzeros in column k of L
786  for(typename MatrixType::InnerIterator it(ap,k); it; ++it)
787  {
788  Index i = it.index();
789  if(i <= k)
790  {
791  y[i] += internal::conj(it.value()); /* scatter A(i,k) into Y (sum duplicates) */
792  Index len;
793  for(len = 0; tags[i] != k; i = m_parent[i])
794  {
795  pattern[len++] = i; /* L(k,i) is nonzero */
796  tags[i] = k; /* mark i as visited */
797  }
798  while(len > 0)
799  pattern[--top] = pattern[--len];
800  }
801  }
802 
803  /* compute numerical values kth row of L (a sparse triangular solve) */
804 
805  RealScalar d = internal::real(y[k]) * m_shiftScale + m_shiftOffset; // get D(k,k), apply the shift function, and clear Y(k)
806  y[k] = 0.0;
807  for(; top < size; ++top)
808  {
809  Index i = pattern[top]; /* pattern[top:n-1] is pattern of L(:,k) */
810  Scalar yi = y[i]; /* get and clear Y(i) */
811  y[i] = 0.0;
812 
813  /* the nonzero entry L(k,i) */
814  Scalar l_ki;
815  if(DoLDLT)
816  l_ki = yi / m_diag[i];
817  else
818  yi = l_ki = yi / Lx[Lp[i]];
819 
820  Index p2 = Lp[i] + m_nonZerosPerCol[i];
821  Index p;
822  for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p)
823  y[Li[p]] -= internal::conj(Lx[p]) * yi;
824  d -= internal::real(l_ki * internal::conj(yi));
825  Li[p] = k; /* store L(k,i) in column form of L */
826  Lx[p] = l_ki;
827  ++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */
828  }
829  if(DoLDLT)
830  {
831  m_diag[k] = d;
832  if(d == RealScalar(0))
833  {
834  ok = false; /* failure, D(k,k) is zero */
835  break;
836  }
837  }
838  else
839  {
840  Index p = Lp[k] + m_nonZerosPerCol[k]++;
841  Li[p] = k ; /* store L(k,k) = sqrt (d) in column k */
842  if(d <= RealScalar(0)) {
843  ok = false; /* failure, matrix is not positive definite */
844  break;
845  }
846  Lx[p] = internal::sqrt(d) ;
847  }
848  }
849 
850  m_info = ok ? Success : NumericalIssue;
851  m_factorizationIsOk = true;
852 }
853 
854 namespace internal {
855 
856 template<typename Derived, typename Rhs>
857 struct solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
858  : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
859 {
861  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
862 
863  template<typename Dest> void evalTo(Dest& dst) const
864  {
865  dec().derived()._solve(rhs(),dst);
866  }
867 };
868 
869 template<typename Derived, typename Rhs>
870 struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
871  : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
872 {
873  typedef SimplicialCholeskyBase<Derived> Dec;
875 
876  template<typename Dest> void evalTo(Dest& dst) const
877  {
878  dec().derived()._solve_sparse(rhs(),dst);
879  }
880 };
881 
882 } // end namespace internal
883 
884 } // end namespace Eigen
885 
886 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H