← Back to branch summary

~ubuntu-branches/ubuntu/trusty/blender/trusty

« back to all changes in this revision

Viewing changes to extern/Eigen3/Eigen/src/Cholesky/LDLT.h

  • Committer: Package Import Robot
  • Author(s): Jeremy Bicha
  • Date: 2013-03-06 12:08:47 UTC
  • mfrom: (1.5.1) (14.1.8 experimental)
  • Revision ID: package-import@ubuntu.com-20130306120847-frjfaryb2zrotwcg
Tags: 2.66a-1ubuntu1
https://launchpad.net/bugs/1076930
https://launchpad.net/bugs/1089256
https://launchpad.net/bugs/1052743
https://launchpad.net/bugs/999024
* Resynchronize with Debian (LP: #1076930, #1089256, #1052743, #999024,
  #1122888, #1147084)
* debian/control:
  - Lower build-depends on libavcodec-dev since we're not
    doing the libav9 transition in Ubuntu yet

Show diffs side-by-side

added added

removed removed

Lines of Context:
extern/Eigen3/Eigen/src/Cholesky/LDLT.h
1
1
// This file is part of Eigen, a lightweight C++ template library
2
2
// for linear algebra.
3
3
//
4
 
// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
 
4
// Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
5
5
// Copyright (C) 2009 Keir Mierle <mierle@gmail.com>
6
6
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
7
 
//
8
 
// Eigen is free software; you can redistribute it and/or
9
 
// modify it under the terms of the GNU Lesser General Public
10
 
// License as published by the Free Software Foundation; either
11
 
// version 3 of the License, or (at your option) any later version.
12
 
//
13
 
// Alternatively, you can redistribute it and/or
14
 
// modify it under the terms of the GNU General Public License as
15
 
// published by the Free Software Foundation; either version 2 of
16
 
// the License, or (at your option) any later version.
17
 
//
18
 
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
19
 
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20
 
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
21
 
// GNU General Public License for more details.
22
 
//
23
 
// You should have received a copy of the GNU Lesser General Public
24
 
// License and a copy of the GNU General Public License along with
25
 
// Eigen. If not, see <http://www.gnu.org/licenses/>.
 
7
// Copyright (C) 2011 Timothy E. Holy <tim.holy@gmail.com >
 
8
//
 
9
// This Source Code Form is subject to the terms of the Mozilla
 
10
// Public License v. 2.0. If a copy of the MPL was not distributed
 
11
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
26
12
 
27
13
#ifndef EIGEN_LDLT_H
28
14
#define EIGEN_LDLT_H
29
15
 
 
16
namespace Eigen { 
 
17
 
30
18
namespace internal {
31
19
template<typename MatrixType, int UpLo> struct LDLT_Traits;
32
20
}
33
21
 
34
 
/** \ingroup cholesky_Module
 
22
/** \ingroup Cholesky_Module
35
23
  *
36
24
  * \class LDLT
37
25
  *
38
26
  * \brief Robust Cholesky decomposition of a matrix with pivoting
39
27
  *
40
28
  * \param MatrixType the type of the matrix of which to compute the LDL^T Cholesky decomposition
 
29
  * \param UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper.
 
30
  *             The other triangular part won't be read.
41
31
  *
42
32
  * Perform a robust Cholesky decomposition of a positive semidefinite or negative semidefinite
43
33
  * matrix \f$ A \f$ such that \f$ A =  P^TLDL^*P \f$, where P is a permutation matrix, L
48
38
  * on D also stabilizes the computation.
49
39
  *
50
40
  * Remember that Cholesky decompositions are not rank-revealing. Also, do not use a Cholesky
51
 
        * decomposition to determine whether a system of equations has a solution.
 
41
  * decomposition to determine whether a system of equations has a solution.
52
42
  *
53
43
  * \sa MatrixBase::ldlt(), class LLT
54
44
  */
55
 
 /* THIS PART OF THE DOX IS CURRENTLY DISABLED BECAUSE INACCURATE BECAUSE OF BUG IN THE DECOMPOSITION CODE
56
 
  * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
57
 
  * the strict lower part does not have to store correct values.
58
 
  */
59
45
template<typename _MatrixType, int _UpLo> class LDLT
60
46
{
61
47
  public:
98
84
        m_isInitialized(false)
99
85
    {}
100
86
 
 
87
    /** \brief Constructor with decomposition
 
88
      *
 
89
      * This calculates the decomposition for the input \a matrix.
 
90
      * \sa LDLT(Index size)
 
91
      */
101
92
    LDLT(const MatrixType& matrix)
102
93
      : m_matrix(matrix.rows(), matrix.cols()),
103
94
        m_transpositions(matrix.rows()),
107
98
      compute(matrix);
108
99
    }
109
100
 
 
101
    /** Clear any existing decomposition
 
102
     * \sa rankUpdate(w,sigma)
 
103
     */
 
104
    void setZero()
 
105
    {
 
106
      m_isInitialized = false;
 
107
    }
 
108
 
110
109
    /** \returns a view of the upper triangular matrix U */
111
110
    inline typename Traits::MatrixU matrixU() const
112
111
    {
130
129
    }
131
130
 
132
131
    /** \returns the coefficients of the diagonal matrix D */
133
 
    inline Diagonal<const MatrixType> vectorD(void) const
 
132
    inline Diagonal<const MatrixType> vectorD() const
134
133
    {
135
134
      eigen_assert(m_isInitialized && "LDLT is not initialized.");
136
135
      return m_matrix.diagonal();
137
136
    }
138
137
 
139
138
    /** \returns true if the matrix is positive (semidefinite) */
140
 
    inline bool isPositive(void) const
 
139
    inline bool isPositive() const
141
140
    {
142
141
      eigen_assert(m_isInitialized && "LDLT is not initialized.");
143
142
      return m_sign == 1;
196
195
 
197
196
    LDLT& compute(const MatrixType& matrix);
198
197
 
 
198
    template <typename Derived>
 
199
    LDLT& rankUpdate(const MatrixBase<Derived>& w,RealScalar alpha=1);
 
200
 
199
201
    /** \returns the internal LDLT decomposition matrix
200
202
      *
201
203
      * TODO: document the storage layout
211
213
    inline Index rows() const { return m_matrix.rows(); }
212
214
    inline Index cols() const { return m_matrix.cols(); }
213
215
 
 
216
    /** \brief Reports whether previous computation was successful.
 
217
      *
 
218
      * \returns \c Success if computation was succesful,
 
219
      *          \c NumericalIssue if the matrix.appears to be negative.
 
220
      */
 
221
    ComputationInfo info() const
 
222
    {
 
223
      eigen_assert(m_isInitialized && "LDLT is not initialized.");
 
224
      return Success;
 
225
    }
 
226
 
214
227
  protected:
215
228
 
216
229
    /** \internal
249
262
      return true;
250
263
    }
251
264
 
252
 
    RealScalar cutoff = 0, biggest_in_corner;
 
265
    RealScalar cutoff(0), biggest_in_corner;
253
266
 
254
267
    for (Index k = 0; k < size; ++k)
255
268
    {
317
330
 
318
331
    return true;
319
332
  }
 
333
 
 
334
  // Reference for the algorithm: Davis and Hager, "Multiple Rank
 
335
  // Modifications of a Sparse Cholesky Factorization" (Algorithm 1)
 
336
  // Trivial rearrangements of their computations (Timothy E. Holy)
 
337
  // allow their algorithm to work for rank-1 updates even if the
 
338
  // original matrix is not of full rank.
 
339
  // Here only rank-1 updates are implemented, to reduce the
 
340
  // requirement for intermediate storage and improve accuracy
 
341
  template<typename MatrixType, typename WDerived>
 
342
  static bool updateInPlace(MatrixType& mat, MatrixBase<WDerived>& w, typename MatrixType::RealScalar sigma=1)
 
343
  {
 
344
    using internal::isfinite;
 
345
    typedef typename MatrixType::Scalar Scalar;
 
346
    typedef typename MatrixType::RealScalar RealScalar;
 
347
    typedef typename MatrixType::Index Index;
 
348
 
 
349
    const Index size = mat.rows();
 
350
    eigen_assert(mat.cols() == size && w.size()==size);
 
351
 
 
352
    RealScalar alpha = 1;
 
353
 
 
354
    // Apply the update
 
355
    for (Index j = 0; j < size; j++)
 
356
    {
 
357
      // Check for termination due to an original decomposition of low-rank
 
358
      if (!(isfinite)(alpha))
 
359
        break;
 
360
 
 
361
      // Update the diagonal terms
 
362
      RealScalar dj = real(mat.coeff(j,j));
 
363
      Scalar wj = w.coeff(j);
 
364
      RealScalar swj2 = sigma*abs2(wj);
 
365
      RealScalar gamma = dj*alpha + swj2;
 
366
 
 
367
      mat.coeffRef(j,j) += swj2/alpha;
 
368
      alpha += swj2/dj;
 
369
 
 
370
 
 
371
      // Update the terms of L
 
372
      Index rs = size-j-1;
 
373
      w.tail(rs) -= wj * mat.col(j).tail(rs);
 
374
      if(gamma != 0)
 
375
        mat.col(j).tail(rs) += (sigma*conj(wj)/gamma)*w.tail(rs);
 
376
    }
 
377
    return true;
 
378
  }
 
379
 
 
380
  template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
 
381
  static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w, typename MatrixType::RealScalar sigma=1)
 
382
  {
 
383
    // Apply the permutation to the input w
 
384
    tmp = transpositions * w;
 
385
 
 
386
    return ldlt_inplace<Lower>::updateInPlace(mat,tmp,sigma);
 
387
  }
320
388
};
321
389
 
322
390
template<> struct ldlt_inplace<Upper>
327
395
    Transpose<MatrixType> matt(mat);
328
396
    return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
329
397
  }
 
398
 
 
399
  template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
 
400
  static EIGEN_STRONG_INLINE bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& tmp, WType& w, typename MatrixType::RealScalar sigma=1)
 
401
  {
 
402
    Transpose<MatrixType> matt(mat);
 
403
    return ldlt_inplace<Lower>::update(matt, transpositions, tmp, w.conjugate(), sigma);
 
404
  }
330
405
};
331
406
 
332
407
template<typename MatrixType> struct LDLT_Traits<MatrixType,Lower>
333
408
{
334
 
  typedef TriangularView<MatrixType, UnitLower> MatrixL;
335
 
  typedef TriangularView<typename MatrixType::AdjointReturnType, UnitUpper> MatrixU;
336
 
  inline static MatrixL getL(const MatrixType& m) { return m; }
337
 
  inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); }
 
409
  typedef const TriangularView<const MatrixType, UnitLower> MatrixL;
 
410
  typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitUpper> MatrixU;
 
411
  static inline MatrixL getL(const MatrixType& m) { return m; }
 
412
  static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
338
413
};
339
414
 
340
415
template<typename MatrixType> struct LDLT_Traits<MatrixType,Upper>
341
416
{
342
 
  typedef TriangularView<typename MatrixType::AdjointReturnType, UnitLower> MatrixL;
343
 
  typedef TriangularView<MatrixType, UnitUpper> MatrixU;
344
 
  inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); }
345
 
  inline static MatrixU getU(const MatrixType& m) { return m; }
 
417
  typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitLower> MatrixL;
 
418
  typedef const TriangularView<const MatrixType, UnitUpper> MatrixU;
 
419
  static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); }
 
420
  static inline MatrixU getU(const MatrixType& m) { return m; }
346
421
};
347
422
 
348
423
} // end namespace internal
367
442
  return *this;
368
443
}
369
444
 
 
445
/** Update the LDLT decomposition:  given A = L D L^T, efficiently compute the decomposition of A + sigma w w^T.
 
446
 * \param w a vector to be incorporated into the decomposition.
 
447
 * \param sigma a scalar, +1 for updates and -1 for "downdates," which correspond to removing previously-added column vectors. Optional; default value is +1.
 
448
 * \sa setZero()
 
449
  */
 
450
template<typename MatrixType, int _UpLo>
 
451
template<typename Derived>
 
452
LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::rankUpdate(const MatrixBase<Derived>& w,typename NumTraits<typename MatrixType::Scalar>::Real sigma)
 
453
{
 
454
  const Index size = w.rows();
 
455
  if (m_isInitialized)
 
456
  {
 
457
    eigen_assert(m_matrix.rows()==size);
 
458
  }
 
459
  else
 
460
  {    
 
461
    m_matrix.resize(size,size);
 
462
    m_matrix.setZero();
 
463
    m_transpositions.resize(size);
 
464
    for (Index i = 0; i < size; i++)
 
465
      m_transpositions.coeffRef(i) = i;
 
466
    m_temporary.resize(size);
 
467
    m_sign = sigma>=0 ? 1 : -1;
 
468
    m_isInitialized = true;
 
469
  }
 
470
 
 
471
  internal::ldlt_inplace<UpLo>::update(m_matrix, m_transpositions, m_temporary, w, sigma);
 
472
 
 
473
  return *this;
 
474
}
 
475
 
370
476
namespace internal {
371
477
template<typename _MatrixType, int _UpLo, typename Rhs>
372
478
struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
481
587
  return LDLT<PlainObject>(derived());
482
588
}
483
589
 
 
590
} // end namespace Eigen
 
591
 
484
592
#endif // EIGEN_LDLT_H