Source

`tests/purrc/linalg_cr_test.cpp`

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

3// Compile-run unit tests for the `linalg` module: one test per purr-callable

4// function. Each writes a tiny .purr that calls a single linalg routine on a

5// small fixed matrix/vector, compiles it with purrc, runs it under the cheatah

6// runtime, and asserts the exact stdout. Results are chosen to be integer-valued

7// (or exactly representable) so the io.print formatting is deterministic.

8// Complements the in-process unit tests (stdlib/tests/linalg_routines_test.cpp),

9// the per-module system test (StdlibE2E.Linalg), and the least-squares system

10// test (linalg_lsq_e2e_test.cpp).

11#include "e2e_harness.hpp"

13// ---- Matrix and vector products ----

15TEST(LinalgCompileRun, Dot) {

16 e2e::expect_e2e("linalg_dot", R"PURR(import io

17import ndarray

18import linalg

19io.print(linalg.dot(ndarray.array([1.0, 2.0, 3.0]), ndarray.array([4.0, 5.0, 6.0])))

20)PURR", "32\n");

21}

23TEST(LinalgCompileRun, Vdot) {

24 e2e::expect_e2e("linalg_vdot", R"PURR(import io

25import ndarray

26import linalg

27io.print(linalg.vdot(ndarray.array([1.0, 2.0, 3.0]), ndarray.array([4.0, 5.0, 6.0])))

28)PURR", "32\n");

29}

31TEST(LinalgCompileRun, Inner) {

32 e2e::expect_e2e("linalg_inner", R"PURR(import io

33import ndarray

34import linalg

35io.print(linalg.inner(ndarray.array([1.0, 2.0, 3.0]), ndarray.array([4.0, 5.0, 6.0])))

36)PURR", "32\n");

37}

39TEST(LinalgCompileRun, Outer) {

40 e2e::expect_e2e("linalg_outer", R"PURR(import io

41import ndarray

42import linalg

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

44)PURR", "[[3, 4], [6, 8]]\n");

45}

47TEST(LinalgCompileRun, Matmul) {

48 e2e::expect_e2e("linalg_matmul", R"PURR(import io

49import ndarray

50import linalg

51let a = ndarray.reshape(ndarray.array([1.0, 2.0, 3.0, 4.0]), [2, 2])

52let b = ndarray.reshape(ndarray.array([5.0, 6.0, 7.0, 8.0]), [2, 2])

53io.print(ndarray.to_string(linalg.matmul(a, b)))

54)PURR", "[[19, 22], [43, 50]]\n");

55}

57TEST(LinalgCompileRun, MatrixPower) {

58 e2e::expect_e2e("linalg_matrix_power", R"PURR(import io

59import ndarray

60import linalg

61let a = ndarray.reshape(ndarray.array([2.0, 0.0, 0.0, 3.0]), [2, 2])

62io.print(ndarray.to_string(linalg.matrix_power(a, 3)))

63)PURR", "[[8, 0], [0, 27]]\n");

64}

66TEST(LinalgCompileRun, Kron) {

67 e2e::expect_e2e("linalg_kron", R"PURR(import io

68import ndarray

69import linalg

70let a = ndarray.reshape(ndarray.array([1.0, 0.0, 0.0, 1.0]), [2, 2])

71let b = ndarray.reshape(ndarray.array([1.0, 2.0, 3.0, 4.0]), [2, 2])

72io.print(ndarray.to_string(linalg.kron(a, b)))

73)PURR", "[[1, 2, 0, 0], [3, 4, 0, 0], [0, 0, 1, 2], [0, 0, 3, 4]]\n");

74}

76// ---- Decompositions ----

78TEST(LinalgCompileRun, Cholesky) {

79 e2e::expect_e2e("linalg_cholesky", R"PURR(import io

80import ndarray

81import linalg

82let a = ndarray.reshape(ndarray.array([4.0, 0.0, 0.0, 9.0]), [2, 2])

83io.print(ndarray.to_string(linalg.cholesky(a)))

84)PURR", "[[2, 0], [0, 3]]\n");

85}

87TEST(LinalgCompileRun, Qr) {

88 // Diagonal input -> Householder gives diagonal R with the reflector's sign

89 // convention (negated diagonal), which is exact and deterministic.

90 e2e::expect_e2e("linalg_qr", R"PURR(import io

91import ndarray

92import linalg

93let a = ndarray.reshape(ndarray.array([6.0, 0.0, 0.0, 5.0]), [2, 2])

94let f = linalg.qr(a)

95io.print(ndarray.to_string(f.r))

96)PURR", "[[-6, 0], [0, -5]]\n");

97}

