Factorisation uses only Int64

Author: filikat

highs/ipm/hipo/auxiliary/Auxiliary.cpp

@@ -4,129 +4,129 @@
 
 namespace hipo {
 
-void counts2Ptr(std::vector<Int>& ptr, std::vector<Int>& w) {
+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.
 
-  Int temp_nz{};
-  Int n = w.size();
-  for (Int j = 0; j < n; ++j) {
+  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<Int>& perm, std::vector<Int>& iperm) {
+void inversePerm(const std::vector<Int64>& perm, std::vector<Int64>& 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 (Int i = 0; i < perm.size(); ++i) {
+  for (Int64 i = 0; i < perm.size(); ++i) {
     iperm[perm[i]] = i;
   }
 }
 
-void subtreeSize(const std::vector<Int>& parent, std::vector<Int>& sizes) {
+void subtreeSize(const std::vector<Int64>& parent, std::vector<Int64>& sizes) {
   // Compute sizes of subtrees of the tree given by parent
 
-  Int n = parent.size();
+  Int64 n = parent.size();
   sizes.assign(n, 1);
 
-  for (Int i = 0; i < n; ++i) {
-    Int k = parent[i];
+  for (Int64 i = 0; i < n; ++i) {
+    Int64 k = parent[i];
     if (k != -1) sizes[k] += sizes[i];
   }
 }
 
-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<Int64>& ptr, const std::vector<Int64>& rows,
+               std::vector<Int64>& ptrT, std::vector<Int64>& rowsT) {
   // Compute the transpose of the matrix and return it in rowsT and ptrT
 
-  Int n = ptr.size() - 1;
+  Int64 n = ptr.size() - 1;
 
-  std::vector<Int> work(n);
+  std::vector<Int64> work(n);
 
   // count the entries in each row into work
-  for (Int i = 0; i < ptr.back(); ++i) {
+  for (Int64 i = 0; i < ptr.back(); ++i) {
     ++work[rows[i]];
   }
 
   // sum row sums to obtain pointers
   counts2Ptr(ptrT, work);
 
-  for (Int j = 0; j < n; ++j) {
-    for (Int el = ptr[j]; el < ptr[j + 1]; ++el) {
-      Int i = rows[el];
+  for (Int64 j = 0; j < n; ++j) {
+    for (Int64 el = ptr[j]; el < ptr[j + 1]; ++el) {
+      Int64 i = rows[el];
 
       // entry (i,j) becomes entry (j,i)
-      Int pos = work[i]++;
+      Int64 pos = work[i]++;
       rowsT[pos] = j;
     }
   }
 }
 
-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 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) {
   // Compute the transpose of the matrix and return it in rowsT, ptrT and valT
 
-  Int n = ptr.size() - 1;
+  Int64 n = ptr.size() - 1;
 
-  std::vector<Int> work(n);
+  std::vector<Int64> work(n);
 
   // count the entries in each row into work
-  for (Int i = 0; i < ptr.back(); ++i) {
+  for (Int64 i = 0; i < ptr.back(); ++i) {
     ++work[rows[i]];
   }
 
   // sum row sums to obtain pointers
   counts2Ptr(ptrT, work);
 
-  for (Int j = 0; j < n; ++j) {
-    for (Int el = ptr[j]; el < ptr[j + 1]; ++el) {
-      Int i = rows[el];
+  for (Int64 j = 0; j < n; ++j) {
+    for (Int64 el = ptr[j]; el < ptr[j + 1]; ++el) {
+      Int64 i = rows[el];
 
       // entry (i,j) becomes entry (j,i)
-      Int pos = work[i]++;
+      Int64 pos = work[i]++;
       rowsT[pos] = j;
       valT[pos] = val[el];
     }
   }
 }
 
-void childrenLinkedList(const std::vector<Int>& parent, std::vector<Int>& head,
-                        std::vector<Int>& next) {
+void childrenLinkedList(const std::vector<Int64>& parent,
+                        std::vector<Int64>& head, std::vector<Int64>& next) {
   // Create linked lists of children in elimination tree.
   // parent gives the dependencies of the tree,
   // head[node] is the first child of node,
   // next[head[node]] is the second child,
   // next[next[head[node]]] is the third child...
   // until -1 is reached.
 
-  Int n = parent.size();
+  Int64 n = parent.size();
   head.assign(n, -1);
   next.assign(n, -1);
-  for (Int node = n - 1; node >= 0; --node) {
+  for (Int64 node = n - 1; node >= 0; --node) {
     if (parent[node] == -1) continue;
     next[node] = head[parent[node]];
     head[parent[node]] = node;
   }
 }
 
-void reverseLinkedList(std::vector<Int>& head, std::vector<Int>& next) {
+void reverseLinkedList(std::vector<Int64>& head, std::vector<Int64>& next) {
   // Reverse the linked list of children of each node.
   // If a node has children (a -> b -> c -> -1), the reverse list contains
   // children (c -> b -> a -> -1).
 
-  const Int n = head.size();
+  const Int64 n = head.size();
 
-  for (Int node = 0; node < n; ++node) {
-    Int prev_node = -1;
-    Int curr_node = head[node];
-    Int next_node = -1;
+  for (Int64 node = 0; node < n; ++node) {
+    Int64 prev_node = -1;
+    Int64 curr_node = head[node];
+    Int64 next_node = -1;
 
     while (curr_node != -1) {
       next_node = next[curr_node];
@@ -139,18 +139,18 @@ void reverseLinkedList(std::vector<Int>& head, std::vector<Int>& next) {
   }
 }
 
-void dfsPostorder(Int node, Int& start, std::vector<Int>& head,
-                  const std::vector<Int>& next, std::vector<Int>& order) {
+void dfsPostorder(Int64 node, Int64& start, std::vector<Int64>& head,
+                  const std::vector<Int64>& next, std::vector<Int64>& 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.
 
-  std::stack<Int> stack;
+  std::stack<Int64> stack;
   stack.push(node);
 
   while (!stack.empty()) {
-    const Int current = stack.top();
-    const Int child = head[current];
+    const Int64 current = stack.top();
+    const Int64 child = head[current];
 
     if (child == -1) {
       // no children left to order,
@@ -166,9 +166,9 @@ void dfsPostorder(Int node, Int& start, std::vector<Int>& head,
   }
 }
 
-void processEdge(Int j, Int i, const std::vector<Int>& first,
-                 std::vector<Int>& maxfirst, std::vector<Int>& delta,
-                 std::vector<Int>& prevleaf, std::vector<Int>& ancestor) {
+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) {
   // Process edge of skeleton matrix.
   // Taken from Tim Davis "Direct Methods for Sparse Linear Systems".
 
@@ -181,21 +181,21 @@ void processEdge(Int j, Int i, const std::vector<Int>& first,
   maxfirst[i] = first[j];
 
   // previous leaf of ith row subtree
-  Int jprev = prevleaf[i];
+  Int64 jprev = prevleaf[i];
 
   // A(i,j) is in the skeleton matrix
   delta[j]++;
 
   if (jprev != -1) {
     // find least common ancestor of jprev and j
-    Int q = jprev;
+    Int64 q = jprev;
     while (q != ancestor[q]) {
       q = ancestor[q];
     }
 
     // path compression
-    Int sparent;
-    for (Int s = jprev; s != q; s = sparent) {
+    Int64 sparent;
+    for (Int64 s = jprev; s != q; s = sparent) {
       sparent = ancestor[s];
       ancestor[s] = q;
     }
@@ -208,18 +208,18 @@ void processEdge(Int j, Int i, const std::vector<Int>& first,
   prevleaf[i] = j;
 }
 
-double getDiagStart(Int n, Int k, Int nb, Int n_blocks, std::vector<Int>& start,
-                    bool triang) {
+Int64 getDiagStart(Int64 n, Int64 k, Int64 nb, Int64 n_blocks,
+                   std::vector<Int64>& start, bool triang) {
   // start position of diagonal blocks for blocked dense formats
   start.assign(n_blocks, 0);
-  for (Int i = 1; i < n_blocks; ++i) {
+  for (Int64 i = 1; i < n_blocks; ++i) {
     start[i] = start[i - 1] + nb * (n - (i - 1) * nb);
     if (triang) start[i] -= nb * (nb - 1) / 2;
   }
 
-  Int jb = std::min(nb, k - (n_blocks - 1) * nb);
-  double result = (double)start.back() + (double)(n - (n_blocks - 1) * nb) * jb;
-  if (triang) result -= (double)jb * (jb - 1) / 2;
+  Int64 jb = std::min(nb, k - (n_blocks - 1) * nb);
+  Int64 result = start.back() + (n - (n_blocks - 1) * nb) * jb;
+  if (triang) result -= jb * (jb - 1) / 2;
   return result;
 }
 

highs/ipm/hipo/auxiliary/Auxiliary.h

@@ -11,37 +11,37 @@
 
 namespace hipo {
 
-void counts2Ptr(std::vector<Int>& ptr, std::vector<Int>& w);
-void inversePerm(const std::vector<Int>& perm, std::vector<Int>& iperm);
-void subtreeSize(const std::vector<Int>& parent, std::vector<Int>& sizes);
-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<Int>& parent, std::vector<Int>& head,
-                        std::vector<Int>& next);
-void reverseLinkedList(std::vector<Int>& head, std::vector<Int>& next);
-void dfsPostorder(Int node, Int& start, std::vector<Int>& head,
-                  const std::vector<Int>& next, std::vector<Int>& order);
-void processEdge(Int j, Int i, const std::vector<Int>& first,
-                 std::vector<Int>& maxfirst, std::vector<Int>& delta,
-                 std::vector<Int>& prevleaf, std::vector<Int>& ancestor);
-double getDiagStart(Int n, Int k, Int nb, Int n_blocks, std::vector<Int>& start,
-                    bool triang = false);
+void counts2Ptr(std::vector<Int64>& ptr, std::vector<Int64>& w);
+void inversePerm(const std::vector<Int64>& perm, std::vector<Int64>& 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 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);
+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<Int>& perm) {
+void permuteVector(std::vector<T>& v, const std::vector<Int64>& perm) {
   // Permute vector v according to permutation perm.
   std::vector<T> temp_v(v);
-  for (Int i = 0; i < v.size(); ++i) v[i] = temp_v[perm[i]];
+  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<Int>& iperm) {
+void permuteVectorInverse(std::vector<T>& v, const std::vector<Int64>& iperm) {
   // Permute vector v according to inverse permutation iperm.
   std::vector<T> temp_v(v);
-  for (Int i = 0; i < v.size(); ++i) v[iperm[i]] = temp_v[i];
+  for (Int64 i = 0; i < v.size(); ++i) v[iperm[i]] = temp_v[i];
 }
 
 template <typename T>