This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
View | Details | Raw Unified | Return to bug 608
Collapse All | Expand All

(-)a/Eigen/src/Cholesky/LDLT.h (-16 / +17 lines)
Lines 254-270 template<> struct ldlt_inplace<Lower> Link Here
254
    typedef typename MatrixType::Index Index;
254
    typedef typename MatrixType::Index Index;
255
    eigen_assert(mat.rows()==mat.cols());
255
    eigen_assert(mat.rows()==mat.cols());
256
    const Index size = mat.rows();
256
    const Index size = mat.rows();
257
257
258
    if (size <= 1)
258
    if (size <= 1)
259
    {
259
    {
260
      transpositions.setIdentity();
260
      transpositions.setIdentity();
261
      if(sign)
261
      if(sign)
262
        *sign = real(mat.coeff(0,0))>0 ? 1:-1;
262
        *sign = math::real(mat.coeff(0,0))>0 ? 1:-1;
263
      return true;
263
      return true;
264
    }
264
    }
265
265
266
    RealScalar cutoff(0), biggest_in_corner;
266
    RealScalar cutoff(0), biggest_in_corner;
267
267
268
    for (Index k = 0; k < size; ++k)
268
    for (Index k = 0; k < size; ++k)
269
    {
269
    {
270
      // Find largest diagonal element
270
      // Find largest diagonal element
Lines 273-324 template<> struct ldlt_inplace<Lower> Link Here
273
      index_of_biggest_in_corner += k;
273
      index_of_biggest_in_corner += k;
274
274
275
      if(k == 0)
275
      if(k == 0)
276
      {
276
      {
277
        // The biggest overall is the point of reference to which further diagonals
277
        // The biggest overall is the point of reference to which further diagonals
278
        // are compared; if any diagonal is negligible compared
278
        // are compared; if any diagonal is negligible compared
279
        // to the largest overall, the algorithm bails.
279
        // to the largest overall, the algorithm bails.
280
        cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);
280
        cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);
281
282
        if(sign)
283
          *sign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1;
284
      }
285
      else if(sign)
286
      {
287
        // LDLT is not guaranteed to work for indefinite matrices, but let's try to get the sign right
288
        int newSign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0;
289
        if(newSign != *sign)
290
          *sign = 0;
291
      }
281
      }
292
282
293
      // Finish early if the matrix is not full rank.
283
      // Finish early if the matrix is not full rank.
294
      if(biggest_in_corner < cutoff)
284
      if(biggest_in_corner < cutoff)
295
      {
285
      {
296
        for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
286
        for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
287
        if(sign) *sign = 0;
297
        break;
288
        break;
298
      }
289
      }
299
290
300
      transpositions.coeffRef(k) = index_of_biggest_in_corner;
291
      transpositions.coeffRef(k) = index_of_biggest_in_corner;
301
      if(k != index_of_biggest_in_corner)
292
      if(k != index_of_biggest_in_corner)
302
      {
293
      {
303
        // apply the transposition while taking care to consider only
294
        // apply the transposition while taking care to consider only
304
        // the lower triangular part
295
        // the lower triangular part
305
        Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element
296
        Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element
306
        mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
297
        mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
307
        mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
298
        mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
308
        std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
299
        std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
309
        for(int i=k+1;i<index_of_biggest_in_corner;++i)
300
        for(int i=k+1;i<index_of_biggest_in_corner;++i)
310
        {
301
        {
311
          Scalar tmp = mat.coeffRef(i,k);
302
          Scalar tmp = mat.coeffRef(i,k);
312
          mat.coeffRef(i,k) = conj(mat.coeffRef(index_of_biggest_in_corner,i));
303
          mat.coeffRef(i,k) = math::conj(mat.coeffRef(index_of_biggest_in_corner,i));
313
          mat.coeffRef(index_of_biggest_in_corner,i) = conj(tmp);
304
          mat.coeffRef(index_of_biggest_in_corner,i) = math::conj(tmp);
314
        }
305
        }
315
        if(NumTraits<Scalar>::IsComplex)
306
        if(NumTraits<Scalar>::IsComplex)
316
          mat.coeffRef(index_of_biggest_in_corner,k) = conj(mat.coeff(index_of_biggest_in_corner,k));
307
          mat.coeffRef(index_of_biggest_in_corner,k) = math::conj(mat.coeff(index_of_biggest_in_corner,k));
317
      }
308
      }
318
309
319
      // partition the matrix:
310
      // partition the matrix:
320
      //       A00 |  -  |  -
311
      //       A00 |  -  |  -
321
      // lu  = A10 | A11 |  -
312
      // lu  = A10 | A11 |  -
322
      //       A20 | A21 | A22
313
      //       A20 | A21 | A22
323
      Index rs = size - k - 1;
