Add TypeTensor check for Hermitian positive-definiteness

In the real-valued case, the check reduces gracefully to symmetric
positive-definiteness. We allow the user to specify a tolerance so
that matrices which fail to satisfy the checks by a "small" amount
(due to e.g. floating point error) can still be considered Hermitian
positive definite.

Author: jwpeterson

include/numerics/type_tensor.h

@@ -436,6 +436,13 @@ class TypeTensor
   void write_unformatted (std::ostream & out_stream,
                           const bool newline = true) const;
 
+  /**
+   * Returns true if the TypeTensor is Hermitian (or symmetric, when
+   * T==Real) and positive-definite to within the provided relative
+   * tolerance, false otherwise.
+   */
+  bool is_hpd(Real rel_tol = TOLERANCE*TOLERANCE) const;
+
 protected:
 
   /**

src/numerics/type_tensor.C

@@ -63,6 +63,59 @@ void TypeTensor<T>::write_unformatted (std::ostream & out_stream,
     out_stream << '\n';
 }
 
+template <typename T>
+bool TypeTensor<T>::is_hpd(Real rel_tol) const
+{
+  // Convenient reference to this object, to be used for matrix
+  // operations.
+  const auto & A = *this;
+
+  // The norm of this matrix, needed for relative tolerance checks.
+  auto A_norm = A.norm();
+
+  // The zero matrix is only positive semi-definite
+  if (A_norm == 0)
+    return false;
+
+  // An absolute tolerance (based on the user's provided relative
+  // tolerance) to be used in the floating point comparisons below.
+  auto abs_tol = rel_tol * A_norm;
+
+  // Form the complex conjugate transpose of A
+  TypeTensor<T> A_conjugate_transpose;
+  for (unsigned int i=0; i<LIBMESH_DIM; i++)
+    for (unsigned int j=0; j<LIBMESH_DIM; j++)
+      A_conjugate_transpose(i,j) = libmesh_conj(A(j,i));
+
+  // Check if Hermitian
+  if ((A - A_conjugate_transpose).norm() > abs_tol)
+    return false;
+
+  // If we made it here, then we are Hermitian, so now we just need to
+  // check if we are positive-definite by checking that all principal
+  // minors are positive. Note: the determinant of a Hermitian matrix
+  // and all principal minors are real-valued. Since we already
+  // checked that the matrix is Hermitian above, we don't bother to
+  // check the complex parts are also zero now.
+  if (std::real(A.det()) < -abs_tol)
+    return false;
+
+  // For 3x3, also check the upper 2x2 determinant
+#if LIBMESH_DIM > 2
+  if (std::real(A(0,0)*A(1,1) - A(0,1)*A(1,0)) < -abs_tol)
+    return false;
+#endif
+
+  // For 3x3 and 2x2, also check the 1x1 determinant
+#if LIBMESH_DIM > 1
+  if (std::real(A(0,0)) < -abs_tol)
+    return false;
+#endif
+
+  // If we made it here, then the matrix is Hermitian
+  // positive-definite.
+  return true;
+}
 
 
 template <>