diff options
author | 2009-08-06 16:41:54 +0200 | |
---|---|---|
committer | 2009-08-06 16:41:54 +0200 | |
commit | 3ac01b1400dccefcb2db85ab1caf88fc3bd379db (patch) | |
tree | 74ff5edd66389c3c24966e9d635ffe96a0ab0107 /Eigen/src/LU | |
parent | e82e30862aa4af91ee2d33e554406f20ee7b4155 (diff) |
compilation fix in EigenSolver,
bugfix in PartialLU
Diffstat (limited to 'Eigen/src/LU')
-rw-r--r-- | Eigen/src/LU/PartialLU.h | 31 |
1 files changed, 15 insertions, 16 deletions
diff --git a/Eigen/src/LU/PartialLU.h b/Eigen/src/LU/PartialLU.h index 1cca1c7db..10d39055a 100644 --- a/Eigen/src/LU/PartialLU.h +++ b/Eigen/src/LU/PartialLU.h @@ -245,11 +245,10 @@ struct ei_partial_lu_impl if(k<rows-1) { - lu.col(k).end(rows-k-1) /= lu.coeff(k,k); - - // TODO implement a fast rank one update routine - for(int col = k + 1; col < size; ++col) - lu.col(col).end(rows-k-1) -= lu.col(k).end(rows-k-1) * lu.coeff(k,col); + int rrows = rows-k-1; + int rsize = size-k-1; + lu.col(k).end(rrows) /= lu.coeff(k,k); + lu.corner(BottomRight,rrows,rsize) -= (lu.col(k).end(rrows) * lu.row(k).end(rsize)).lazy(); } } } @@ -271,7 +270,6 @@ struct ei_partial_lu_impl { MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols); MatrixType lu(lu1,0,0,rows,cols); - const int size = std::min(rows,cols); @@ -294,24 +292,25 @@ struct ei_partial_lu_impl nb_transpositions = 0; for(int k = 0; k < size; k+=blockSize) { - int bs = std::min(size-k,blockSize); - int ps = size - k; - int rs = size - k - bs; + int bs = std::min(size-k,blockSize); // actual size of the block + int trows = rows - k - bs; // trailing rows + int tsize = size - k - bs; // trailing size + // partition the matrix: // A00 | A01 | A02 // lu = A10 | A11 | A12 // A20 | A21 | A22 - BlockType A_0(lu,0,0,size,k); - BlockType A_2(lu,0,k+bs,size,rs); + BlockType A_0(lu,0,0,rows,k); + BlockType A_2(lu,0,k+bs,rows,tsize); BlockType A11(lu,k,k,bs,bs); - BlockType A12(lu,k,k+bs,bs,rs); - BlockType A21(lu,k+bs,k,rs,bs); - BlockType A22(lu,k+bs,k+bs,rs,rs); + BlockType A12(lu,k,k+bs,bs,tsize); + BlockType A21(lu,k+bs,k,trows,bs); + BlockType A22(lu,k+bs,k+bs,trows,tsize); int nb_transpositions_in_panel; // recursively calls the blocked LU algorithm with a very small // blocking size: - blocked_lu(ps, bs, &lu.coeffRef(k,k), luStride, + blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride, row_transpositions+k, nb_transpositions_in_panel, 16); nb_transpositions_in_panel += nb_transpositions_in_panel; @@ -322,7 +321,7 @@ struct ei_partial_lu_impl A_0.row(i).swap(A_0.row(piv)); } - if(rs) + if(trows) { // apply permutations to A_2 for(int i=k;i<k+bs; ++i) |