|
MOAB
4.9.3pre
|
|
MOAB
4.9.3pre
|
#include <PartialPivLU.h>
Public Types | |
| typedef Map< Matrix< Scalar, Dynamic, Dynamic, StorageOrder > > | MapLU |
| typedef Block< MapLU, Dynamic, Dynamic > | MatrixType |
| typedef Block< MatrixType, Dynamic, Dynamic > | BlockType |
| typedef MatrixType::RealScalar | RealScalar |
Static Public Member Functions | |
| static Index | unblocked_lu (MatrixType &lu, PivIndex *row_transpositions, PivIndex &nb_transpositions) |
| static Index | blocked_lu (Index rows, Index cols, Scalar *lu_data, Index luStride, PivIndex *row_transpositions, PivIndex &nb_transpositions, Index maxBlockSize=256) |
This is the blocked version of fullpivlu_unblocked()
Definition at line 289 of file PartialPivLU.h.
| typedef Block<MatrixType,Dynamic,Dynamic> Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex >::BlockType |
Definition at line 298 of file PartialPivLU.h.
| typedef Map<Matrix<Scalar, Dynamic, Dynamic, StorageOrder> > Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex >::MapLU |
Definition at line 296 of file PartialPivLU.h.
| typedef Block<MapLU, Dynamic, Dynamic> Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex >::MatrixType |
Definition at line 297 of file PartialPivLU.h.
| typedef MatrixType::RealScalar Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex >::RealScalar |
Definition at line 299 of file PartialPivLU.h.
| static Index Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex >::blocked_lu | ( | Index | rows, |
| Index | cols, | ||
| Scalar * | lu_data, | ||
| Index | luStride, | ||
| PivIndex * | row_transpositions, | ||
| PivIndex & | nb_transpositions, | ||
| Index | maxBlockSize = 256 |
||
| ) | [inline, static] |
performs the LU decomposition in-place of the matrix represented by the variables rows, cols, lu_data, and lu_stride using a recursive, blocked algorithm.
In addition, this function returns the row transpositions in the vector row_transpositions which must have a size equal to the number of columns of the matrix lu, and an integer nb_transpositions which returns the actual number of transpositions.
Definition at line 372 of file PartialPivLU.h.
{
MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols);
MatrixType lu(lu1,0,0,rows,cols);
const Index size = (std::min)(rows,cols);
// if the matrix is too small, no blocking:
if(size<=16)
{
return unblocked_lu(lu, row_transpositions, nb_transpositions);
}
// automatically adjust the number of subdivisions to the size
// of the matrix so that there is enough sub blocks:
Index blockSize;
{
blockSize = size/8;
blockSize = (blockSize/16)*16;
blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize);
}
nb_transpositions = 0;
Index first_zero_pivot = -1;
for(Index k = 0; k < size; k+=blockSize)
{
Index bs = (std::min)(size-k,blockSize); // actual size of the block
Index trows = rows - k - bs; // trailing rows
Index tsize = size - k - bs; // trailing size
// partition the matrix:
// A00 | A01 | A02
// lu = A_0 | A_1 | A_2 = A10 | A11 | A12
// A20 | A21 | A22
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,tsize);
BlockType A21(lu,k+bs,k,trows,bs);
BlockType A22(lu,k+bs,k+bs,trows,tsize);
PivIndex nb_transpositions_in_panel;
// recursively call the blocked LU algorithm on [A11^T A21^T]^T
// with a very small blocking size:
Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
row_transpositions+k, nb_transpositions_in_panel, 16);
if(ret>=0 && first_zero_pivot==-1)
first_zero_pivot = k+ret;
nb_transpositions += nb_transpositions_in_panel;
// update permutations and apply them to A_0
for(Index i=k; i<k+bs; ++i)
{
Index piv = (row_transpositions[i] += k);
A_0.row(i).swap(A_0.row(piv));
}
if(trows)
{
// apply permutations to A_2
for(Index i=k;i<k+bs; ++i)
A_2.row(i).swap(A_2.row(row_transpositions[i]));
// A12 = A11^-1 A12
A11.template triangularView<UnitLower>().solveInPlace(A12);
A22.noalias() -= A21 * A12;
}
}
return first_zero_pivot;
}
| static Index Eigen::internal::partial_lu_impl< Scalar, StorageOrder, PivIndex >::unblocked_lu | ( | MatrixType & | lu, |
| PivIndex * | row_transpositions, | ||
| PivIndex & | nb_transpositions | ||
| ) | [inline, static] |
performs the LU decomposition in-place of the matrix lu using an unblocked algorithm.
In addition, this function returns the row transpositions in the vector row_transpositions which must have a size equal to the number of columns of the matrix lu, and an integer nb_transpositions which returns the actual number of transpositions.
Definition at line 311 of file PartialPivLU.h.
{
typedef scalar_score_coeff_op<Scalar> Scoring;
typedef typename Scoring::result_type Score;
const Index rows = lu.rows();
const Index cols = lu.cols();
const Index size = (std::min)(rows,cols);
nb_transpositions = 0;
Index first_zero_pivot = -1;
for(Index k = 0; k < size; ++k)
{
Index rrows = rows-k-1;
Index rcols = cols-k-1;
Index row_of_biggest_in_col;
Score biggest_in_corner
= lu.col(k).tail(rows-k).unaryExpr(Scoring()).maxCoeff(&row_of_biggest_in_col);
row_of_biggest_in_col += k;
row_transpositions[k] = PivIndex(row_of_biggest_in_col);
if(biggest_in_corner != Score(0))
{
if(k != row_of_biggest_in_col)
{
lu.row(k).swap(lu.row(row_of_biggest_in_col));
++nb_transpositions;
}
// FIXME shall we introduce a safe quotient expression in cas 1/lu.coeff(k,k)
// overflow but not the actual quotient?
lu.col(k).tail(rrows) /= lu.coeff(k,k);
}
else if(first_zero_pivot==-1)
{
// the pivot is exactly zero, we record the index of the first pivot which is exactly 0,
// and continue the factorization such we still have A = PLU
first_zero_pivot = k;
}
if(k<rows-1)
lu.bottomRightCorner(rrows,rcols).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rcols);
}
return first_zero_pivot;
}
1.7.6.1