-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathmat.h
More file actions
9085 lines (7691 loc) · 292 KB
/
Copy pathmat.h
File metadata and controls
9085 lines (7691 loc) · 292 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#ifndef MAT_H_
#define MAT_H_
#define MAT_VERSION_MAJOR 0
#define MAT_VERSION_MINOR 1
#define MAT_VERSION_PATCH 1
#define MAT_VERSION "0.1.1"
#include <stdbool.h>
#include <stddef.h>
#include <stdint.h>
// Column-major storage is always enabled (required for optimal BLAS performance)
#define MAT_COLUMN_MAJOR
// ARM NEON detection
#if defined(__ARM_NEON) || defined(__ARM_NEON__) || defined(_M_ARM64)
#define MAT_HAS_ARM_NEON
#include <arm_neon.h>
#endif
// x86 AVX2 detection
#if defined(__AVX2__)
#define MAT_HAS_AVX2
#include <immintrin.h>
#endif
// Architecture indicator macro - shows which SIMD backend is active
// Usage: printf("libmat using: %s\n", MAT_ARCH);
#if defined(MAT_HAS_ARM_NEON)
#define MAT_ARCH "NEON"
#elif defined(MAT_HAS_AVX2)
#define MAT_ARCH "AVX2"
#else
#define MAT_ARCH "SCALAR"
#endif
/* SIMD Dispatch Architecture */
// OpenMP detection (auto-enabled when compiled with -fopenmp)
#if defined(_OPENMP)
#include <omp.h>
#define MAT_HAS_OPENMP
#ifndef MAT_OMP_THRESHOLD
#define MAT_OMP_THRESHOLD \
(1024 * 1024) // Skip parallelization for small matrices
#endif
#endif
#ifndef MATDEF
#define MATDEF
#endif
// Mark experimental functions (warns on use)
#if defined(__GNUC__)
#define MAT_EXPERIMENTAL \
__attribute__((warning("experimental: does not guarantee correct results")))
#elif defined(_MSC_VER)
#define MAT_EXPERIMENTAL \
__declspec(deprecated("experimental: does not guarantee correct results"))
#else
#define MAT_EXPERIMENTAL
#endif
// Mark unimplemented functions (errors on use)
#if defined(__GNUC__)
#define MAT_NOT_IMPLEMENTED __attribute__((error("not implemented")))
#elif defined(_MSC_VER)
#define MAT_NOT_IMPLEMENTED \
__declspec(deprecated("not implemented - will not link"))
#else
#define MAT_NOT_IMPLEMENTED
#endif
// Overridable allocator macros
// Define these before including mat.h to use custom allocators
// (e.g., arenas)
#ifndef MAT_MALLOC
#define MAT_MALLOC(sz) malloc(sz)
#endif
#ifndef MAT_CALLOC
#define MAT_CALLOC(n, sz) calloc(n, sz)
#endif
#ifndef MAT_FREE
#define MAT_FREE(p) free(p)
#endif
#ifndef MAT_FREE_MAT
#define MAT_FREE_MAT(m) \
do { \
MAT_FREE((m)->data); \
MAT_FREE(m); \
} while (0)
#endif
#ifndef MAT_FREE_PERM
#define MAT_FREE_PERM(p) \
do { \
MAT_FREE((p)->data); \
MAT_FREE(p); \
} while (0)
#endif
// Scratch arena for temporary allocations in hot paths
// Define MAT_NO_SCRATCH to disable (uses malloc/free instead)
// Define MAT_SCRATCH_SIZE to override the default size
#ifndef MAT_SCRATCH_SIZE
#define MAT_SCRATCH_SIZE (4 * 1024 * 1024) // 4MB default
#endif
#ifndef MAT_NO_SCRATCH
typedef struct {
char *buf;
size_t offset;
size_t size;
} MatArena;
#endif
// Transpose flags for GEMM operations (BLAS-style)
typedef enum {
MAT_NO_TRANS = 0, // Use matrix as-is
MAT_TRANS = 1 // Use transpose of matrix
} mat_trans_t;
#ifdef MAT_STRIP_PREFIX
#define mat mat_mat
#define empty mat_empty
#define zeros mat_zeros
#define ones mat_ones
#define eye mat_reye
#define reye mat_reye
#define deep_copy mat_deep_copy
#define rdeep_copy mat_rdeep_copy
#define t mat_t
#define rt mat_rt
#define reshape mat_reshape
#define rreshape mat_rreshape
#define diag mat_diag
#define diag_from mat_diag_from
#define vec_from mat_vec_from
#define free_mat mat_free_mat
#define hadamard mat_hadamard
#define rhadamard mat_rhadamard
#define add_scalar mat_add_scalar
#define radd_scalar mat_radd_scalar
#define add_many mat_add_many
#define radd_many mat_radd_many
#endif // MAT_STRIP_PREFIX
// Initialization macros
#define mat_new(cols, ...) \
mat_from(sizeof((mat_elem_t[][cols])__VA_ARGS__) / sizeof(mat_elem_t[cols]), \
cols, (mat_elem_t *)((mat_elem_t[][cols])__VA_ARGS__))
#define mat_set(out, ...) mat_init(out, (mat_elem_t[])__VA_ARGS__)
#define mat_vnew(...) \
mat_vec_from(sizeof((mat_elem_t[])__VA_ARGS__) / sizeof(mat_elem_t), \
(mat_elem_t[])__VA_ARGS__)
#define mat_rnew(...) \
mat_from(1, sizeof((mat_elem_t[])__VA_ARGS__) / sizeof(mat_elem_t), \
(mat_elem_t[])__VA_ARGS__)
// Element type (float or double precision)
#ifdef MAT_DOUBLE_PRECISION
typedef double mat_elem_t;
#ifndef MAT_DEFAULT_EPSILON
#define MAT_DEFAULT_EPSILON 1e-9
#endif
#define MAT_FABS fabs
#define MAT_SQRT sqrt
#define MAT_HUGE_VAL HUGE_VAL
#else
typedef float mat_elem_t;
#ifndef MAT_DEFAULT_EPSILON
#define MAT_DEFAULT_EPSILON 1e-6f
#endif
#define MAT_FABS fabsf
#define MAT_SQRT sqrtf
#define MAT_HUGE_VAL HUGE_VALF
#endif
// NEON SIMD macros (only defined when targeting ARM with NEON)
#ifdef MAT_HAS_ARM_NEON
#ifdef MAT_DOUBLE_PRECISION
// NEON double precision (2 doubles per 128-bit register)
#define MAT_NEON_TYPE float64x2_t
#define MAT_NEON_UTYPE uint64x2_t
#define MAT_NEON_WIDTH 2
#define MAT_NEON_LOAD vld1q_f64
#define MAT_NEON_STORE vst1q_f64
#define MAT_NEON_DUP vdupq_n_f64
#define MAT_NEON_DUP_U vdupq_n_u64
#define MAT_NEON_FMA vfmaq_f64
#define MAT_NEON_FMA_LANE vfmaq_laneq_f64
#define MAT_NEON_FMS vfmsq_f64
#define MAT_NEON_ADD vaddq_f64
#define MAT_NEON_ADDV vaddvq_f64
#define MAT_NEON_ABS vabsq_f64
#define MAT_NEON_MAX vmaxq_f64
#define MAT_NEON_MAXV vmaxvq_f64
#define MAT_NEON_ABD vabdq_f64
#define MAT_NEON_CGT vcgtq_f64
#define MAT_NEON_CEQ vceqq_f64
#define MAT_NEON_ORR_U vorrq_u64
#define MAT_NEON_AND_U vandq_u64
#define MAT_NEON_MVN_U(x) veorq_u64(x, vdupq_n_u64(0xFFFFFFFFFFFFFFFFULL))
#define MAT_NEON_MAXV_U(x) (vgetq_lane_u64(x, 0) | vgetq_lane_u64(x, 1))
#define MAT_NEON_ADDV_U(x) (vgetq_lane_u64(x, 0) + vgetq_lane_u64(x, 1))
#define MAT_NEON_ADD_U vaddq_u64
#define MAT_NEON_MUL vmulq_f64
#define MAT_NEON_SUB vsubq_f64
#define MAT_NEON_ZIP1 vzip1q_f64
#define MAT_NEON_ZIP2 vzip2q_f64
#define MAT_NEON_GET_LANE(v, n) vgetq_lane_f64(v, n)
#else
// NEON single precision (4 floats per 128-bit register)
#define MAT_NEON_TYPE float32x4_t
#define MAT_NEON_UTYPE uint32x4_t
#define MAT_NEON_WIDTH 4
#define MAT_NEON_LOAD vld1q_f32
#define MAT_NEON_STORE vst1q_f32
#define MAT_NEON_DUP vdupq_n_f32
#define MAT_NEON_DUP_U vdupq_n_u32
#define MAT_NEON_FMA vfmaq_f32
#define MAT_NEON_FMA_LANE vfmaq_laneq_f32
#define MAT_NEON_FMS vfmsq_f32
#define MAT_NEON_ADD vaddq_f32
#define MAT_NEON_ADDV vaddvq_f32
#define MAT_NEON_ABS vabsq_f32
#define MAT_NEON_MAX vmaxq_f32
#define MAT_NEON_MAXV vmaxvq_f32
#define MAT_NEON_ABD vabdq_f32
#define MAT_NEON_CGT vcgtq_f32
#define MAT_NEON_CEQ vceqq_f32
#define MAT_NEON_ORR_U vorrq_u32
#define MAT_NEON_AND_U vandq_u32
#define MAT_NEON_MVN_U(x) vmvnq_u32(x)
#define MAT_NEON_MAXV_U(x) vmaxvq_u32(x)
#define MAT_NEON_ADDV_U(x) vaddvq_u32(x)
#define MAT_NEON_ADD_U vaddq_u32
#define MAT_NEON_MUL vmulq_f32
#define MAT_NEON_SUB vsubq_f32
#define MAT_NEON_ZIP1 vzip1q_f32
#define MAT_NEON_ZIP2 vzip2q_f32
#define MAT_NEON_GET_LANE(v, n) vgetq_lane_f32(v, n)
#endif
// Accumulator macros for overflow-safe reductions (always use double)
#define MAT_ACC_TYPE float64x2_t
#define MAT_ACC_ZERO vdupq_n_f64(0)
#define MAT_ACC_ADD vaddq_f64
#define MAT_ACC_ADDV vaddvq_f64
#define MAT_ACC_FMA vfmaq_f64
#ifdef MAT_DOUBLE_PRECISION
#define MAT_ACC_WIDTH 2
#define MAT_ACC_LOAD_SQ(acc, ptr) \
do { \
float64x2_t _v = vld1q_f64(ptr); \
acc = vfmaq_f64(acc, _v, _v); \
} while (0)
#else
#define MAT_ACC_WIDTH 2
#define MAT_ACC_LOAD_SQ(acc, ptr) \
do { \
float32x2_t _v = vld1_f32(ptr); \
float64x2_t _d = vcvt_f64_f32(_v); \
acc = vfmaq_f64(acc, _d, _d); \
} while (0)
#endif
#endif // __ARM_NEON
/* AVX2 SIMD Macros (TODO) */
#ifdef MAT_HAS_AVX2
// AVX2 double precision (4 doubles per 256-bit register)
// AVX2 single precision (8 floats per 256-bit register)
// TODO: Add AVX2 abstraction macros similar to NEON
// When implementing, use _mm256_* intrinsics from <immintrin.h>
#endif // MAT_HAS_AVX2
// Control visibility of internal implementations
#ifdef MAT_EXPOSE_INTERNALS
#define MAT_INTERNAL_STATIC
#else
#ifdef __GNUC__
#define MAT_INTERNAL_STATIC static __attribute__((unused))
#else
#define MAT_INTERNAL_STATIC static
#endif
#endif
#define identity reye
#ifndef MAT_LOG_LEVEL
#define MAT_LOG_LEVEL 0
#endif
// Logging (set MAT_LOG_LEVEL > 0 to enable)
#if MAT_LOG_LEVEL > 0
#include <stdio.h>
#ifndef MAT_LOG_OUTPUT_ERR
#define MAT_LOG_OUTPUT_ERR(msg) fprintf(stderr, "%s\n", msg)
#endif
#ifndef MAT_LOG_OUTPUT
#define MAT_LOG_OUTPUT(msg) fprintf(stdout, "%s\n", msg)
#endif
#endif
#if MAT_LOG_LEVEL >= 1
#define MAT_LOG_ERROR(msg) MAT_LOG_OUTPUT_ERR("[ERROR] " msg)
#else
#define MAT_LOG_ERROR(msg)
#endif
#if MAT_LOG_LEVEL >= 2
#define MAT_LOG_WARN(msg) MAT_LOG_OUTPUT("[WARN] " msg)
#else
#define MAT_LOG_WARN(msg)
#endif
#if MAT_LOG_LEVEL >= 3
#define MAT_LOG_INFO(msg) MAT_LOG_OUTPUT("[INFO] " msg)
#else
#define MAT_LOG_INFO(msg)
#endif
// Assertions (define MAT_ASSERT to override)
#ifndef MAT_ASSERT
#include <assert.h>
#define MAT_ASSERT(x) assert(x)
#endif
// Type-specific assertion macros
#ifndef MAT_ASSERT_MAT
#define MAT_ASSERT_MAT(m) \
do { \
MAT_ASSERT((m) != NULL); \
MAT_ASSERT((m)->data != NULL); \
MAT_ASSERT((m)->rows > 0 && (m)->cols > 0); \
} while (0)
#endif
#ifndef MAT_ASSERT_DIM
#define MAT_ASSERT_DIM(rows, cols) \
do { \
MAT_ASSERT((rows) > 0); \
MAT_ASSERT((cols) > 0); \
} while (0)
#endif
#ifndef MAT_ASSERT_SQUARE
#define MAT_ASSERT_SQUARE(m) \
do { \
MAT_ASSERT_MAT(m); \
MAT_ASSERT((m)->rows == (m)->cols); \
} while (0)
#endif
#ifdef __cplusplus
extern "C" {
#endif
// Core matrix type
typedef struct {
size_t rows;
size_t cols;
mat_elem_t *data;
} Mat;
// Matrix dimensions tuple
typedef struct {
size_t x;
size_t y;
} MatSize;
// Vectors in this library are implemented as a Mat
// of size nx1 (column vector) by default.
// Row vectors are just transposed column vectors.
typedef Mat Vec;
// Permutation vector (array of indices)
typedef struct {
size_t *data;
size_t size;
} Perm;
// Storage order: Column-major (Fortran/BLAS style, column-contiguous)
// data[j * rows + i] accesses element at row i, column j
#define MAT_IDX(m, i, j) ((j) * (m)->rows + (i))
// Element access macros - use these instead of direct indexing
#define MAT_AT(m, i, j) ((m)->data[MAT_IDX(m, i, j)])
#define MAT_SET(m, i, j, v) ((m)->data[MAT_IDX(m, i, j)] = (v))
/* Construction & Memory */
// Allocate matrix struct without data buffer. Use for custom memory management.
MATDEF Mat *mat_empty(size_t rows, size_t cols);
// Allocate zero-initialized matrix.
MATDEF Mat *mat_mat(size_t rows, size_t cols);
// Create matrix from array of values (row-major order).
MATDEF Mat *mat_from(size_t rows, size_t cols, const mat_elem_t *values);
// Initialize existing matrix with values (row-major order).
MATDEF void mat_init(Mat *out, const mat_elem_t *values);
// Free matrix and its data buffer.
MATDEF void mat_free_mat(Mat *m);
// Create matrix filled with zeros.
MATDEF Mat *mat_zeros(size_t rows, size_t cols);
// Create matrix filled with ones.
MATDEF Mat *mat_ones(size_t rows, size_t cols);
// Create a matrix filled with any given number
MATDEF void mat_fill(Mat *out, mat_elem_t value);
// Set matrix to identity (must be square).
MATDEF void mat_eye(Mat *out);
// Create identity matrix of given dimension.
MATDEF Mat *mat_reye(size_t dim);
// Allocate zero-initialized column vector.
MATDEF Vec *mat_vec(size_t dim);
// Allocate a zero-initialized row vector.
MATDEF Vec *mat_row_vec(size_t dim);
// Create column vector from array of values.
MATDEF Vec *mat_vec_from(size_t dim, const mat_elem_t *values);
// Allocate permutation of given size.
MATDEF Perm *mat_perm(size_t n);
// Free permutation and its data.
MATDEF void mat_free_perm(Perm *p);
// Set permutation to identity [0, 1, 2, ...].
MATDEF void mat_perm_identity(Perm *p);
// Convert permutation to explicit n×n permutation matrix.
MATDEF Mat *mat_perm_mat(const Perm *p);
// Shallow copy (copies struct, shares data pointer). Use mat_rdeep_copy for
// full copy.
MATDEF Mat *mat_copy(const Mat *m);
// Deep copy src into pre-allocated out.
MATDEF void mat_deep_copy(Mat *out, const Mat *src);
// Allocate and return deep copy.
MATDEF Mat *mat_rdeep_copy(const Mat *m);
/* Accessors & Info */
// Get element at (row, col). Bounds checked via MAT_ASSERT.
MATDEF mat_elem_t mat_at(const Mat *mat, size_t row, size_t col);
// Set element at (row, col). Bounds checked via MAT_ASSERT.
MATDEF void mat_set_at(Mat *m, size_t row, size_t col, mat_elem_t value);
// Get matrix dimensions as {rows, cols} struct.
MATDEF MatSize mat_size(const Mat *m);
// Print matrix to stdout in MATLAB-like format.
MATDEF void mat_print(const Mat *m);
/* Comparison */
// Exact equality. Returns true if all elements are bitwise equal.
MATDEF bool mat_equals(const Mat *a, const Mat *b);
// Approximate equality. Returns true if max|a-b| < epsilon.
MATDEF bool mat_equals_tol(const Mat *a, const Mat *b, mat_elem_t epsilon);
/* Element-wise Unary */
// NOTE: Transcendental functions (exp, log, sin, cos, atan2) use scalar libc
// calls. For SIMD-optimized versions, consider libraries like SLEEF or Eigen.
// out[i] = |a[i]|
MATDEF void mat_abs(Mat *out, const Mat *a);
// out[i] = sqrt(a[i]). Undefined for negative inputs.
MATDEF void mat_sqrt(Mat *out, const Mat *a);
// out[i] = e^a[i]
MATDEF void mat_exp(Mat *out, const Mat *a);
// out[i] = ln(a[i]). Undefined for non-positive inputs.
MATDEF void mat_log(Mat *out, const Mat *a);
// out[i] = log10(a[i]). Undefined for non-positive inputs.
MATDEF void mat_log10(Mat *out, const Mat *a);
// out[i] = sin(a[i]), a in radians.
MATDEF void mat_sin(Mat *out, const Mat *a);
// out[i] = cos(a[i]), a in radians.
MATDEF void mat_cos(Mat *out, const Mat *a);
// out[i] = a[i]^exp. Undefined for negative base with fractional exponent.
MATDEF void mat_pow(Mat *out, const Mat *a, mat_elem_t exp);
// out[i] = clamp(a[i], min_val, max_val).
MATDEF void mat_clip(Mat *out, const Mat *a, mat_elem_t min_val,
mat_elem_t max_val);
// Element-wise Binary
// out[i] = a[i] / b[i]. No division-by-zero check.
MATDEF void mat_div(Mat *out, const Mat *a, const Mat *b);
// out[i] = atan2(y[i], x[i]). Result in radians, range [-pi, pi].
MATDEF void mat_atan2(Mat *out, const Mat *y, const Mat *x);
/* Scalar Operations */
// out[i] *= k. In-place scaling. SIMD-optimized.
MATDEF void mat_scale(Mat *out, mat_elem_t k);
// Return m scaled by k. Allocates new matrix.
MATDEF Mat *mat_rscale(const Mat *m, mat_elem_t k);
// Normalize vector in place: v = v / ||v||. Returns the norm.
// If norm < MAT_DEFAULT_EPSILON, vector is unchanged and 0 is returned.
// SIMD-optimized.
MATDEF mat_elem_t mat_normalize(Mat *v);
// out[i] += k. In-place scalar addition.
MATDEF void mat_add_scalar(Mat *out, mat_elem_t k);
// Return m + k. Allocates new matrix.
MATDEF Mat *mat_radd_scalar(const Mat *m, mat_elem_t k);
// Matrix Arithmetic
// out = a + b. SIMD-optimized.
MATDEF void mat_add(Mat *out, const Mat *a, const Mat *b);
// Return a + b. Allocates new matrix.
MATDEF Mat *mat_radd(const Mat *a, const Mat *b);
// out = a - b. SIMD-optimized.
MATDEF void mat_sub(Mat *out, const Mat *a, const Mat *b);
// Return a - b. Allocates new matrix.
MATDEF Mat *mat_rsub(const Mat *a, const Mat *b);
// out = sum of count matrices (variadic). Modifies out in-place.
MATDEF void mat_add_many(Mat *out, size_t count, ...);
// Return sum of count matrices (variadic). Allocates new matrix.
MATDEF Mat *mat_radd_many(size_t count, ...);
// Matrix Products
// out = a * b (matrix multiplication). SIMD-optimized.
// Dimensions: a(m,k) * b(k,n) = out(m,n).
MATDEF void mat_mul(Mat *out, const Mat *a, const Mat *b);
// Return a * b. Allocates new matrix.
MATDEF Mat *mat_rmul(const Mat *a, const Mat *b);
// out[i] = a[i] * b[i] (element-wise/Hadamard product). SIMD-optimized.
MATDEF void mat_hadamard(Mat *out, const Mat *a, const Mat *b);
// Return element-wise a * b. Allocates new matrix.
MATDEF Mat *mat_rhadamard(const Mat *a, const Mat *b);
// Return sum(v1[i] * v2[i]) (dot/inner product). SIMD-optimized.
MATDEF mat_elem_t mat_dot(const Vec *v1, const Vec *v2);
// out = v1 x v2 (cross product). Vectors must be 3D.
MATDEF void mat_cross(Vec *out, const Vec *v1, const Vec *v2);
// out = v1 * v2^T (outer product). out(m,n) where v1 is m-dim, v2 is n-dim.
// SIMD-optimized.
MATDEF void mat_outer(Mat *out, const Vec *v1, const Vec *v2);
// Return x^T * A * y (bilinear form). x is m-dim, A is m x n, y is n-dim.
MATDEF mat_elem_t mat_bilinear(const Vec *x, const Mat *A, const Vec *y);
// Return x^T * A * x (quadratic form). x is n-dim, A is n x n.
MATDEF mat_elem_t mat_quadform(const Vec *x, const Mat *A);
// Fused Operations (BLAS-like)
// y = alpha * x + y (AXPY). SIMD-optimized.
MATDEF void mat_axpy(Vec *y, mat_elem_t alpha, const Vec *x);
// y = alpha * A * x + beta * y (GEMV). SIMD-optimized.
MATDEF void mat_gemv(Vec *y, mat_elem_t alpha, const Mat *A, const Vec *x,
mat_elem_t beta);
// y = alpha * A^T * x + beta * y (GEMV transposed). SIMD-optimized.
MATDEF void mat_gemv_t(Vec *y, mat_elem_t alpha, const Mat *A, const Vec *x,
mat_elem_t beta);
// A = A + alpha * x * y^T (GER/rank-1 update).
MATDEF void mat_ger(Mat *A, mat_elem_t alpha, const Vec *x, const Vec *y);
// A = alpha * x * x^T + A (SYR). SIMD-optimized.
// uplo: 'L' for lower, 'U' for upper triangle.
MATDEF void mat_syr(Mat *A, mat_elem_t alpha, const Vec *x, char uplo);
// C = alpha * A * B + beta * C (GEMM). SIMD-optimized.
MATDEF void mat_gemm(Mat *C, mat_elem_t alpha, const Mat *A, const Mat *B,
mat_elem_t beta);
// C = alpha * A * A^T + beta * C (SYRK). SIMD-optimized.
// uplo: 'L' for lower, 'U' for upper triangle.
MATDEF void mat_syrk(Mat *C, const Mat *A, mat_elem_t alpha, mat_elem_t beta,
char uplo);
// C = alpha * A^T * A + beta * C (SYRK transposed). SIMD-optimized.
// uplo: 'L' for lower, 'U' for upper triangle.
MATDEF void mat_syrk_t(Mat *C, const Mat *A, mat_elem_t alpha, mat_elem_t beta,
char uplo);
// C = alpha * A * B^T + alpha * B * A^T + beta * C (SYR2K). SIMD-optimized.
// Symmetric rank-2k update. uplo: 'L' for lower, 'U' for upper triangle.
MATDEF void mat_syr2k(Mat *C, const Mat *A, const Mat *B, mat_elem_t alpha,
mat_elem_t beta, char uplo);
// Structure Operations
// out = m^T (transpose). out must be pre-allocated with swapped dimensions.
MATDEF void mat_t(Mat *out, const Mat *m);
// Return m^T. Allocates new matrix.
MATDEF Mat *mat_rt(const Mat *m);
// Reshape out in-place. Total elements must remain constant.
MATDEF void mat_reshape(Mat *out, size_t rows, size_t cols);
// Return reshaped copy of m. Allocates new matrix.
MATDEF Mat *mat_rreshape(const Mat *m, size_t rows, size_t cols);
// out = [a, b] (horizontal concatenation). a and b must have same row count.
MATDEF void mat_hcat(Mat *out, const Mat *a, const Mat *b);
// out = [a; b] (vertical concatenation). a and b must have same column count.
MATDEF void mat_vcat(Mat *out, const Mat *a, const Mat *b);
// Extract row as column vector. Allocates new vector.
MATDEF Vec *mat_row(const Mat *m, size_t row);
// Extract column as column vector. Allocates new vector.
MATDEF Vec *mat_col(const Mat *m, size_t col);
// Return a view of row as Vec (no allocation, no copy).
// The returned Vec points to data owned by m; do not free it.
MATDEF Vec mat_row_view(const Mat *m, size_t row);
// Extract submatrix m[row_start:row_end, col_start:col_end]. Allocates new
// matrix. Indices are inclusive start, exclusive end.
MATDEF Mat *mat_slice(const Mat *m, size_t row_start, size_t row_end,
size_t col_start, size_t col_end);
// Copy src into m starting at (row_start, col_start).
MATDEF void mat_slice_set(Mat *m, size_t row_start, size_t col_start,
const Mat *src);
// Diagonal Operations
// Extract main diagonal as column vector. Allocates new vector.
MATDEF Vec *mat_diag(const Mat *m);
// Create diagonal matrix from values. Returns dim x dim matrix.
MATDEF Mat *mat_diag_from(size_t dim, const mat_elem_t *values);
// Reduction Operations
// Return sum of all elements. SIMD-optimized.
MATDEF mat_elem_t mat_sum(const Mat *a);
// Return mean of all elements.
MATDEF mat_elem_t mat_mean(const Mat *a);
// Return minimum element value. SIMD-optimized.
MATDEF mat_elem_t mat_min(const Mat *a);
// Return maximum element value. SIMD-optimized.
MATDEF mat_elem_t mat_max(const Mat *a);
// Sum along axis. axis=0: sum columns (out has rows elements).
// axis=1: sum rows (out has cols elements).
MATDEF void mat_sum_axis(Vec *out, const Mat *a, int axis);
// Return flat index of minimum element.
MATDEF size_t mat_argmin(const Mat *a);
// Return flat index of maximum element.
MATDEF size_t mat_argmax(const Mat *a);
// Return population standard deviation.
MATDEF mat_elem_t mat_std(const Mat *a);
// Norms
// General p-norm: (sum |a_i|^p)^(1/p). Uses pow(), slow for large matrices.
MATDEF mat_elem_t mat_norm(const Mat *a, size_t p);
// L2 norm. Alias for mat_norm_fro.
MATDEF mat_elem_t mat_norm2(const Mat *a);
// Infinity norm: max |a_ij|. SIMD-optimized.
MATDEF mat_elem_t mat_norm_max(const Mat *a);
// Frobenius norm: sqrt(sum a_ij^2). SIMD-optimized.
// For float32: accumulates in double to prevent overflow/underflow.
// For float64: no overflow protection (same as fast). Blue's scaling may be
// added in the future to handle extreme values.
MATDEF mat_elem_t mat_norm_fro(const Mat *a);
// Frobenius norm, fast version. SIMD-optimized.
// For float32: ~2x faster than safe, but no overflow protection.
// For float64: same as safe (no higher precision available).
// Overflows if any |a_ij|^2 exceeds type max (~1e19 for float, ~1e154 for
// double).
MATDEF mat_elem_t mat_norm_fro_fast(const Mat *a);
// Matrix Properties
// Return trace (sum of diagonal elements). Matrix must be square.
MATDEF mat_elem_t mat_trace(const Mat *a);
// Return determinant. Matrix must be square. Uses LU decomposition.
MATDEF mat_elem_t mat_det(const Mat *a);
// Return count of non-zero elements.
MATDEF mat_elem_t mat_nnz(const Mat *a);
// Decomposition
// QR decomposition via Householder reflections.
// A = Q * R where Q is orthogonal (m x m), R is upper triangular (m x n).
// Q and R must be pre-allocated with correct dimensions.
MATDEF void mat_qr(const Mat *A, Mat *Q, Mat *R);
// QR decomposition (R factor only) - faster when Q is not needed.
// Useful for least squares, rank determination, etc.
// R must be pre-allocated with dimensions (m x n).
MATDEF void mat_qr_r(const Mat *A, Mat *R);
// Householder reflection: compute v and tau such that H*x = beta*e1
// where H = I - tau*v*v^T is orthogonal.
// v is modified in-place from x (v[0] = 1, rest normalized).
// Returns beta (the resulting first element after reflection).
MATDEF mat_elem_t mat_householder(Vec *v, mat_elem_t *tau, const Vec *x);
// Apply Householder reflection from left: A = H*A = A - tau*v*(v^T*A)
// v[0] is assumed to be 1 (implicit).
MATDEF void mat_householder_left(Mat *A, const Vec *v, mat_elem_t tau);
// Apply Householder reflection from right: A = A*H = A - tau*(A*v)*v^T
// v[0] is assumed to be 1 (implicit).
MATDEF void mat_householder_right(Mat *A, const Vec *v, mat_elem_t tau);
// LU decomposition with full pivoting.
// P * A * Q = L * U where P, Q are permutations (full pivoting).
// L is lower triangular with 1s on diagonal, U is upper triangular.
// L and U must be pre-allocated with dimensions n x n.
// P and Q must be pre-allocated permutations of size n.
// Returns the number of row+column swaps (useful for determinant sign).
MATDEF int mat_lu(const Mat *A, Mat *L, Mat *U, Perm *p, Perm *q);
// P * A = L * U where P is row permutation (partial pivoting).
// Faster than mat_lu, sufficient for determinant, solve, and inverse.
// Returns the number of row swaps (useful for determinant sign).
MATDEF int mat_plu(const Mat *A, Mat *L, Mat *U, Perm *p);
// Cholesky decomposition (A = L * L^T, A must be symmetric positive definite).
// Returns 0 on success, -1 if matrix is not positive definite.
MATDEF int mat_chol(const Mat *A, Mat *L);
// Singular value decomposition (A = U * S * Vt).
// Uses one-sided Jacobi algorithm.
// U is m x m orthogonal, S is min(m,n) vector of singular values (descending),
// Vt is n x n orthogonal (V transposed).
MATDEF void mat_svd(const Mat *A, Mat *U, Vec *S, Mat *Vt);
// Matrix inverse using LU decomposition.
MATDEF void mat_inv(Mat *out, const Mat *A);
// Moore-Penrose pseudoinverse via SVD.
// out must be n x m for input A of size m x n.
// Tolerance for rank determination: max(m,n) * max(S) * epsilon.
MATDEF void mat_pinv(Mat *out, const Mat *A);
// Matrix rank via SVD.
// Returns the number of singular values above tolerance.
// Tolerance: max(m,n) * max(S) * epsilon.
MATDEF size_t mat_rank(const Mat *A);
// Condition number via SVD.
// Returns sigma_max / sigma_min.
// For singular matrices, returns infinity.
MATDEF mat_elem_t mat_cond(const Mat *A);
// Eigendecomposition.
MAT_NOT_IMPLEMENTED MATDEF void mat_eig(const Mat *A, Vec *eigenvalues,
Mat *eigenvectors);
// Eigenvalues only (faster, no eigenvectors).
// Uses Hessenberg reduction followed by implicit QR iteration.
// Eigenvalues are returned in arbitrary order.
// For complex eigenvalues (conjugate pairs), the real part is stored.
MATDEF void mat_eigvals(Vec *out, const Mat *A);
// Eigenvalues of symmetric matrix (faster than mat_eigvals for symmetric input).
// Uses tridiagonal reduction + implicit QR iteration. O(n) per QR step vs O(n^2).
// Eigenvalues are returned in arbitrary order.
MATDEF void mat_eigvals_sym(Vec *out, const Mat *A);
// Solve Ax = b for x. A must be square and non-singular.
// Uses LU decomposition with partial pivoting.
MATDEF void mat_solve(Vec *x, const Mat *A, const Vec *b);
// Solve Ax = b for x where A is symmetric positive definite.
// Uses Cholesky decomposition (~2x faster than mat_solve for SPD matrices).
// Returns 0 on success, -1 if A is not positive definite.
MATDEF int mat_solve_spd(Vec *x, const Mat *A, const Vec *b);
// Triangular solvers (TRSV operations).
// Solve Lx = b where L is lower triangular.
MATDEF void mat_solve_tril(Vec *x, const Mat *L, const Vec *b);
// Solve Lx = b where L is unit lower triangular (implicit 1s on diagonal).
MATDEF void mat_solve_tril_unit(Vec *x, const Mat *L, const Vec *b);
// Solve Ux = b where U is upper triangular.
MATDEF void mat_solve_triu(Vec *x, const Mat *U, const Vec *b);
// Solve L^T x = b where L is lower triangular (uses L directly, no transpose).
MATDEF void mat_solve_trilt(Vec *x, const Mat *L, const Vec *b);
// Least squares solution.
MAT_NOT_IMPLEMENTED MATDEF void mat_lstsq(Vec *x, const Mat *A, const Vec *b);
// Kronecker product.
MAT_NOT_IMPLEMENTED MATDEF void mat_kron(Mat *out, const Mat *A, const Mat *B);
// 2D convolution.
MAT_NOT_IMPLEMENTED MATDEF void mat_conv2d(Mat *out, const Mat *A,
const Mat *kernel);
#ifdef __cplusplus
}
#endif
#endif // MAT_H_
#ifdef MAT_IMPLEMENTATION
#include <math.h>
#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
// Scratch arena for temporary allocations
#ifndef MAT_NO_SCRATCH
static MatArena mat_scratch_ = {0};
static inline void mat_scratch_init_(void) {
if (mat_scratch_.buf == NULL) {
mat_scratch_.buf = (char *)MAT_MALLOC(MAT_SCRATCH_SIZE);
mat_scratch_.size = MAT_SCRATCH_SIZE;
mat_scratch_.offset = 0;
}
}
static inline void *mat_scratch_alloc_(size_t bytes) {
mat_scratch_init_();
bytes = (bytes + 15) & ~15; // align to 16 bytes for NEON
if (mat_scratch_.offset + bytes > mat_scratch_.size) {
// Fall back to heap if arena is full
return MAT_MALLOC(bytes);
}
void *ptr = mat_scratch_.buf + mat_scratch_.offset;
mat_scratch_.offset += bytes;
return ptr;
}
static inline void mat_scratch_reset_(void) { mat_scratch_.offset = 0; }
#else
// No scratch arena - use malloc/free directly
static inline void *mat_scratch_alloc_(size_t bytes) {
return MAT_MALLOC(bytes);
}
static inline void mat_scratch_free_(void *ptr) { MAT_FREE(ptr); }
#endif
/* BLAS-like Kernels */
/* Kernel dispatch macro - selects architecture-specific implementation */
#ifdef MAT_HAS_ARM_NEON
#define MAT_DISPATCH(name) mat__##name##_neon
#elif defined(MAT_HAS_AVX2)
#define MAT_DISPATCH(name) mat__##name##_avx2
#else
#define MAT_DISPATCH(name) mat__##name##_scalar
#endif
// AXPY: y[0:n] += alpha * x[0:n]
MAT_INTERNAL_STATIC void mat__axpy_scalar(mat_elem_t *y, mat_elem_t alpha,
const mat_elem_t *x, size_t n) {
for (size_t i = 0; i < n; i++) {
y[i] += alpha * x[i];
}
}
#ifdef MAT_HAS_ARM_NEON
MAT_INTERNAL_STATIC void mat__axpy_neon(mat_elem_t *y, mat_elem_t alpha,
const mat_elem_t *x, size_t n) {
MAT_NEON_TYPE valpha = MAT_NEON_DUP(alpha);
size_t i = 0;
for (; i + MAT_NEON_WIDTH * 4 <= n; i += MAT_NEON_WIDTH * 4) {
MAT_NEON_TYPE vy0 = MAT_NEON_LOAD(&y[i]);
MAT_NEON_TYPE vy1 = MAT_NEON_LOAD(&y[i + MAT_NEON_WIDTH]);
MAT_NEON_TYPE vy2 = MAT_NEON_LOAD(&y[i + MAT_NEON_WIDTH * 2]);
MAT_NEON_TYPE vy3 = MAT_NEON_LOAD(&y[i + MAT_NEON_WIDTH * 3]);
MAT_NEON_TYPE vx0 = MAT_NEON_LOAD(&x[i]);
MAT_NEON_TYPE vx1 = MAT_NEON_LOAD(&x[i + MAT_NEON_WIDTH]);
MAT_NEON_TYPE vx2 = MAT_NEON_LOAD(&x[i + MAT_NEON_WIDTH * 2]);
MAT_NEON_TYPE vx3 = MAT_NEON_LOAD(&x[i + MAT_NEON_WIDTH * 3]);
vy0 = MAT_NEON_FMA(vy0, vx0, valpha);
vy1 = MAT_NEON_FMA(vy1, vx1, valpha);
vy2 = MAT_NEON_FMA(vy2, vx2, valpha);
vy3 = MAT_NEON_FMA(vy3, vx3, valpha);
MAT_NEON_STORE(&y[i], vy0);
MAT_NEON_STORE(&y[i + MAT_NEON_WIDTH], vy1);
MAT_NEON_STORE(&y[i + MAT_NEON_WIDTH * 2], vy2);
MAT_NEON_STORE(&y[i + MAT_NEON_WIDTH * 3], vy3);
}
for (; i + MAT_NEON_WIDTH <= n; i += MAT_NEON_WIDTH) {
MAT_NEON_TYPE vy = MAT_NEON_LOAD(&y[i]);
MAT_NEON_TYPE vx = MAT_NEON_LOAD(&x[i]);
vy = MAT_NEON_FMA(vy, vx, valpha);
MAT_NEON_STORE(&y[i], vy);
}
for (; i < n; i++) {
y[i] += alpha * x[i];
}
}
#endif
MAT_INTERNAL_STATIC void mat__axpy_kernel(mat_elem_t *y, mat_elem_t alpha,
const mat_elem_t *x, size_t n) {
MAT_DISPATCH(axpy)(y, alpha, x, n);
}
// SCAL: y[0:n] *= alpha
MAT_INTERNAL_STATIC void mat__scal_scalar(mat_elem_t *y, mat_elem_t alpha, size_t n) {
for (size_t i = 0; i < n; i++) {
y[i] *= alpha;
}
}
#ifdef MAT_HAS_ARM_NEON
MAT_INTERNAL_STATIC void mat__scal_neon(mat_elem_t *y, mat_elem_t alpha, size_t n) {
MAT_NEON_TYPE valpha = MAT_NEON_DUP(alpha);
size_t i = 0;
for (; i + MAT_NEON_WIDTH * 4 <= n; i += MAT_NEON_WIDTH * 4) {
MAT_NEON_TYPE vy0 = MAT_NEON_LOAD(&y[i]);
MAT_NEON_TYPE vy1 = MAT_NEON_LOAD(&y[i + MAT_NEON_WIDTH]);
MAT_NEON_TYPE vy2 = MAT_NEON_LOAD(&y[i + MAT_NEON_WIDTH * 2]);
MAT_NEON_TYPE vy3 = MAT_NEON_LOAD(&y[i + MAT_NEON_WIDTH * 3]);
MAT_NEON_STORE(&y[i], MAT_NEON_MUL(vy0, valpha));
MAT_NEON_STORE(&y[i + MAT_NEON_WIDTH], MAT_NEON_MUL(vy1, valpha));
MAT_NEON_STORE(&y[i + MAT_NEON_WIDTH * 2], MAT_NEON_MUL(vy2, valpha));
MAT_NEON_STORE(&y[i + MAT_NEON_WIDTH * 3], MAT_NEON_MUL(vy3, valpha));
}
for (; i + MAT_NEON_WIDTH <= n; i += MAT_NEON_WIDTH) {
MAT_NEON_TYPE vy = MAT_NEON_LOAD(&y[i]);
MAT_NEON_STORE(&y[i], MAT_NEON_MUL(vy, valpha));
}
for (; i < n; i++) {
y[i] *= alpha;
}
}
#endif
MAT_INTERNAL_STATIC void mat__scal_kernel(mat_elem_t *y, mat_elem_t alpha, size_t n) {
if (alpha == 1) return;
if (alpha == 0) {