265 |
{ |
{ |
266 |
TinyMatrix<N,N,T> M; |
TinyMatrix<N,N,T> M; |
267 |
TinyMatrix<N,N,T> lu; |
TinyMatrix<N,N,T> lu; |
268 |
LU(*this, lu); |
TinyMatrix<N,N,T> P; |
269 |
|
LU(*this, lu, P); |
270 |
|
|
271 |
TinyVector<N,T> v; |
TinyVector<N,T> v; |
272 |
TinyVector<N,T> u; |
TinyVector<N,T> u; |
275 |
for (size_t j=0; j<N; ++j) |
for (size_t j=0; j<N; ++j) |
276 |
v[j] = (i==j)?1:0; |
v[j] = (i==j)?1:0; |
277 |
|
|
278 |
u = lu.applyLU(v); |
u = lu.applyLU(P*v); |
279 |
for (size_t j=0; j<N; ++j) |
for (size_t j=0; j<N; ++j) |
280 |
M(j,i) = u[j]; |
M(j,i) = u[j]; |
281 |
} |
} |
521 |
} |
} |
522 |
} |
} |
523 |
|
|
524 |
/*! Returns the LU factorisation of a matrix |
/*! Computes P, L and U such that PA=LU. |
525 |
|
L and U are stored in a single matrix. L non zero coefficients are i<j. |
526 |
\remark As only sens for square matrix |
\remark As only sens for square matrix |
527 |
*/ |
*/ |
528 |
friend void LU(const TinyMatrix<N,N,T>& A, |
friend void LU(const TinyMatrix<N,N,T>& A, |
529 |
TinyMatrix<N,N,T>& lu) |
TinyMatrix<N,N,T>& lu, |
530 |
|
TinyMatrix<N,N,T>& P) |
531 |
{ |
{ |
532 |
|
// Reinitialize P to Identity |
533 |
|
for (size_t i=0; i<N; i++) |
534 |
|
for (size_t j=0; j<N; j++) { |
535 |
|
P(i,j) = (i==j); |
536 |
|
} |
537 |
|
|
538 |
lu = A; |
lu = A; |
539 |
|
|
540 |
for (int k=0; k<N-1; k++) { |
for (int k=0; k<N-1; k++) { |
541 |
#warning Do not proceed to swaping anymore when needed. |
// Checkes for Pivot. |
542 |
// int maxLine = k; |
int maxLine = k; |
543 |
// for (int i=k+1; i<N; i++) { |
for (int i=k+1; i<N; i++) { |
544 |
// maxLine = (std::abs(lu(i,k))<std::abs(lu(maxLine,k)))?maxLine:i; |
maxLine = (std::abs(lu(i,k))<std::abs(lu(maxLine,k)))?maxLine:i; |
545 |
// } |
} |
546 |
// lu.swapLines(k,maxLine); |
lu.swapLines(k,maxLine); |
547 |
|
P.swapLines(k,maxLine); |
548 |
|
|
549 |
for (int i=k+1; i<N; i++) { |
for (int i=k+1; i<N; i++) { |
550 |
lu(i,k) /= lu(k,k); |
lu(i,k) /= lu(k,k); |
551 |
for (int j=k+1; j<N; j++) |
for (int j=k+1; j<N; j++) |
586 |
const TinyMatrix<N,N,T>& A) |
const TinyMatrix<N,N,T>& A) |
587 |
{ |
{ |
588 |
TinyMatrix<N,N,T> lu; |
TinyMatrix<N,N,T> lu; |
589 |
LU(A,lu); |
TinyMatrix<N,N,T> P; |
590 |
|
LU(A,lu,P); |
591 |
|
|
592 |
return lu.applyLU(b); |
return lu.applyLU(P*b); |
593 |
} |
} |
594 |
|
|
595 |
//! Constructor does nothing. |
//! Constructor does nothing. |