Determinant of 3x3 with std::fma

The determinant for a 3x3 matrix in blaze/math/expressions/DMatDetExpr.h can be subject to fatal floating point cancellation. An example is given here: https://pharr.org/matt/blog/2019/11/03/difference-of-floats

This could be solved by using the proposed function DifferenceOfProducts to calculate the determinant of a 2x2 matrix and build the determinant of a 3x3 matrix with

/*! \cond BLAZE_INTERNAL */
/*!\brief Computation of the determinant of the given dense \f$ 3 \times 3 \f$ matrix.
// \ingroup dense_matrix
//
// \param dm The given dense matrix.
// \return The determinant of the given matrix.
//
// This function computes the determinant of the given dense \f$ 3 \times 3 \f$ matrix via the
// rule of Sarrus.
*/
template< typename MT  // Type of the dense matrix
        , bool SO >    // Storage order of the dense matrix
inline ElementType_t<MT> det3x3( const DenseMatrix<MT,SO>& dm )
{
   BLAZE_INTERNAL_ASSERT( (*dm).rows()    == 3UL, "Invalid number of rows detected"    );
   BLAZE_INTERNAL_ASSERT( (*dm).columns() == 3UL, "Invalid number of columns detected" );

   CompositeType_t<MT> A( *dm );

  // Here det2 is the determinant of a 2x2 matrix
    return   det2
             (
                A(0,0),det2(A(1,1),A(2,2),A(1,2),A(2,1)),
                A(0,1),det2(A(1,0),A(2,2),A(1,2),A(2,0))
             )
           + A(0,2)*det2(A(1,0),A(2,1),A(1,1),A(2,0));
}


//*************************************************************************************************
/*! \cond BLAZE_INTERNAL */
/*!\brief Computation of the determinant of the given dense \f$ 2 \times 2 \f$ matrix.
// \ingroup dense_matrix
//
// \param dm The given dense matrix.
// \return The determinant of the given matrix.
//
// This function computes the determinant of the given dense \f$ 2 \times 2 \f$ matrix via the
// rule of Sarrus.
*/
template< typename MT  // Type of the dense matrix
        , bool SO >    // Storage order of the dense matrix
inline ElementType_t<MT> det2x2( const DenseMatrix<MT,SO>& dm )
{
   BLAZE_INTERNAL_ASSERT( (*dm).rows()    == 2UL, "Invalid number of rows detected"    );
   BLAZE_INTERNAL_ASSERT( (*dm).columns() == 2UL, "Invalid number of columns detected" );

    CompositeType_t<MT> A( *dm );

    // For a matrix of the form A = a  b
    //                              c  d
    double cb = A(1,0) * A(0,1);
    float err = std::fma(-A(1,0), A(0,1), cb);
    float dop = std::fma(A(0,0), A(1,1), -cb);
    return dop + err;
}

‌

Comments (1)