314
      Index rs = size - k - 1;
324
      Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
315
      Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
Lines 329-344 template<> struct ldlt_inplace<Lower> Link Here
329
      {
320
      {
330
        temp.head(k) = mat.diagonal().head(k).asDiagonal() * A10.adjoint();
321
        temp.head(k) = mat.diagonal().head(k).asDiagonal() * A10.adjoint();
331
        mat.coeffRef(k,k) -= (A10 * temp.head(k)).value();
322
        mat.coeffRef(k,k) -= (A10 * temp.head(k)).value();
332
        if(rs>0)
323
        if(rs>0)
333
          A21.noalias() -= A20 * temp.head(k);
324
          A21.noalias() -= A20 * temp.head(k);
334
      }
325
      }
335
      if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
326
      if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
336
        A21 /= mat.coeffRef(k,k);
327
        A21 /= mat.coeffRef(k,k);
328
      
329
      if(sign)
330
      {
331
        // LDLT is not guaranteed to work for indefinite matrices, but let's try to get the sign right
332
        int newSign = math::real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0;
333
        if(k == 0)
334
          *sign = newSign;
335
        else if(*sign != newSign)
336
          *sign = 0;
337
      }
337
    }
338
    }
338
339
339
    return true;
340
    return true;
340
  }
341
  }
341
342
342
  // Reference for the algorithm: Davis and Hager, "Multiple Rank
343
  // Reference for the algorithm: Davis and Hager, "Multiple Rank
343
  // Modifications of a Sparse Cholesky Factorization" (Algorithm 1)
344
  // Modifications of a Sparse Cholesky Factorization" (Algorithm 1)
344
  // Trivial rearrangements of their computations (Timothy E. Holy)
345
  // Trivial rearrangements of their computations (Timothy E. Holy)
Lines 362-391 template<> struct ldlt_inplace<Lower> Link Here
362
    // Apply the update
363
    // Apply the update
363
    for (Index j = 0; j < size; j++)
364
    for (Index j = 0; j < size; j++)
364
    {
365
    {
365
      // Check for termination due to an original decomposition of low-rank
366
      // Check for termination due to an original decomposition of low-rank
366
      if (!(isfinite)(alpha))
367
      if (!(isfinite)(alpha))
367
        break;
368
        break;
368
369
369
      // Update the diagonal terms
370
      // Update the diagonal terms
370
      RealScalar dj = real(mat.coeff(j,j));
371
      RealScalar dj = math::real(mat.coeff(j,j));
371
      Scalar wj = w.coeff(j);
372
      Scalar wj = w.coeff(j);
372
      RealScalar swj2 = sigma*abs2(wj);
373
      RealScalar swj2 = sigma*abs2(wj);
373
      RealScalar gamma = dj*alpha + swj2;
374
      RealScalar gamma = dj*alpha + swj2;
374
375
375
      mat.coeffRef(j,j) += swj2/alpha;
376
      mat.coeffRef(j,j) += swj2/alpha;
376
      alpha += swj2/dj;
377
      alpha += swj2/dj;
377
378
378
379
379
      // Update the terms of L
380
      // Update the terms of L
380
      Index rs = size-j-1;
381
      Index rs = size-j-1;
381
      w.tail(rs) -= wj * mat.col(j).tail(rs);
382
      w.tail(rs) -= wj * mat.col(j).tail(rs);
382
      if(gamma != 0)
383
      if(gamma != 0)
383
        mat.col(j).tail(rs) += (sigma*conj(wj)/gamma)*w.tail(rs);
384
        mat.col(j).tail(rs) += (sigma*math::conj(wj)/gamma)*w.tail(rs);
384
    }
385
    }
385
    return true;
386
    return true;
386
  }
387
  }
387
388
388
  template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
389
  template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
389
  static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w, const typename MatrixType::RealScalar& sigma=1)
390
  static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w, const typename MatrixType::RealScalar& sigma=1)
390
  {
391
  {
391
    // Apply the permutation to the input w
392
    // Apply the permutation to the input w
(-)a/Eigen/src/Cholesky/LLT.h (-3 / +3 lines)
Lines 227-243 static typename MatrixType::Index llt_ra Link Here
227
    }
227
    }
228
  }
228
  }
229
  else
229
  else
