defaultmatrixhelper.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_GEOMETRY_UTILITY_DEFAULTMATRIXHELPER_HH
#define DUNE_GEOMETRY_UTILITY_DEFAULTMATRIXHELPER_HH
 
#include <cmath>
 
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
 
namespace Dune {
namespace Impl {
 
// FieldMatrixHelper
// -----------------
 
template< class ct >
struct FieldMatrixHelper
{
  using ctype = ct;
 
  template< int m, int n >
  [[ deprecated("Use A.mv(x,y) instead.") ]]
  static void Ax ( const FieldMatrix< ctype, m, n > &A, const FieldVector< ctype, n > &x, FieldVector< ctype, m > &ret )
  {
    A.mv(x,ret);
  }
 
  template< int m, int n >
  [[ deprecated("Use A.mtv(x,y) instead.") ]]
  static void ATx ( const FieldMatrix< ctype, m, n > &A, const FieldVector< ctype, m > &x, FieldVector< ctype, n > &ret )
  {
    A.mtv(x,ret);
  }
 
  template< int m, int n, int p >
  [[ deprecated("Use FMatrixHelp::multMatrix(A,B,ret)") ]]
  static void AB ( const FieldMatrix< ctype, m, n > &A, const FieldMatrix< ctype, n, p > &B, FieldMatrix< ctype, m, p > &ret )
  {
    FMatrixHelp::multMatrix(A,B,ret);
  }
 
  template< int m, int n, int p >
  static void ATBT ( const FieldMatrix< ctype, m, n > &A, const FieldMatrix< ctype, p, m > &B, FieldMatrix< ctype, n, p > &ret )
  {
    for( int i = 0; i < n; ++i )
    {
      for( int j = 0; j < p; ++j )
      {
        ret[ i ][ j ] = ctype( 0 );
        for( int k = 0; k < m; ++k )
          ret[ i ][ j ] += A[ k ][ i ] * B[ j ][ k ];
      }
    }
  }
 
  template< int m, int n >
  static void ATA_L ( const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, n > &ret )
  {
    for( int i = 0; i < n; ++i )
    {
      for( int j = 0; j <= i; ++j )
      {
        ret[ i ][ j ] = ctype( 0 );
        for( int k = 0; k < m; ++k )
          ret[ i ][ j ] += A[ k ][ i ] * A[ k ][ j ];
      }
    }
  }
 
  template< int m, int n >
  [[ deprecated("Use FMatrixHelp::multTransposedMatrix(A,ret)") ]]
  static void ATA ( const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, n > &ret )
  {
    return FMatrixHelp::multTransposedMatrix(A,ret);
  }
 
  template< int m, int n >
  static void AAT_L ( const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, m, m > &ret )
  {
    for( int i = 0; i < m; ++i )
    {
      for( int j = 0; j <= i; ++j )
      {
        ctype &retij = ret[ i ][ j ];
        retij = A[ i ][ 0 ] * A[ j ][ 0 ];
        for( int k = 1; k < n; ++k )
          retij += A[ i ][ k ] * A[ j ][ k ];
      }
    }
  }
 
  template< int m, int n >
  static void AAT ( const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, m, m > &ret )
  {
    for( int i = 0; i < m; ++i )
    {
      for( int j = 0; j < i; ++j )
      {
        ret[ i ][ j ] = ctype( 0 );
        for( int k = 0; k < n; ++k )
          ret[ i ][ j ] += A[ i ][ k ] * A[ j ][ k ];
        ret[ j ][ i ] = ret[ i ][ j ];
      }
      ret[ i ][ i ] = ctype( 0 );
      for( int k = 0; k < n; ++k )
        ret[ i ][ i ] += A[ i ][ k ] * A[ i ][ k ];
    }
  }
 
