The _pivots array should use PetscBLASInt when PETSc is available.

include/numerics/dense_matrix.h

@@ -24,6 +24,12 @@
 #include "libmesh/libmesh_common.h"
 #include "libmesh/dense_matrix_base.h"
 
+// For the definition of PetscBLASInt.
+#if (LIBMESH_HAVE_PETSC)
+# include "libmesh/petsc_macro.h"
+# include <petscblaslapack.h>
+#endif
+
 // C++ includes
 #include <vector>
 #include <algorithm>
@@ -646,11 +652,16 @@ class DenseMatrix : public DenseMatrixBase<T>
                    DenseMatrix<T> * VR = libmesh_nullptr);
 
   /**
-   * This array is used to store pivot indices.  May be used
-   * by whatever factorization is currently active, clients of
-   * the class should not rely on it for any reason.
+   * Array used to store pivot indices.  May be used by whatever
+   * factorization is currently active, clients of the class should
+   * not rely on it for any reason.
    */
-  std::vector<int> _pivots;
+#if (LIBMESH_HAVE_PETSC && LIBMESH_USE_REAL_NUMBERS)
+  typedef PetscBLASInt pivot_index_t;
+#else
+  typedef int pivot_index_t;
+#endif
+  std::vector<pivot_index_t> _pivots;
 
   /**
    * Companion function to _lu_decompose_lapack().  Do not use
@@ -733,8 +744,7 @@ DenseMatrix<T>::DenseMatrix(const unsigned int new_m,
   use_blas_lapack(false),
 #endif
   _val(),
-  _decomposition_type(NONE),
-  _pivots()
+  _decomposition_type(NONE)
 {
   this->resize(new_m,new_n);
 }

src/numerics/dense_matrix.C

@@ -665,7 +665,7 @@ void DenseMatrix<T>::_lu_back_substitute (const DenseVector<T> & b,
   for (unsigned int i=0; i<n_cols; ++i)
     {
       // Swap
-      if (_pivots[i] != static_cast<int>(i))
+      if (_pivots[i] != static_cast<pivot_index_t>(i))
         std::swap( z(i), z(_pivots[i]) );
 
       x(i) = z(i);
@@ -731,7 +731,7 @@ void DenseMatrix<T>::_lu_decompose ()
       // If the max was found in a different row, interchange rows.
       // Here we interchange the *entire* row, in Gaussian elimination
       // you would only interchange the subrows A(i,j) and A(p(i),j), for j>i
-      if (_pivots[i] != static_cast<int>(i))
+      if (_pivots[i] != static_cast<pivot_index_t>(i))
         {
           for (unsigned int j=0; j<n_rows; ++j)
             std::swap( A(i,j), A(_pivots[i], j) );
@@ -881,12 +881,12 @@ T DenseMatrix<T>::det ()
   for (unsigned int i=0; i<this->m(); i++)
     {
       if (this->_decomposition_type==LU)
-        if (_pivots[i] != static_cast<int>(i))
+        if (_pivots[i] != static_cast<pivot_index_t>(i))
           n_interchanges++;
 
       // Lapack pivots are 1-based!
       if (this->_decomposition_type==LU_BLAS_LAPACK)
-        if (_pivots[i] != static_cast<int>(i+1))
+        if (_pivots[i] != static_cast<pivot_index_t>(i+1))
           n_interchanges++;
 
       determinant *= (*this)(i,i);