Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-08 13:23:51 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Break LUP factorization when invalid matrix found

Browse Source 
Singular matrices cannot be LUP factorized, so this condition is
detected and returned in a `valid` flag. However, when the detection
happened, the rest of the computation continued to happen. This could
lead to floating point exceptions, which some applications do not like.
So, when an invalid array is detected, return immediately.
This commit is contained in:
Kenneth Moreland 2021-07-09 13:24:51 -06:00
parent 6447b17303
commit f74a2239eb
1 changed files with 14 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
vtkm/Matrix.h
View File

@ -273,6 +273,7 @@ VTKM_EXEC_CONT void MatrixLUPFactorFindPivot(vtkm::Matrix<T, Size, Size>& A,
if (maxValue < vtkm::Epsilon<T>())
{
valid = false;
return;
}
if (maxRowIndex != topCornerIndex)
@ -295,15 +296,14 @@ VTKM_EXEC_CONT void MatrixLUPFactorFindPivot(vtkm::Matrix<T, Size, Size>& A,
// Used with MatrixLUPFactor
template <typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT void MatrixLUPFactorFindUpperTriangleElements(vtkm::Matrix<T, Size, Size>& A,
vtkm::IdComponent topCornerIndex)
vtkm::IdComponent topCornerIndex,
bool& valid)
{
// Compute values for upper triangle on row topCornerIndex
if (A(topCornerIndex, topCornerIndex) == 0)
{
for (vtkm::IdComponent colIndex = topCornerIndex + 1; colIndex < Size; colIndex++)
{
A(topCornerIndex, colIndex) = std::numeric_limits<T>::quiet_NaN();
}
valid = false;
return;
}
else
{
@ -374,7 +374,15 @@ VTKM_EXEC_CONT void MatrixLUPFactor(vtkm::Matrix<T, Size, Size>& A,
for (vtkm::IdComponent rowIndex = 0; rowIndex < Size; rowIndex++)
{
MatrixLUPFactorFindPivot(A, permutation, rowIndex, inversionParity, valid);
MatrixLUPFactorFindUpperTriangleElements(A, rowIndex);
if (!valid)
{
break;
}
MatrixLUPFactorFindUpperTriangleElements(A, rowIndex, valid);
if (!valid)
{
break;
}
}
}