Matrix is formed with Int, factorisation uses Int64

Author: filikat

highs/io/HighsIO.cpp

@@ -29,14 +29,20 @@ void highsLogHeader(const HighsLogOptions& log_options,
                githash_text.c_str(), kHighsCopyrightStatement.c_str());
 
 #ifdef HIPO
+  std::string blas_model =
+      hipo::getIntegerModelString(hipo::getBlasIntegerModel());
+
 #ifdef BLAS_LIBRARIES
-  highsLogUser(
-      log_options, HighsLogType::kInfo, "Using blas: %s - %s\n", BLAS_LIBRARIES,
-      hipo::getIntegerModelString(hipo::getBlasIntegerModel()).c_str());
+  highsLogUser(log_options, HighsLogType::kInfo, "Using blas: %s - %s\n",
+               BLAS_LIBRARIES, blas_model.c_str());
+#else
+#ifdef HIPO_USES_OPENBLAS
+  highsLogUser(log_options, HighsLogType::kInfo, "Using blas: OpenBLAS - %s\n",
+               blas_model.c_str());
 #else
-  highsLogUser(
-      log_options, HighsLogType::kInfo, "Using blas: unknown - %s\n",
-      hipo::getIntegerModelString(hipo::getBlasIntegerModel()).c_str());
+  highsLogUser(log_options, HighsLogType::kInfo, "Using blas: unknown - %s\n",
+               blas_model.c_str());
+#endif
 #endif
 #endif
 }

highs/ipm/IpxWrapper.cpp

