Source

`tests/purrc/linalg_sys_test.cpp`

2// Original work; see ACKNOWLEDGMENTS.md for the open-source ideas we build upon.

3// System-level end-to-end test for the cheatah `linalg` stdlib module.

4//

5// Unlike the per-function compile-run tests (linalg_cr_test.cpp), this is ONE

6// cohesive linear-algebra program that exercises EVERY purr-callable routine in

7// stdlib/linalg/routines.hpp in a single run: the 7 matrix/vector products

8// (dot, vdot, inner, outer, matmul, matrix_power, kron), the 3 decompositions

9// (cholesky, qr, svd), the 4 eigen routines (eig, eigvals, eigh, eigvalsh), the

10// 6 "numbers" routines (norm, cond, det, matrix_rank, slogdet, trace), and the

11// 4 solve/invert routines (solve, lstsq, inv, pinv) — 24 functions total.

12//

13// It builds a few fixed operands (two vectors, an SPD diagonal matrix, a general

14// 2x2, and the 2x2 identity) and runs a sequence of operations over them. All

15// printed values are integer-valued or exactly representable so io.print's

16// formatting is deterministic. Struct-returning routines print a derived single

17// field (qr().r, svd().s, eig().values, eigh().values, slogdet().sign).

18#include "e2e_harness.hpp"

20TEST(StdlibE2E, Linalg) {

21 e2e::expect_e2e("linalg_sys", R"PURR(import io

22import ndarray

23import linalg

25# Fixed building blocks: a couple of small vectors and matrices.

26let u = ndarray.array([1.0, 2.0, 3.0])

27let v = ndarray.array([4.0, 5.0, 6.0])

29# A symmetric positive-definite diagonal matrix and a general 2x2.

30let d = ndarray.reshape(ndarray.array([4.0, 0.0, 0.0, 9.0]), [2, 2])

31let g = ndarray.reshape(ndarray.array([4.0, 3.0, 6.0, 3.0]), [2, 2])

32let id2 = ndarray.reshape(ndarray.array([1.0, 0.0, 0.0, 1.0]), [2, 2])

34# ---- products ----

35io.print(linalg.dot(u, v))

36io.print(linalg.vdot(u, v))

37io.print(linalg.inner(u, v))

38io.print(ndarray.to_string(linalg.outer(ndarray.array([1.0, 2.0]), ndarray.array([3.0, 4.0]))))

39io.print(ndarray.to_string(linalg.matmul(d, id2)))

40io.print(ndarray.to_string(linalg.matrix_power(d, 2)))

41io.print(ndarray.to_string(linalg.kron(id2, ndarray.reshape(ndarray.array([1.0, 2.0, 3.0, 4.0]), [2, 2]))))

43# ---- decompositions (print derived scalars / single fields) ----

44io.print(ndarray.to_string(linalg.cholesky(d)))

45io.print(ndarray.to_string(linalg.qr(d).r))

46io.print(ndarray.to_string(linalg.svd(d).s))

48# ---- eigenvalues ----

49io.print(ndarray.to_string(linalg.eig(d).values))

50io.print(ndarray.to_string(linalg.eigvals(d)))

51io.print(ndarray.to_string(linalg.eigh(d).values))

52io.print(ndarray.to_string(linalg.eigvalsh(d)))

54# ---- norms and numbers ----

55io.print(linalg.norm(ndarray.array([3.0, 4.0])))

56io.print(linalg.cond(d))

57io.print(linalg.det(g))

58io.print(linalg.matrix_rank(d))

59io.print(linalg.slogdet(d).sign)

60io.print(linalg.trace(g))

62# ---- solving and inverting ----

63io.print(ndarray.to_string(linalg.solve(g, ndarray.array([10.0, 12.0]))))

64io.print(ndarray.to_string(linalg.lstsq(id2, ndarray.reshape(ndarray.array([5.0, 7.0]), [2, 1]))))

65io.print(ndarray.to_string(linalg.inv(d)))

66io.print(ndarray.to_string(linalg.pinv(d)))

67)PURR",

68 "32\n"

69 "32\n"

70 "32\n"

71 "[[3, 4], [6, 8]]\n"

72 "[[4, 0], [0, 9]]\n"

73 "[[16, 0], [0, 81]]\n"

74 "[[1, 2, 0, 0], [3, 4, 0, 0], [0, 0, 1, 2], [0, 0, 3, 4]]\n"

75 "[[2, 0], [0, 3]]\n"

76 "[[-4, 0], [0, -9]]\n"

77 "[9, 4]\n" // svd().s

78 "[9+0j, 4+0j]\n" // eig().values — general -> complex spectrum

79 "[9+0j, 4+0j]\n" // eigvals() — general -> complex spectrum

80 "[9, 4]\n" // eigh().values — Hermitian -> real

81 "[9, 4]\n" // eigvalsh() — Hermitian -> real

82 "5\n"

83 "2.25\n"

84 "-6\n"

85 "2\n"

86 "1\n"

87 "7\n"

88 "[1, 2]\n"

89 "[[5], [7]]\n"

90 "[[0.25, 0], [0, 0.111111]]\n"

91 "[[0.25, 0], [0, 0.111111]]\n");

92}

94// Complex linear algebra exercised together end-to-end: build complex vectors and a

95// matrix, take bilinear (dot) and Hermitian (vdot) inner products, the conjugate

96// transpose, and a complex matmul. M·Mᴴ is Hermitian; vdot(a,a) is the real ‖a‖².

97TEST(StdlibE2E, LinalgComplex) {

98 e2e::expect_e2e("linalg_complex_sys", R"PURR(import io

99import ndarray

100import linalg

101

102let a = ndarray.complex(ndarray.array([1.0, 3.0]), ndarray.array([2.0, -1.0]))

103let b = ndarray.complex(ndarray.array([0.0, 2.0]), ndarray.array([1.0, 0.0]))

104

let M = ndarray.reshape(ndarray.complex(ndarray.array([1.0, 2.0, 0.0, 3.0]), ndarray.array([1.0, 0.0, 0.0, -1.0])), [2, 2])

105

let H = ndarray.reshape(ndarray.complex(ndarray.array([2.0, 1.0, 1.0, 3.0]), ndarray.array([0.0, 1.0, -1.0, 0.0])), [2, 2])

106

107io.print(linalg.dot(a, b))

108io.print(linalg.vdot(a, b))

109io.print(linalg.vdot(a, a))

110io.print(ndarray.to_string(linalg.conj_transpose(M)))

111io.print(ndarray.to_string(linalg.matmul(M, linalg.conj_transpose(M))))

112io.print(ndarray.to_string(linalg.eigvalsh(H)))

113io.print(ndarray.to_string(linalg.eigh(H).values))

114)PURR",

115 "4-1j\n"

116 "8+3j\n"

117 "15+0j\n"

118 "[[1-1j, 0+0j], [2+0j, 3+1j]]\n"

119 "[[6+0j, 6+2j], [6-2j, 10+0j]]\n"

120 "[4, 1]\n"

121 "[4, 1]\n");

122}