  // [[ expects: L is lower triangular ]]
  template< int n >
  static void Lx ( const FieldMatrix< ctype, n, n > &L, const FieldVector< ctype, n > &x, FieldVector< ctype, n > &ret )
  {
    for( int i = 0; i < n; ++i )
    {
      ret[ i ] = ctype( 0 );
      for( int j = 0; j <= i; ++j )
        ret[ i ] += L[ i ][ j ] * x[ j ];
    }
  }
 
  // [[ expects: L is lower triangular ]]
  template< int n >
  static void LTx ( const FieldMatrix< ctype, n, n > &L, const FieldVector< ctype, n > &x, FieldVector< ctype, n > &ret )
  {
    for( int i = 0; i < n; ++i )
    {
      ret[ i ] = ctype( 0 );
      for( int j = i; j < n; ++j )
        ret[ i ] += L[ j ][ i ] * x[ j ];
    }
  }
 
  // [[ expects: L is lower triangular ]]
  template< int n >
  static void LTL ( const FieldMatrix< ctype, n, n > &L, FieldMatrix< ctype, n, n > &ret )
  {
    for( int i = 0; i < n; ++i )
    {
      for( int j = 0; j < i; ++j )
      {
        ret[ i ][ j ] = ctype( 0 );
        for( int k = i; k < n; ++k )
          ret[ i ][ j ] += L[ k ][ i ] * L[ k ][ j ];
        ret[ j ][ i ] = ret[ i ][ j ];
      }
      ret[ i ][ i ] = ctype( 0 );
      for( int k = i; k < n; ++k )
        ret[ i ][ i ] += L[ k ][ i ] * L[ k ][ i ];
    }
  }
 
  // [[ expects: L is lower triangular ]]
  template< int n >
  static void LLT ( const FieldMatrix< ctype, n, n > &L, FieldMatrix< ctype, n, n > &ret )
  {
    for( int i = 0; i < n; ++i )
    {
      for( int j = 0; j < i; ++j )
      {
        ret[ i ][ j ] = ctype( 0 );
        for( int k = 0; k <= j; ++k )
          ret[ i ][ j ] += L[ i ][ k ] * L[ j ][ k ];
        ret[ j ][ i ] = ret[ i ][ j ];
      }
      ret[ i ][ i ] = ctype( 0 );
      for( int k = 0; k <= i; ++k )
        ret[ i ][ i ] += L[ i ][ k ] * L[ i ][ k ];
    }
  }
 
  // [[ expects: A is spd ]]
  template< int n >
  static bool cholesky_L ( const FieldMatrix< ctype, n, n > &A, FieldMatrix< ctype, n, n > &ret, const bool checkSingular = false )
  {
    using std::sqrt;
    for( int i = 0; i < n; ++i )
    {
      ctype &rii = ret[ i ][ i ];
 
      ctype xDiag = A[ i ][ i ];
      for( int j = 0; j < i; ++j )
        xDiag -= ret[ i ][ j ] * ret[ i ][ j ];
 
      // in some cases A can be singular, e.g. when checking local for
      // outside points during checkInside
      if( checkSingular && ! ( xDiag > ctype( 0 )) )
        return false ;
 
      // otherwise this should be true always
      assert( xDiag > ctype( 0 ) );
      rii = sqrt( xDiag );
 
      ctype invrii = ctype( 1 ) / rii;
      for( int k = i+1; k < n; ++k )
      {
        ctype x = A[ k ][ i ];
        for( int j = 0; j < i; ++j )
          x -= ret[ i ][ j ] * ret[ k ][ j ];
        ret[ k ][ i ] = invrii * x;
      }
    }
 
    // return true for meaning A is non-singular
    return true;
  }
 
  // [[ expects: L is lower triangular ]]
  template< int n >
  static ctype detL ( const FieldMatrix< ctype, n, n > &L )
  {
    ctype det( 1 );
    for( int i = 0; i < n; ++i )
      det *= L[ i ][ i ];
    return det;
  }
 