99TEST(LinalgCompileRun, Svd) {

100 e2e::expect_e2e("linalg_svd", R"PURR(import io

101import ndarray

102import linalg

103let a = ndarray.reshape(ndarray.array([4.0, 0.0, 0.0, 9.0]), [2, 2])

104let f = linalg.svd(a)

105io.print(ndarray.to_string(f.s))

106)PURR", "[9, 4]\n");

107}

108

109TEST(LinalgCompileRun, Svdvals) {

110 e2e::expect_e2e("linalg_svdvals", R"PURR(import io

111import ndarray

112import linalg

113let a = ndarray.reshape(ndarray.array([4.0, 0.0, 0.0, 9.0]), [2, 2])

114io.print(ndarray.to_string(linalg.svdvals(a)))

115)PURR", "[9, 4]\n");

116}

117

118// ---- Matrix eigenvalues ----

119

120TEST(LinalgCompileRun, Eig) {

121 e2e::expect_e2e("linalg_eig", R"PURR(import io

122import ndarray

123import linalg

124let a = ndarray.reshape(ndarray.array([2.0, 0.0, 0.0, 3.0]), [2, 2])

125let e = linalg.eig(a)

126io.print(ndarray.to_string(e.values))

127)PURR", "[3+0j, 2+0j]\n"); // general eig -> complex spectrum (real parts here)

128}

129

130TEST(LinalgCompileRun, Eigvals) {

131 e2e::expect_e2e("linalg_eigvals", R"PURR(import io

132import ndarray

133import linalg

134let a = ndarray.reshape(ndarray.array([2.0, 0.0, 0.0, 3.0]), [2, 2])

135io.print(ndarray.to_string(linalg.eigvals(a)))

136)PURR", "[3+0j, 2+0j]\n");

137}

138

139TEST(LinalgCompileRun, EigvalsComplex) {

140 // A real rotation matrix [[0,-1],[1,0]] has eigenvalues ±i — printed Python-style.

141 e2e::expect_e2e("linalg_eigvals_complex", R"PURR(import io

142import ndarray

143import linalg

144let a = ndarray.reshape(ndarray.array([0.0, -1.0, 1.0, 0.0]), [2, 2])

145io.print(ndarray.to_string(linalg.eigvals(a)))

146)PURR", "[0+1j, 0-1j]\n");

147}

148

149TEST(LinalgCompileRun, Eigh) {

150 e2e::expect_e2e("linalg_eigh", R"PURR(import io

151import ndarray

152import linalg

153let a = ndarray.reshape(ndarray.array([2.0, 0.0, 0.0, 5.0]), [2, 2])

154let e = linalg.eigh(a)

155io.print(ndarray.to_string(e.values))

156)PURR", "[5, 2]\n");

157}

158

159TEST(LinalgCompileRun, Eigvalsh) {

160 e2e::expect_e2e("linalg_eigvalsh", R"PURR(import io

161import ndarray

162import linalg

163let a = ndarray.reshape(ndarray.array([2.0, 0.0, 0.0, 5.0]), [2, 2])

164io.print(ndarray.to_string(linalg.eigvalsh(a)))

165)PURR", "[5, 2]\n");

166}

167

168// ---- Norms and other numbers ----

169

170TEST(LinalgCompileRun, Norm) {

171 e2e::expect_e2e("linalg_norm", R"PURR(import io

172import ndarray

173import linalg

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

175)PURR", "5\n");

176}

177

178TEST(LinalgCompileRun, Cond) {

179 e2e::expect_e2e("linalg_cond", R"PURR(import io

180import ndarray

181import linalg

182let a = ndarray.reshape(ndarray.array([2.0, 0.0, 0.0, 8.0]), [2, 2])

183io.print(linalg.cond(a))

184)PURR", "4\n");

185}

186

187TEST(LinalgCompileRun, Det) {

188 e2e::expect_e2e("linalg_det", R"PURR(import io

189import ndarray

190import linalg

191let a = ndarray.reshape(ndarray.array([1.0, 2.0, 3.0, 4.0]), [2, 2])

192io.print(linalg.det(a))

193)PURR", "-2\n");

194}

195

196TEST(LinalgCompileRun, MatrixRank) {

197 // Rank-1 matrix (second row = 2x first) -> numerical rank 1.

198 e2e::expect_e2e("linalg_matrix_rank", R"PURR(import io

199import ndarray

200import linalg

201let a = ndarray.reshape(ndarray.array([1.0, 2.0, 2.0, 4.0]), [2, 2])

202io.print(linalg.matrix_rank(a))

203)PURR", "1\n");

204}

205

