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 |