transforming states with transposed matrix is 2,580% faster

Author: ZocoLini

benches/bench_transform_states.c

@@ -8,14 +8,16 @@
 #define CUQDYN_CONF "../tests/data/nfkb_cuqdyn_config.xml"
 
 void bench_transform_states();
+void bench_transform_states_2();
 
 int main()
 {
     CuqDynContext context = init_cuqdyn_context_from_file(CUQDYN_CONF);
     CuqdynConf *conf = get_cuqdyn_conf(context);
 
     BenchFunc funcs[] = {
-        bench_transform_states
+        bench_transform_states,
+        bench_transform_states_2
     };
 
     BenchOptions *options = create_bench_options("Transform States", funcs, sizeof(funcs) / sizeof(funcs[0]));
@@ -25,49 +27,43 @@ int main()
     return 0;
 }
 
+#define STATES_ROWS 1000
+#define STATES_COLS 16
+
 void bench_transform_states()
 {
-    #define STATES_ROWS 28
-    #define STATES_COLS 16
-
-    sunrealtype states_values[STATES_COLS][STATES_ROWS] = {
-        0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
-        1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        4, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        5, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        6, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        7, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        8, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        9, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        10, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        11, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
-        1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        4, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        5, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        6, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        7, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        8, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        9, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        10, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-        11, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
-    };
+    sunrealtype states_values[STATES_COLS][STATES_ROWS] = { 0 };
 
     SUNMatrix states = NewDenseMatrix(STATES_ROWS, STATES_COLS);
 
-    for (int i = 0; i < STATES_ROWS; ++i) 
+    for (int i = 0; i < STATES_ROWS; ++i)
     {
-        for (int j = 0; j < STATES_COLS; ++j) 
+        for (int j = 0; j < STATES_COLS; ++j)
         {
-            SM_ELEMENT_D(states, i, j) = states_values[j][i];      
+            SM_ELEMENT_D(states, i, j) = states_values[j][i];
         }
     }
 
     SUNMatrix transformed_states = transform_states(states);
 
     SUNMatDestroy(transformed_states);
 }
+
+void bench_transform_states_2()
+{
+    sunrealtype states_values[STATES_ROWS][STATES_COLS] = { 0 };
+
+    SUNMatrix states = NewDenseMatrix(STATES_COLS, STATES_ROWS);
+
+    for (int i = 0; i < STATES_ROWS; ++i)
+    {
+        for (int j = 0; j < STATES_COLS; ++j)
+        {
+            SM_ELEMENT_D(states, j, i) = states_values[i][j];
+        }
+    }
+
+    SUNMatrix transformed_states = transform_states_2(states);
+
+    SUNMatDestroy(transformed_states);
+}

benches/bench_wrapper.cpp

@@ -45,9 +45,9 @@ extern "C" void run_benchmark(BenchOptions *options)
         PAPI_cleanup_eventset(EventSet);
         PAPI_destroy_eventset(&EventSet);
 
-        std::cout << "#### Cache analysis for " << options->name << " ###\n"
-                  << "L1 cache misses: " << values[0] << ", L2 cache misses: " << values[1] << "\n"
-                  << "L1 cache hits: " << values[2] << ", L2 cache hits: " << values[3] << "\n";
+        // std::cout << "#### Cache analysis for " << options->name << " ###\n"
+        //          << "L1 cache misses: " << values[0] << ", L2 cache misses: " << values[1] << "\n"
+        //          << "L1 cache hits: " << values[2] << ", L2 cache hits: " << values[3] << "\n";
     }
 
     PAPI_shutdown();

modules/cuqdyn-c/include/states_transformer.h

@@ -10,4 +10,6 @@
 /// transformation expr is defined in the config XML.
 SUNMatrix transform_states(SUNMatrix states);
 
+SUNMatrix transform_states_2(SUNMatrix states);
+
 #endif //STATES_TRANSFORMER_H

modules/cuqdyn-c/src/states_transformer.c

@@ -2,10 +2,13 @@
 // Created by borja on 12/06/25.
 //
 #include <nvector/nvector_serial.h>
+#include <string.h>
 #include <sunmatrix/sunmatrix_dense.h>
 
 #include "config.h"
 #include "matlab.h"
+#include "sundials/sundials_matrix.h"
+#include "sundials/sundials_types.h"
 
 extern void eval_states_transformer_expr(sunrealtype *input, sunrealtype *output, CuqDynContext context);
 
@@ -20,10 +23,10 @@ SUNMatrix transform_states(SUNMatrix states)
 
     const int rows = SM_ROWS_D(states);
     const int cols = SM_COLUMNS_D(states);
-    
+
     SUNMatrix transformed_result = NewDenseMatrix(rows, conf->states_transformer.count + 1); // + 1 adds the time column
     sunrealtype *output = malloc(conf->states_transformer.count * sizeof(sunrealtype));
-    
+
     for (int i = 0; i < rows; ++i)
     {
         N_Vector input = copy_matrix_row(states, i, 1, cols);
@@ -38,8 +41,38 @@ SUNMatrix transform_states(SUNMatrix states)
 
         N_VDestroy(input);
     }
-    
+
     free(output);
     SUNMatDestroy(states);
     return transformed_result;
-}
+}
+
+SUNMatrix transform_states_2(SUNMatrix transposed_states)
+{
+    CuqdynConf *conf = get_cuqdyn_conf(get_cuqdyn_context());
+
+    if (conf->states_transformer.count == 0)
+    {
+        return transposed_states;
+    }
+
+    const int rows = SM_ROWS_D(transposed_states);
+    const int cols = SM_COLUMNS_D(transposed_states);
+
+    SUNMatrix transposed_result = NewDenseMatrix(conf->states_transformer.count + 1, rows);
+
+    for (int i = 0; i < rows; ++i)
+    {
+        // The first element is the time point
+        sunrealtype *input = &SM_COLUMN_D(transposed_states, i)[1];
+        sunrealtype *dest = &SM_COLUMN_D(transposed_result, i)[1];
+
+        eval_states_transformer_expr(input, dest, get_cuqdyn_context());
+
+        dest[0] = input[0]; // Copy the time point to the first element of the transformed state
+    }
+
+    SUNMatDestroy(transposed_states);
+
+    return transposed_result;
+}