@@ -1521,9 +1521,9 @@ HighsStatus reportHipoStatus(const HighsOptions& options,
   else if (status == hipo::kStatusError) {
     highsLogUser(options.log_options, HighsLogType::kError,
                  "Hipo: Internal error\n");
-  } else if (status == hipo::kStatusOoM) {
+  } else if (status == hipo::kStatusOverflow) {
     highsLogUser(options.log_options, HighsLogType::kError,
-                 "Hipo: Out of memory\n");
+                 "Hipo: Integer overflow\n");
   } else if (status == hipo::kStatusErrorAnalyse) {
     highsLogUser(options.log_options, HighsLogType::kError,
                  "Hipo: Error in analyse phase\n");

highs/ipm/hipo/auxiliary/AutoDetect.cpp

@@ -2,14 +2,15 @@
 
 #include <stdint.h>
 
-#define IDXTYPEWIDTH 64
+// declare Metis with 32-bit integers for the test
+#define IDXTYPEWIDTH 32
 #include "metis.h"
 
 // Weird tricks to detect the integer type width used by BLAS and Metis.
 // These are technically undefined behaviour, because they rely on using a
-// function declaration that involves int64_t, while the actual implementation
-// may use a different type. Their behaviour may depend on the endianness of the
-// CPU (?).
+// function declaration that involves a certain integer type, while the actual
+// implementation may use a different one. Their behaviour may depend on the
+// endianness of the CPU (?).
 
 namespace hipo {
 
@@ -80,10 +81,9 @@ IntegerModel getMetisIntegerModel() {
     idx_t options[METIS_NOPTIONS];
     for (int i = 0; i < METIS_NOPTIONS; ++i) options[i] = -1;
 
-    // if Metis is using 32-bit internally, this should set ptype to 0 and
-    // objtype to 1, which should trigger an error. If it uses 64-bit then
-    // everything should be fine.
-    options[METIS_OPTION_PTYPE] = 1;
+    // if 32 bits are used, this sets iptype to 2, otherwise it sets objtype to
+    // a wrong value
+    options[METIS_OPTION_IPTYPE] = 2;
 
     idx_t n = 3;
     idx_t ptr[4] = {0, 2, 4, 6};
@@ -96,9 +96,9 @@ IntegerModel getMetisIntegerModel() {
       if (perm[0] != 0 || perm[1] != 1 || perm[2] != 2)
         metis_model = IntegerModel::unknown;
       else
-        metis_model = IntegerModel::ilp64;
+        metis_model = IntegerModel::lp64;
     } else if (status == METIS_ERROR_INPUT) {
-      metis_model = IntegerModel::lp64;
+      metis_model = IntegerModel::ilp64;
     } else {
       metis_model = IntegerModel::unknown;
     }

highs/ipm/hipo/auxiliary/Auxiliary.cpp

@@ -4,27 +4,12 @@
 
 namespace hipo {
 
-void counts2Ptr(std::vector<Int64>& ptr, std::vector<Int64>& w) {
-  // Given the column counts in the vector w (of size n),
-  // compute the column pointers in the vector ptr (of size n+1),
-  // and copy the first n pointers back into w.
-
-  Int64 temp_nz{};
-  Int64 n = w.size();
-  for (Int64 j = 0; j < n; ++j) {
-    ptr[j] = temp_nz;
-    temp_nz += w[j];
-    w[j] = ptr[j];
-  }
-  ptr[n] = temp_nz;
-}
-
-void inversePerm(const std::vector<Int64>& perm, std::vector<Int64>& iperm) {
+void inversePerm(const std::vector<Int>& perm, std::vector<Int>& iperm) {
   // Given the permutation perm, produce the inverse permutation iperm.
   // perm[i] : i-th entry to use in the new order.
   // iperm[i]: where entry i is located in the new order.
 
-  for (Int64 i = 0; i < perm.size(); ++i) {
+  for (Int i = 0; i < perm.size(); ++i) {
     iperm[perm[i]] = i;
   }
 }
@@ -41,56 +26,56 @@ void subtreeSize(const std::vector<Int64>& parent, std::vector<Int64>& sizes) {
   }
 }
 
-void transpose(const std::vector<Int64>& ptr, const std::vector<Int64>& rows,
-               std::vector<Int64>& ptrT, std::vector<Int64>& rowsT) {
+void transpose(const std::vector<Int>& ptr, const std::vector<Int>& rows,
+               std::vector<Int>& ptrT, std::vector<Int>& rowsT) {
   // Compute the transpose of the matrix and return it in rowsT and ptrT
 
-  Int64 n = ptr.size() - 1;
+  Int n = ptr.size() - 1;
 
-  std::vector<Int64> work(n);
+  std::vector<Int> work(n);
 
   // count the entries in each row into work
-  for (Int64 i = 0; i < ptr.back(); ++i) {
+  for (Int i = 0; i < ptr.back(); ++i) {
     ++work[rows[i]];
   }
 
   // sum row sums to obtain pointers
   counts2Ptr(ptrT, work);
 
-  for (Int64 j = 0; j < n; ++j) {
-    for (Int64 el = ptr[j]; el < ptr[j + 1]; ++el) {
-      Int64 i = rows[el];
+  for (Int j = 0; j < n; ++j) {
+    for (Int el = ptr[j]; el < ptr[j + 1]; ++el) {
+      Int i = rows[el];
 
       // entry (i,j) becomes entry (j,i)
-      Int64 pos = work[i]++;
+      Int pos = work[i]++;
       rowsT[pos] = j;
     }
   }
 }
 
-void transpose(const std::vector<Int64>& ptr, const std::vector<Int64>& rows,
-               const std::vector<double>& val, std::vector<Int64>& ptrT,
-               std::vector<Int64>& rowsT, std::vector<double>& valT) {
+void transpose(const std::vector<Int>& ptr, const std::vector<Int>& rows,
+               const std::vector<double>& val, std::vector<Int>& ptrT,
+               std::vector<Int>& rowsT, std::vector<double>& valT) {
   // Compute the transpose of the matrix and return it in rowsT, ptrT and valT
 
-  Int64 n = ptr.size() - 1;
+  Int n = ptr.size() - 1;
 
-  std::vector<Int64> work(n);
+  std::vector<Int> work(n);
 
   // count the entries in each row into work
-  for (Int64 i = 0; i < ptr.back(); ++i) {
+  for (Int i = 0; i < ptr.back(); ++i) {
     ++work[rows[i]];
   }
 
   // sum row sums to obtain pointers
   counts2Ptr(ptrT, work);
 
-  for (Int64 j = 0; j < n; ++j) {
-    for (Int64 el = ptr[j]; el < ptr[j + 1]; ++el) {
-      Int64 i = rows[el];
+  for (Int j = 0; j < n; ++j) {
+    for (Int el = ptr[j]; el < ptr[j + 1]; ++el) {
+      Int i = rows[el];
 
       // entry (i,j) becomes entry (j,i)
-      Int64 pos = work[i]++;
+      Int pos = work[i]++;
       rowsT[pos] = j;
       valT[pos] = val[el];
     }
@@ -140,7 +125,7 @@ void reverseLinkedList(std::vector<Int64>& head, std::vector<Int64>& next) {
 }
 
 void dfsPostorder(Int64 node, Int64& start, std::vector<Int64>& head,
-                  const std::vector<Int64>& next, std::vector<Int64>& order) {
+                  const std::vector<Int64>& next, std::vector<Int>& order) {
   // Perform depth first search starting from root node and order the nodes
   // starting from the value start. head and next contain the linked list of
   // children.

highs/ipm/hipo/auxiliary/Auxiliary.h

@@ -11,34 +11,49 @@
 
 namespace hipo {
 
-void counts2Ptr(std::vector<Int64>& ptr, std::vector<Int64>& w);
-void inversePerm(const std::vector<Int64>& perm, std::vector<Int64>& iperm);
+void inversePerm(const std::vector<Int>& perm, std::vector<Int>& iperm);
 void subtreeSize(const std::vector<Int64>& parent, std::vector<Int64>& sizes);
-void transpose(const std::vector<Int64>& ptr, const std::vector<Int64>& rows,
-               std::vector<Int64>& ptrT, std::vector<Int64>& rowsT);
-void transpose(const std::vector<Int64>& ptr, const std::vector<Int64>& rows,
-               const std::vector<double>& val, std::vector<Int64>& ptrT,
-               std::vector<Int64>& rowsT, std::vector<double>& valT);
+void transpose(const std::vector<Int>& ptr, const std::vector<Int>& rows,
+               std::vector<Int>& ptrT, std::vector<Int>& rowsT);
+void transpose(const std::vector<Int>& ptr, const std::vector<Int>& rows,
+               const std::vector<double>& val, std::vector<Int>& ptrT,
+               std::vector<Int>& rowsT, std::vector<double>& valT);
 void childrenLinkedList(const std::vector<Int64>& parent,
                         std::vector<Int64>& head, std::vector<Int64>& next);
 void reverseLinkedList(std::vector<Int64>& head, std::vector<Int64>& next);
 void dfsPostorder(Int64 node, Int64& start, std::vector<Int64>& head,
-                  const std::vector<Int64>& next, std::vector<Int64>& order);
+                  const std::vector<Int64>& next, std::vector<Int>& order);
 void processEdge(Int64 j, Int64 i, const std::vector<Int64>& first,
                  std::vector<Int64>& maxfirst, std::vector<Int64>& delta,
                  std::vector<Int64>& prevleaf, std::vector<Int64>& ancestor);
 Int64 getDiagStart(Int64 n, Int64 k, Int64 nb, Int64 n_blocks,
                    std::vector<Int64>& start, bool triang = false);
 
 template <typename T>
-void permuteVector(std::vector<T>& v, const std::vector<Int64>& perm) {
+void counts2Ptr(std::vector<T>& ptr, std::vector<T>& w) {
+  // Given the column counts in the vector w (of size n),
+  // compute the column pointers in the vector ptr (of size n+1),
+  // and copy the first n pointers back into w.
+
+  T temp_nz{};
+  T n = w.size();
+  for (T j = 0; j < n; ++j) {
+    ptr[j] = temp_nz;
+    temp_nz += w[j];
+    w[j] = ptr[j];
+  }
+  ptr[n] = temp_nz;
+}
+
+template <typename T>
+void permuteVector(std::vector<T>& v, const std::vector<Int>& perm) {
   // Permute vector v according to permutation perm.
   std::vector<T> temp_v(v);
   for (Int64 i = 0; i < v.size(); ++i) v[i] = temp_v[perm[i]];
 }
 
 template <typename T>
-void permuteVectorInverse(std::vector<T>& v, const std::vector<Int64>& iperm) {
+void permuteVectorInverse(std::vector<T>& v, const std::vector<Int>& iperm) {
   // Permute vector v according to inverse permutation iperm.
   std::vector<T> temp_v(v);
   for (Int64 i = 0; i < v.size(); ++i) v[iperm[i]] = temp_v[i];

highs/ipm/hipo/auxiliary/IntConfig.h

@@ -9,17 +9,17 @@ namespace hipo {
 // Generic integer type from HiGHS
 typedef HighsInt Int;
 
-// Integer type for indices of matrices in HiPO and factorisation
+// Integer type for factorisation
 typedef int64_t Int64;
 
-// The matrix (AS or NE) size must fit into Int.
-// Type Int64 is used only for the nonzeros of the matrix and during the
-// factorisation.
+// The matrix (AS or NE) is formed using Int, so it must have fewer than
+// kHighsIInf nonzero entries. Metis works with the same type as Int, so it must
+// be compiled accordingly.
 //
-// For NE, AS, factorisations:
-// - ptr and rows are std::vector<Index>.
-// - ptr can be accessed with Int, rows and val must be accessed with Index.
-// - rows[i] can be stored as Int, ptr[i] must be stored as Index.
+// The factorisation uses Int64 everywhere, apart from where it interfaces with
+// the matrix stored using Int.
+// BLAS is 32-bit, so the vectors used by BLAS must be addressable with 32-bit
+// integers.
 //
 
 }  // namespace hipo