Implement "declamping" for SuiteSparse Matrices, fix formatting for

C-style casts

Author: BenBrock

.clang-format

@@ -4,4 +4,5 @@ PointerAlignment: Left
 ColumnLimit: 80
 AlwaysBreakTemplateDeclarations: Yes
 AllowShortFunctionsOnASingleLine: Empty
+SpaceAfterCStyleCast: true
 ---

examples/check_equivalence.c

@@ -98,16 +98,16 @@ bsp_fdataset_info_t bsp_parse_fdataset_string(char* str) {
 
   if (split == -1) {
     bsp_fdataset_info_t info;
-    info.fname = (char*)malloc(sizeof(char) * (len + 1));
+    info.fname = (char*) malloc(sizeof(char) * (len + 1));
     strcpy(info.fname, str);
     info.dataset = NULL;
     return info;
   } else {
     bsp_fdataset_info_t info;
-    info.fname = (char*)malloc(sizeof(char) * (split + 1));
+    info.fname = (char*) malloc(sizeof(char) * (split + 1));
     strncpy(info.fname, str, split);
     info.fname[split] = '\0';
-    info.dataset = (char*)malloc(sizeof(char) * (len - split));
+    info.dataset = (char*) malloc(sizeof(char) * (len - split));
     strcpy(info.dataset, &str[split + 1]);
     return info;
   }

examples/mtx2bsp.c

@@ -9,6 +9,11 @@ int main(int argc, char** argv) {
     return 1;
   }
 
+  // NOTE: this binary automatically performs "declamping" of floating point
+  //       values greater/less than 1e308/-1e308 in order to restore inf values
+  //       in SuiteSparse Matrix Collection matrices.
+  bool perform_suitesparse_declamping = true;
+
   char* input_fname = argv[1];
   char* output_fname = argv[2];
 
@@ -41,6 +46,10 @@ int main(int argc, char** argv) {
   bsp_matrix_t matrix = bsp_mmread(input_fname);
   printf(" === Done reading. ===\n");
 
+  if (perform_suitesparse_declamping) {
+    bsp_matrix_declamp_values(matrix);
+  }
+
   matrix = bsp_matrix_minimize_values(matrix);
 
   bsp_print_matrix_info(matrix);

include/binsparse/array.h

@@ -108,7 +108,7 @@ void bsp_destroy_array_t(bsp_array_t array) {
       data[index] = value;                                                     \
     } else if (array.type == BSP_BINT8) {                                      \
       int8_t* data = array.data;                                               \
-      data[index] = ((size_t)value) % 2;                                       \
+      data[index] = ((size_t) value) % 2;                                      \
     } else if (array.type == BSP_COMPLEX_FLOAT32) {                            \
       float _Complex* data = array.data;                                       \
       data[index] = value;                                                     \

include/binsparse/binsparse.h

@@ -3,6 +3,7 @@
 #define BINSPARSE_VERSION "0.1"
 
 #include <binsparse/array.h>
+#include <binsparse/detail/detail.h>
 #include <binsparse/generate.h>
 #include <binsparse/hdf5_wrapper.h>
 #include <binsparse/matrix.h>

include/binsparse/detail/declamp_values.h

@@ -0,0 +1,25 @@
+#pragma once
+
+#include <binsparse/matrix.h>
+#include <math.h>
+
+// "Declamp" a matrix from the SuiteSparse Matrix Collection.
+// SuiteSparse Matrix Collection clamps values at 1e308 and -1e308
+// for printing.  These are almost all infinities or negative infinites.
+// Here, we "declamp" these values to restore them to infinity.
+// This allows the `bsp_matrix_minimize_values` to properly minimize
+// matrices that have infinity values.
+void bsp_matrix_declamp_values(bsp_matrix_t matrix) {
+  const double HUGE_DOUBLE = 1e308;
+  if (matrix.values.type == BSP_FLOAT64) {
+    double* values = (double*) matrix.values.data;
+
+    for (size_t i = 0; i < matrix.values.size; i++) {
+      if (values[i] >= HUGE_DOUBLE) {
+        values[i] = INFINITY;
+      } else if (values[i] <= -HUGE_DOUBLE) {
+        values[i] = -INFINITY;
+      }
+    }
+  }
+}

include/binsparse/detail/detail.h

@@ -0,0 +1,3 @@
+#pragma once
+
+#include <binsparse/detail/declamp_values.h>

include/binsparse/matrix_market/coo_sort_tools.h

@@ -21,8 +21,8 @@ static int bsp_compare_int_impl_(size_t x, size_t y) {
 
 static int bsp_coo_comparison_row_sort_operator_impl_(const void* x,
                                                       const void* y) {
-  size_t x_index = *((const size_t*)x);
-  size_t y_index = *((const size_t*)y);
+  size_t x_index = *((const size_t*) x);
+  size_t y_index = *((const size_t*) y);
 
   size_t x_i, x_j;
   size_t y_i, y_j;

include/binsparse/matrix_market/matrix_market_read.h

@@ -289,11 +289,11 @@ bsp_matrix_t bsp_mmread(char* file_path) {
   size_t max_dim =
       (metadata.nrows > metadata.ncols) ? metadata.nrows : metadata.ncols;
 
-  if (max_dim < (size_t)UINT8_MAX) {
+  if (max_dim < (size_t) UINT8_MAX) {
     index_type = BSP_UINT8;
-  } else if (max_dim < (size_t)UINT16_MAX) {
+  } else if (max_dim < (size_t) UINT16_MAX) {
     index_type = BSP_UINT16;
-  } else if (max_dim < (size_t)UINT32_MAX) {
+  } else if (max_dim < (size_t) UINT32_MAX) {
     index_type = BSP_UINT32;
   } else {
     index_type = BSP_UINT64;

include/binsparse/matrix_market/matrix_market_write.h

@@ -69,7 +69,7 @@ void bsp_mmwrite(char* file_path, bsp_matrix_t matrix) {
         bsp_array_read(matrix.indices_0, count, i);
         bsp_array_read(matrix.indices_1, count, j);
         bsp_array_read(matrix.values, count, value);
-        fprintf(f, "%zu %zu %lld\n", i + 1, j + 1, (long long)value);
+        fprintf(f, "%zu %zu %lld\n", i + 1, j + 1, (long long) value);
       } else if (mm_type == BSP_MM_REAL) {
         size_t i, j;
         double value;