206TEST(LinalgCompileRun, Slogdet) {

207 // det = 6 > 0, so the sign is exactly +1 (avoids the float logabsdet).

208 e2e::expect_e2e("linalg_slogdet", R"PURR(import io

209import ndarray

210import linalg

211let a = ndarray.reshape(ndarray.array([2.0, 0.0, 0.0, 3.0]), [2, 2])

212let r = linalg.slogdet(a)

213io.print(r.sign)

214)PURR", "1\n");

215}

216

217TEST(LinalgCompileRun, Trace) {

218 e2e::expect_e2e("linalg_trace", R"PURR(import io

219import ndarray

220import linalg

221let a = ndarray.reshape(ndarray.array([1.0, 2.0, 3.0, 4.0]), [2, 2])

222io.print(linalg.trace(a))

223)PURR", "5\n");

224}

225

226// ---- Solving equations and inverting matrices ----

227

228TEST(LinalgCompileRun, Solve) {

229 // 4x+3y=10, 6x+3y=12 -> x=1, y=2.

230 e2e::expect_e2e("linalg_solve", R"PURR(import io

231import ndarray

232import linalg

233let a = ndarray.reshape(ndarray.array([4.0, 3.0, 6.0, 3.0]), [2, 2])

234io.print(ndarray.to_string(linalg.solve(a, ndarray.array([10.0, 12.0]))))

235)PURR", "[1, 2]\n");

236}

237

238TEST(LinalgCompileRun, Lstsq) {

239 // Identity system with a column-vector RHS -> x = b exactly.

240 e2e::expect_e2e("linalg_lstsq", R"PURR(import io

241import ndarray

242import linalg

243let a = ndarray.reshape(ndarray.array([1.0, 0.0, 0.0, 1.0]), [2, 2])

244let b = ndarray.reshape(ndarray.array([5.0, 7.0]), [2, 1])

245io.print(ndarray.to_string(linalg.lstsq(a, b)))

246)PURR", "[[5], [7]]\n");

247}

248

249TEST(LinalgCompileRun, Inv) {

250 e2e::expect_e2e("linalg_inv", R"PURR(import io

251import ndarray

252import linalg

253let a = ndarray.reshape(ndarray.array([2.0, 0.0, 0.0, 4.0]), [2, 2])

254io.print(ndarray.to_string(linalg.inv(a)))

255)PURR", "[[0.5, 0], [0, 0.25]]\n");

256}

257

258TEST(LinalgCompileRun, Pinv) {

259 e2e::expect_e2e("linalg_pinv", R"PURR(import io

260import ndarray

261import linalg

262let a = ndarray.reshape(ndarray.array([2.0, 0.0, 0.0, 4.0]), [2, 2])

263io.print(ndarray.to_string(linalg.pinv(a)))

264)PURR", "[[0.5, 0], [0, 0.25]]\n");

265}

266

267// ---- complex products (complex inner-product spaces) ----

268

269TEST(LinalgCompileRun, ComplexDot) {

270 e2e::expect_e2e("linalg_complex_dot", R"PURR(import io

271import ndarray

272import linalg

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

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

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

276)PURR", "4-1j\n");

277}

278

279TEST(LinalgCompileRun, ComplexVdot) {

280 e2e::expect_e2e("linalg_complex_vdot", R"PURR(import io

281import ndarray

282import linalg

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

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

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

286)PURR", "8+3j\n");

287}

288

289TEST(LinalgCompileRun, ComplexMatmul) {

290 e2e::expect_e2e("linalg_complex_matmul", R"PURR(import io

291import ndarray

292import linalg

293

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

294

let I = ndarray.reshape(ndarray.complex(ndarray.array([1.0, 0.0, 0.0, 1.0]), ndarray.array([0.0, 0.0, 0.0, 0.0])), [2, 2])

295io.print(ndarray.to_string(linalg.matmul(M, I)))

296)PURR", "[[1+1j, 0+0j], [0+0j, 1+1j]]\n");

297}

298

299TEST(LinalgCompileRun, ConjTranspose) {

300 e2e::expect_e2e("linalg_conj_transpose", R"PURR(import io

301import ndarray

302import linalg

303

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])

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

305)PURR", "[[1-1j, 0+0j], [2+0j, 3+1j]]\n");

306}

307

308TEST(LinalgCompileRun, EighComplex) {

309 // Hermitian [[2, 1+i],[1-i, 3]] -> real eigenvalues 4, 1.

310 e2e::expect_e2e("linalg_eigh_complex", R"PURR(import io

311import ndarray

312import linalg

313

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])

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

315)PURR", "[4, 1]\n");

316}

317

318TEST(LinalgCompileRun, EigvalshComplex) {

319 e2e::expect_e2e("linalg_eigvalsh_complex", R"PURR(import io

320import ndarray

321import linalg

322

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])

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

324)PURR", "[4, 1]\n");

325}