  // [[ expects: L is lower triangular ]]
  template< int n >
  static ctype invL ( FieldMatrix< ctype, n, n > &L )
  {
    ctype det( 1 );
    for( int i = 0; i < n; ++i )
    {
      ctype &lii = L[ i ][ i ];
      det *= lii;
      lii = ctype( 1 ) / lii;
      for( int j = 0; j < i; ++j )
      {
        ctype &lij = L[ i ][ j ];
        ctype x = lij * L[ j ][ j ];
        for( int k = j+1; k < i; ++k )
          x += L[ i ][ k ] * L[ k ][ j ];
        lij = (-lii) * x;
      }
    }
    return det;
  }
 
  // [[ expects: L is lower triangular ]]
  template< int n >
  static void invLx ( FieldMatrix< ctype, n, n > &L, FieldVector< ctype, n > &x )
  {
    for( int i = 0; i < n; ++i )
    {
      for( int j = 0; j < i; ++j )
        x[ i ] -= L[ i ][ j ] * x[ j ];
      x[ i ] /= L[ i ][ i ];
    }
  }
 
  // [[ expects: L is lower triangular ]]
  template< int n >
  static void invLTx ( FieldMatrix< ctype, n, n > &L, FieldVector< ctype, n > &x )
  {
    for( int i = n; i > 0; --i )
    {
      for( int j = i; j < n; ++j )
        x[ i-1 ] -= L[ j ][ i-1 ] * x[ j ];
      x[ i-1 ] /= L[ i-1 ][ i-1 ];
    }
  }
 
  // [[ expects: A is spd ]]
  template< int n >
  static ctype spdDetA ( const FieldMatrix< ctype, n, n > &A )
  {
    FieldMatrix< ctype, n, n > L;
    cholesky_L( A, L );
    return detL( L );
  }
 
  // [[ expects: A is spd ]]
  template< int n >
  static ctype spdInvA ( FieldMatrix< ctype, n, n > &A )
  {
    FieldMatrix< ctype, n, n > L;
    cholesky_L( A, L );
    const ctype det = invL( L );
    LTL( L, A );
    return det;
  }
 
  // [[ expects: A is spd ]]
  template< int n >
  static bool spdInvAx ( FieldMatrix< ctype, n, n > &A, FieldVector< ctype, n > &x, const bool checkSingular = false )
  {
    FieldMatrix< ctype, n, n > L;
    const bool invertible = cholesky_L( A, L, checkSingular );
    if( ! invertible ) return invertible ;
    invLx( L, x );
    invLTx( L, x );
    return invertible;
  }
 
  template< int m, int n >
  static ctype detATA ( const FieldMatrix< ctype, m, n > &A )
  {
    if constexpr( m >= n )
    {
      FieldMatrix< ctype, n, n > ata;
      ATA_L( A, ata );
      return spdDetA( ata );
    }
    else
      return ctype( 0 );
  }
 
  template< int m, int n >
  static ctype sqrtDetAAT ( const FieldMatrix< ctype, m, n > &A )
  {
    using std::abs;
    using std::sqrt;
    // These special cases are here not only for speed reasons:
    // The general implementation aborts if the matrix is almost singular,
    // and the special implementation provide a stable way to handle that case.
    if constexpr( (n == 2) && (m == 2) )
    {
      // Special implementation for 2x2 matrices: faster and more stable
      return abs( A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ] );
    }
    else if constexpr( (n == 3) && (m == 3) )
    {
      // Special implementation for 3x3 matrices
      const ctype v0 = A[ 0 ][ 1 ] * A[ 1 ][ 2 ] - A[ 1 ][ 1 ] * A[ 0 ][ 2 ];
      const ctype v1 = A[ 0 ][ 2 ] * A[ 1 ][ 0 ] - A[ 1 ][ 2 ] * A[ 0 ][ 0 ];
      const ctype v2 = A[ 0 ][ 0 ] * A[ 1 ][ 1 ] - A[ 1 ][ 0 ] * A[ 0 ][ 1 ];
      return abs( v0 * A[ 2 ][ 0 ] + v1 * A[ 2 ][ 1 ] + v2 * A[ 2 ][ 2 ] );
    }
    else if constexpr( (n == 3) && (m == 2) )
    {
      // Special implementation for 2x3 matrices
      const ctype v0 = A[ 0 ][ 0 ] * A[ 1 ][ 1 ] - A[ 0 ][ 1 ] * A[ 1 ][ 0 ];
      const ctype v1 = A[ 0 ][ 0 ] * A[ 1 ][ 2 ] - A[ 1 ][ 0 ] * A[ 0 ][ 2 ];
      const ctype v2 = A[ 0 ][ 1 ] * A[ 1 ][ 2 ] - A[ 0 ][ 2 ] * A[ 1 ][ 1 ];
      return sqrt( v0*v0 + v1*v1 + v2*v2);
    }
    else if constexpr( n >= m )
    {
      // General case
      FieldMatrix< ctype, m, m > aat;
      AAT_L( A, aat );
      return spdDetA( aat );
    }
    else
      return ctype( 0 );
  }
 
