Merge branch 'stclib:master' into modifiedHEADmodified

5 files changed, 302 insertions, 0 deletions

diff --git a/misc/examples/spans/matmult.c b/misc/examples/spans/matmult.c
new file mode 100644
index 00000000..ec992ff9
--- /dev/null
+++ b/misc/examples/spans/matmult.c
@@ -0,0 +1,90 @@
+// https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2023/p2642r2.html
+// C99:
+#include <stdio.h>
+#include <time.h>
+#include <stc/cspan.h>
+
+using_cspan3(Mat, double);
+typedef Mat2 OutMat;
+typedef struct { Mat2 m00, m01, m10, m11; } Partition;
+
+Partition partition(Mat2 A)
+{
+  int32_t M = A.shape[0];
+  int32_t N = A.shape[1];
+  return (Partition){
+    .m00 = cspan_slice(Mat2, &A, {0, M/2}, {0, N/2}),
+    .m01 = cspan_slice(Mat2, &A, {0, M/2}, {N/2, N}),
+    .m10 = cspan_slice(Mat2, &A, {M/2, M}, {0, N/2}),
+    .m11 = cspan_slice(Mat2, &A, {M/2, M}, {N/2, N}),
+  };
+}
+
+// Slow generic implementation
+void base_case_matrix_product(Mat2 A, Mat2 B, OutMat C)
+{
+  for (int j = 0; j < C.shape[1]; ++j) {
+    for (int i = 0; i < C.shape[0]; ++i) {
+      Mat2_value C_ij = 0;
+      for (int k = 0; k < A.shape[1]; ++k) {
+        C_ij += *cspan_at(&A, i,k) * *cspan_at(&B, k,j);
+      }
+      *cspan_at(&C, i,j) += C_ij;
+    }
+  }
+}
+
+void recursive_matrix_product(Mat2 A, Mat2 B, OutMat C)
+{
+  // Some hardware-dependent constant
+  enum {recursion_threshold = 32};
+  if (C.shape[0] <= recursion_threshold || C.shape[1] <= recursion_threshold) {
+    base_case_matrix_product(A, B, C);
+  } else {
+    Partition c = partition(C),
+              a = partition(A),
+              b = partition(B);
+    recursive_matrix_product(a.m00, b.m00, c.m00);
+    recursive_matrix_product(a.m01, b.m10, c.m00);
+    recursive_matrix_product(a.m10, b.m00, c.m10);
+    recursive_matrix_product(a.m11, b.m10, c.m10);
+    recursive_matrix_product(a.m00, b.m01, c.m01);
+    recursive_matrix_product(a.m01, b.m11, c.m01);
+    recursive_matrix_product(a.m10, b.m01, c.m11);
+    recursive_matrix_product(a.m11, b.m11, c.m11);
+  }
+}
+
+
+#define i_type Values
+#define i_val double
+#include <stc/cstack.h>
+#include <stc/crand.h>
+
+int main(void)
+{
+  enum {N = 10, D = 256};
+
+  Values values = {0};
+  for (int i=0; i < N*D*D; ++i)
+      Values_push(&values, (crandf() - 0.5)*4.0);
+
+  double out[D*D];
+  Mat3 data = cspan_md_layout(c_ROWMAJOR, values.data, N, D, D);
+  OutMat c = cspan_md_layout(c_COLMAJOR, out, D, D);
+  Mat2 a = cspan_submd3(&data, 0);
+
+  clock_t t = clock();
+  for (int i=1; i<N; ++i) {
+    Mat2 b = cspan_submd3(&data, i);
+    memset(out, 0, sizeof out);
+    recursive_matrix_product(a, b, c);
+    //base_case_matrix_product(a, b, c);
+  }
+  t = clock() - t;
+
+  double sum = 0.0;
+  c_foreach (i, Mat2, c) sum += *i.ref;
+  printf("sum=%.16g, %f ms\n", sum, (double)t*1000.0/CLOCKS_PER_SEC);
+  Values_drop(&values);
+}
diff --git a/misc/examples/spans/mdspan.c b/misc/examples/spans/mdspan.c
new file mode 100644
index 00000000..630ffddb
--- /dev/null
+++ b/misc/examples/spans/mdspan.c
@@ -0,0 +1,51 @@
+#include <stdio.h>
+#include <stc/cspan.h>
+#include <stdlib.h>
+
+using_cspan3(DSpan, double);
+
+int main(void) {
+    const int nx=5, ny=4, nz=3;
+    double* data = c_new_n(double, nx*ny*nz);
+
+    printf("\nMultidim span ms[5, 4, 3], fortran ordered");
+    DSpan3 ms = cspan_md_layout(c_COLMAJOR, data, nx, ny, nz);
+
+    int idx = 0;
+    for (int i = 0; i < ms.shape[0]; ++i)
+        for (int j = 0; j < ms.shape[1]; ++j)
+            for (int k = 0; k < ms.shape[2]; ++k)
+                *cspan_at(&ms, i, j, k) = ++idx;
+
+    cspan_transpose(&ms);
+
+    printf(", transposed:\n\n");
+    for (int i = 0; i < ms.shape[0]; ++i) {
+        for (int j = 0; j < ms.shape[1]; ++j) {
+            for (int k = 0; k < ms.shape[2]; ++k)
+                printf(" %3g", *cspan_at(&ms, i, j, k));
+            puts("");
+        }
+        puts("");
+    }
+
+    DSpan2 sub;
+
+    puts("Slicing:");
+    printf("\nms[0, :, :] ");
+    sub = cspan_slice(DSpan2, &ms, {0}, {c_ALL}, {c_ALL});
+    c_foreach (i, DSpan2, sub) printf(" %g", *i.ref);
+    puts("");
+
+    printf("\nms[:, 0, :] ");
+    sub = cspan_slice(DSpan2, &ms, {c_ALL}, {0}, {c_ALL});
+    c_foreach (i, DSpan2, sub) printf(" %g", *i.ref);
+    puts("");
+
+    sub = cspan_slice(DSpan2, &ms, {c_ALL}, {c_ALL}, {0});
+    printf("\nms[:, :, 0] ");
+    c_foreach (i, DSpan2, sub) printf(" %g", *i.ref);
+    puts("");
+
+    free(data);
+}
diff --git a/misc/examples/spans/multidim.c b/misc/examples/spans/multidim.c
new file mode 100644
index 00000000..70fda7e2
--- /dev/null
+++ b/misc/examples/spans/multidim.c
@@ -0,0 +1,76 @@
+// Example based on https://en.cppreference.com/w/cpp/container/mdspan
+#define i_val int
+#include <stc/cstack.h>
+#define i_implement
+#include <stc/cspan.h>
+#include <stdio.h>
+
+using_cspan3(ispan, int);
+
+void print2d(ispan2 ms2) {
+    for (int i=0; i < ms2.shape[0]; i++) {
+        for (int j=0; j < ms2.shape[1]; j++)
+            printf(" %3d", *cspan_at(&ms2, i, j));
+        puts("");
+    }
+}
+
+void print3d(ispan3 ms3) {
+    for (int i=0; i < ms3.shape[0]; i++) {
+        for (int j=0; j < ms3.shape[1]; j++) {
+            for (int k=0; k < ms3.shape[2]; k++)
+                printf(" %3d", *cspan_at(&ms3, i, j, k));
+            puts("");
+        }
+        puts("");
+    }
+}
+
+int main(void)
+{
+    cstack_int v = c_init(cstack_int, {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24});
+
+    // Create 1d span from a compatibel container
+    ispan ms1 = cspan_from(&v);
+
+    // Create a 3D mdspan 2 x 3 x 4
+    ispan3 ms3 = cspan_md(v.data, 2, 3, 4);
+
+    puts("ms3:");
+    print3d(ms3);
+
+    // Take a slice of md3
+    ispan3 ss3 = cspan_slice(ispan3, &ms3, {c_ALL}, {1,3}, {1,3});
+    puts("ss3 = ms3[:, 1:3, 1:3]");
+    print3d(ss3);
+
+    puts("Iterate ss3 flat:");
+    c_foreach (i, ispan3, ss3) printf(" %d", *i.ref);
+    puts("");
+
+    // submd3 span reduces rank depending on number of arguments
+    ispan2 ms2 = cspan_submd3(&ms3, 1);
+
+    // Change data on the 2d subspan
+    for (int i=0; i != ms2.shape[0]; i++)
+        for (int j=0; j != ms2.shape[1]; j++)
+            *cspan_at(&ms2, i, j) = (i + 1)*100 + j;
+    
+    puts("\nms2 = ms3[1] with updated data:");    
+    print2d(ms2);
+    puts("");    
+
+    puts("\nOriginal s1 span with updated data:");
+    c_foreach (i, ispan, ms1) printf(" %d", *i.ref);
+    puts("");
+
+    puts("\nOriginal ms3 span with updated data:");
+    print3d(ms3);
+
+    puts("col = ms3[1, :, 2]");
+    ispan col = cspan_slice(ispan, &ms3, {1}, {c_ALL}, {2});
+    c_foreach (i, ispan, col) printf(" %d", *i.ref);
+    puts("");
+
+    cstack_int_drop(&v);
+}
diff --git a/misc/examples/spans/printspan.c b/misc/examples/spans/printspan.c
new file mode 100644
index 00000000..b6999b61
--- /dev/null
+++ b/misc/examples/spans/printspan.c
@@ -0,0 +1,41 @@
+// https://www.modernescpp.com/index.php/c-20-std-span/
+
+#include <stdio.h>
+#define i_key int
+#include <stc/cvec.h>
+
+#define i_key int
+#include <stc/cstack.h>
+
+#include <stc/cspan.h>
+using_cspan(intspan, int);
+
+
+void printMe(intspan container) {
+    printf("%d:", (int)cspan_size(&container));
+    c_foreach (e, intspan, container)
+        printf(" %d", *e.ref);
+    puts("");
+}
+
+
+int main(void)
+{
+    printMe( c_init(intspan, {1, 2, 3, 4}) );
+
+    int arr[] = {1, 2, 3, 4, 5};
+    printMe( (intspan)cspan_from_array(arr) );
+
+    cvec_int vec = c_init(cvec_int, {1, 2, 3, 4, 5, 6});
+    printMe( (intspan)cspan_from(&vec) );
+
+    cstack_int stk = c_init(cstack_int, {1, 2, 3, 4, 5, 6, 7});
+    printMe( (intspan)cspan_from(&stk) );
+
+    intspan spn = c_init(intspan, {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12});
+    printMe( (intspan)cspan_subspan(&spn, 2, 8) );
+
+    // cleanup
+    cvec_int_drop(&vec);
+    cstack_int_drop(&stk);
+}
diff --git a/misc/examples/spans/submdspan.c b/misc/examples/spans/submdspan.c
new file mode 100644
index 00000000..0752dfa1
--- /dev/null
+++ b/misc/examples/spans/submdspan.c
@@ -0,0 +1,44 @@
+// https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2022/p2630r0.html
+// C99:
+#include <stdio.h>
+#include <stc/cspan.h>
+
+using_cspan3(span, double); // shorthand for defining span, span2, span3
+
+// Set all elements of a rank-2 mdspan to zero.
+void zero_2d(span2 grid2d) {
+  (void)c_static_assert(cspan_rank(&grid2d) == 2);
+  for (int i = 0; i < grid2d.shape[0]; ++i) {
+    for (int j = 0; j < grid2d.shape[1]; ++j) {
+      *cspan_at(&grid2d, i,j) = 0;
+    }
+  }
+}
+
+void zero_surface(span3 grid3d) {
+  (void)c_static_assert(cspan_rank(&grid3d) == 3);
+  zero_2d(cspan_slice(span2, &grid3d, {0}, {c_ALL}, {c_ALL}));
+  zero_2d(cspan_slice(span2, &grid3d, {c_ALL}, {0}, {c_ALL}));
+  zero_2d(cspan_slice(span2, &grid3d, {c_ALL}, {c_ALL}, {0}));
+  zero_2d(cspan_slice(span2, &grid3d, {grid3d.shape[0]-1}, {c_ALL}, {c_ALL}));
+  zero_2d(cspan_slice(span2, &grid3d, {c_ALL}, {grid3d.shape[1]-1}, {c_ALL}));
+  zero_2d(cspan_slice(span2, &grid3d, {c_ALL}, {c_ALL}, {grid3d.shape[2]-1}));
+}
+
+int main() {
+  double arr[3*4*5];
+  for (int i=0; i<c_arraylen(arr); ++i) arr[i] = i + 1.0;
+
+  span3 md = cspan_md(arr, 3, 4, 5);
+
+  zero_surface(md);
+
+  for (int i = 0; i < md.shape[0]; i++) {
+    for (int j = 0; j < md.shape[1]; j++) {
+      for (int k = 0; k < md.shape[2]; k++)
+        printf(" %2g", *cspan_at(&md, i,j,k));
+      puts("");
+    }
+    puts("");
+  }
+}