Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://www.mooseframework.org
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #pragma once
      11             : 
      12             : #include "KokkosHeader.h"
      13             : 
      14             : namespace Moose::Kokkos
      15             : {
      16             : namespace Utils
      17             : {
      18             : 
      19             : /**
      20             :  * Returns the sign of a value
      21             :  * @param x The value
      22             :  * @returns The sign of the value
      23             :  */
      24             : template <typename T>
      25             : KOKKOS_INLINE_FUNCTION T
      26    21599928 : sign(T x)
      27             : {
      28    21599928 :   return x >= 0.0 ? 1.0 : -1.0;
      29             : }
      30             : 
      31             : /**
      32             :  * Find a value in an array
      33             :  * @param target The target value to find
      34             :  * @param begin The pointer to the first element of the array
      35             :  * @param end The pointer next to the last element of the array
      36             :  * @returns The pointer to the target element, \p end if the target element was not found
      37             :  */
      38             : template <typename T>
      39             : KOKKOS_INLINE_FUNCTION const T *
      40    14226948 : find(const T & target, const T * const begin, const T * const end)
      41             : {
      42    14226948 :   if (begin == end)
      43           0 :     return end;
      44             : 
      45    14226948 :   auto left = begin;
      46    14226948 :   auto right = end - 1;
      47             : 
      48    37943375 :   while (left <= right)
      49             :   {
      50    37393285 :     auto mid = left + (right - left) / 2;
      51             : 
      52    37393285 :     if (*mid == target)
      53    13676858 :       return mid;
      54    23716427 :     else if (*mid < target)
      55    14276083 :       left = mid + 1;
      56             :     else
      57     9440344 :       right = mid - 1;
      58             :   }
      59             : 
      60      550090 :   return end;
      61             : }
      62             : 
      63             : /**
      64             :  * Perform an in-place linear solve using Cholesky decomposition
      65             :  * Matrix and right-hand-side vector are modified after this call
      66             :  * @param A The row-major matrix
      67             :  * @param x The solution vector
      68             :  * @param b The right-hand-side vector
      69             :  * @param n The system size
      70             :  */
      71             : KOKKOS_INLINE_FUNCTION void
      72      140252 : choleskySolve(Real * const A, Real * const x, Real * const b, const unsigned int n)
      73             : {
      74      979460 :   for (unsigned int i = 0; i < n; ++i)
      75             :   {
      76     3771828 :     for (unsigned int j = 0; j <= i; ++j)
      77             :     {
      78     2932620 :       Real sum = A[j + n * i];
      79             : 
      80     7812640 :       for (unsigned int k = 0; k < j; ++k)
      81     4880020 :         sum -= A[k + n * i] * A[k + n * j];
      82             : 
      83     2932620 :       if (i == j)
      84      839208 :         A[j + n * i] = ::Kokkos::sqrt(sum);
      85             :       else
      86     2093412 :         A[j + n * i] = sum / A[j + n * j];
      87             :     }
      88             :   }
      89             : 
      90      979460 :   for (unsigned int i = 0; i < n; ++i)
      91             :   {
      92      839208 :     Real sum = b[i];
      93             : 
      94     2932620 :     for (unsigned int j = 0; j < i; ++j)
      95     2093412 :       sum -= A[j + n * i] * b[j];
      96             : 
      97      839208 :     b[i] = sum / A[i + n * i];
      98             :   }
      99             : 
     100      979460 :   for (int i = n - 1; i >= 0; --i)
     101             :   {
     102      839208 :     Real sum = b[i];
     103             : 
     104     2932620 :     for (unsigned int j = i + 1; j < n; ++j)
     105     2093412 :       sum -= A[i + n * j] * x[j];
     106             : 
     107      839208 :     x[i] = sum / A[i + n * i];
     108             :   }
     109      140252 : }
     110             : 
     111             : } // namespace Utils
     112             : } // namespace Moose::Kokkos