230
  {
230
  {
231
    temp = vec;
231
    temp = vec;
232
    RealScalar beta = 1;
232
    RealScalar beta = 1;
233
    for(Index j=0; j<n; ++j)
233
    for(Index j=0; j<n; ++j)
234
    {
234
    {
235
      RealScalar Ljj = real(mat.coeff(j,j));
235
      RealScalar Ljj = math::real(mat.coeff(j,j));
236
      RealScalar dj = abs2(Ljj);
236
      RealScalar dj = abs2(Ljj);
237
      Scalar wj = temp.coeff(j);
237
      Scalar wj = temp.coeff(j);
238
      RealScalar swj2 = sigma*abs2(wj);
238
      RealScalar swj2 = sigma*abs2(wj);
239
      RealScalar gamma = dj*beta + swj2;
239
      RealScalar gamma = dj*beta + swj2;
240
240
241
      RealScalar x = dj + swj2/beta;
241
      RealScalar x = dj + swj2/beta;
242
      if (x<=RealScalar(0))
242
      if (x<=RealScalar(0))
243
        return j;
243
        return j;
Lines 246-262 static typename MatrixType::Index llt_ra Link Here
246
      beta += swj2/dj;
246
      beta += swj2/dj;
247
247
248
      // Update the terms of L
248
      // Update the terms of L
249
      Index rs = n-j-1;
249
      Index rs = n-j-1;
250
      if(rs)
250
      if(rs)
251
      {
251
      {
252
        temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
252
        temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
253
        if(gamma != 0)
253
        if(gamma != 0)
254
          mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*conj(wj)/gamma)*temp.tail(rs);
254
          mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*math::conj(wj)/gamma)*temp.tail(rs);
255
      }
255
      }
256
    }
256
    }
257
  }
257
  }
258
  return -1;
258
  return -1;
259
}
259
}
260
260
261
template<typename Scalar> struct llt_inplace<Scalar, Lower>
261
template<typename Scalar> struct llt_inplace<Scalar, Lower>
262
{
262
{
Lines 272-288 template<typename Scalar> struct llt_inp Link Here
272
    for(Index k = 0; k < size; ++k)
272
    for(Index k = 0; k < size; ++k)
273
    {
273
    {
274
      Index rs = size-k-1; // remaining size
274
      Index rs = size-k-1; // remaining size
275
275
276
      Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
276
      Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
277
      Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
277
      Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
278
      Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
278
      Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
279
279
280
      RealScalar x = real(mat.coeff(k,k));
280
      RealScalar x = math::real(mat.coeff(k,k));
281
      if (k>0) x -= A10.squaredNorm();
281
      if (k>0) x -= A10.squaredNorm();
282
      if (x<=RealScalar(0))
282
      if (x<=RealScalar(0))
283
        return k;
283
        return k;
284
      mat.coeffRef(k,k) = x = sqrt(x);
284
      mat.coeffRef(k,k) = x = sqrt(x);
285
      if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
285
      if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
286
      if (rs>0) A21 *= RealScalar(1)/x;
286
      if (rs>0) A21 *= RealScalar(1)/x;
287
    }
287
    }
288
    return -1;
288
    return -1;
(-)a/test/cholesky.cpp (-4 / +12 lines)
Lines 263-282 template<typename MatrixType> void chole Link Here
263
263
264
// LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal.
264
// LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal.
265
// This test checks that LDLT reports correctly that matrix is indefinite. 
265
// This test checks that LDLT reports correctly that matrix is indefinite. 
266
// See http://forum.kde.org/viewtopic.php?f=74&t=106942
266
// See http://forum.kde.org/viewtopic.php?f=74&t=106942
267
template<typename MatrixType> void cholesky_indefinite(const MatrixType& m)
267
template<typename MatrixType> void cholesky_indefinite(const MatrixType& m)
268
{
268
{
269
  eigen_assert(m.rows() == 2 && m.cols() == 2);
269
  eigen_assert(m.rows() == 2 && m.cols() == 2);
270
  MatrixType mat;
270
  MatrixType mat;
271
  mat << 1, 0, 0, -1;
271
  {
272
  LDLT<MatrixType> ldlt(mat);
272
    mat << 1, 0, 0, -1;
273
  VERIFY(!ldlt.isNegative());
273
    LDLT<MatrixType> ldlt(mat);
274
  VERIFY(!ldlt.isPositive());
274
    VERIFY(!ldlt.isNegative());
275
    VERIFY(!ldlt.isPositive());
276
  }
277
  {
278
    mat << 1, 2, 2, 1;
279
    LDLT<MatrixType> ldlt(mat);
280
    VERIFY(!ldlt.isNegative());
281
    VERIFY(!ldlt.isPositive());
282
  }
275
}
283
}
276
284
277
template<typename MatrixType> void cholesky_verify_assert()
285
template<typename MatrixType> void cholesky_verify_assert()
278
{
286
{
279
  MatrixType tmp;
287
  MatrixType tmp;
280
288
281
  LLT<MatrixType> llt;
289
  LLT<MatrixType> llt;
282
  VERIFY_RAISES_ASSERT(llt.matrixL())
290
  VERIFY_RAISES_ASSERT(llt.matrixL())

Return to bug 608