  // => A^{-1}_L A = I
  template< int m, int n >
  static ctype leftInvA ( const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, m > &ret )
  {
    using std::abs;
    if constexpr( (n == 2) && (m == 2) )
    {
      const ctype det = (A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ]);
      const ctype detInv = ctype( 1 ) / det;
      ret[ 0 ][ 0 ] = A[ 1 ][ 1 ] * detInv;
      ret[ 1 ][ 1 ] = A[ 0 ][ 0 ] * detInv;
      ret[ 1 ][ 0 ] = -A[ 1 ][ 0 ] * detInv;
      ret[ 0 ][ 1 ] = -A[ 0 ][ 1 ] * detInv;
      return abs( det );
    }
    else
    {
      FieldMatrix< ctype, n, n > ata;
      ATA_L( A, ata );
      const ctype det = spdInvA( ata );
      ATBT( ata, A, ret );
      return det;
    }
  }
 
  template< int m, int n >
  static bool leftInvAx ( const FieldMatrix< ctype, m, n > &A, const FieldVector< ctype, m > &x, FieldVector< ctype, n > &y )
  {
    static_assert((m >= n), "Matrix has no left inverse.");
    FieldMatrix< ctype, n, n > ata;
    A.mtv(x, y);
    ATA_L( A, ata );
    return spdInvAx( ata, y, true );
  }
 
  template< int m, int n >
  static ctype rightInvA ( const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, m > &ret )
  {
    static_assert((n >= m), "Matrix has no right inverse.");
    using std::abs;
    if constexpr( (n == 2) && (m == 2) )
    {
      const ctype det = (A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ]);
      const ctype detInv = ctype( 1 ) / det;
      ret[ 0 ][ 0 ] = A[ 1 ][ 1 ] * detInv;
      ret[ 1 ][ 1 ] = A[ 0 ][ 0 ] * detInv;
      ret[ 1 ][ 0 ] = -A[ 1 ][ 0 ] * detInv;
      ret[ 0 ][ 1 ] = -A[ 0 ][ 1 ] * detInv;
      return abs( det );
    }
    else
    {
      FieldMatrix< ctype, m , m > aat;
      AAT_L( A, aat );
      const ctype det = spdInvA( aat );
      ATBT( A , aat , ret );
      return det;
    }
  }
 
  template< int m, int n >
  static bool xTRightInvA ( const FieldMatrix< ctype, m, n > &A, const FieldVector< ctype, n > &x, FieldVector< ctype, m > &y )
  {
    static_assert((n >= m), "Matrix has no right inverse.");
    FieldMatrix< ctype, m, m > aat;
    A.mv(x, y);
    AAT_L( A, aat );
    // check whether aat is singular and return true if non-singular
    return spdInvAx( aat, y, true );
  }
};
 
} // namespace Impl
} // namespace Dune
 
#endif // DUNE_GEOMETRY_UTILITY_DEFAULTMATRIXHELPER_HH