// bend-blas, the JS twins: plain loops, correct and slow, for the // interpreter and the tests. An Array is a JS array of numbers. // A failure is io_fail(code): the runtime's message for a code is not // this library's, so the code is the contract on this lane. function bend_fresh(n) { let slots = 1; while (slots < n) slots *= 2; return new Array(slots).fill(0); } function bend_copy(src, n) { const c = bend_fresh(n); for (let i = 0; i < n; i++) c[i] = src[i]; return c; } function bend_short(a, n) { return n > 0x7fffffff || a.length < n; } function bend_gemm(ta, tb, m, n, k, alpha, a, b, beta, c) { const out = bend_copy(c, m * n); for (let i = 0; i < m; i++) for (let j = 0; j < n; j++) { let acc = 0; for (let l = 0; l < k; l++) acc += (ta ? a[l * m + i] : a[i * k + l]) * (tb ? b[j * k + l] : b[l * n + j]); out[i * n + j] = alpha * acc + beta * out[i * n + j]; } return out; } function lapack_sgesv(n, nrhs, a, b) { if (bend_short(a, n * n)) return io_fail(2); if (bend_short(b, n * nrhs)) return io_fail(2); const A = bend_copy(a, n * n), X = bend_copy(b, n * nrhs); for (let k = 0; k < n; k++) { let p = k; for (let i = k + 1; i < n; i++) if (Math.abs(A[i * n + k]) > Math.abs(A[p * n + k])) p = i; if (A[p * n + k] === 0) return io_fail(k + 1); if (p !== k) { for (let j = 0; j < n; j++) { const t = A[k * n + j]; A[k * n + j] = A[p * n + j]; A[p * n + j] = t; } for (let j = 0; j < nrhs; j++) { const t = X[k * nrhs + j]; X[k * nrhs + j] = X[p * nrhs + j]; X[p * nrhs + j] = t; } } for (let i = k + 1; i < n; i++) { const m = A[i * n + k] / A[k * n + k]; for (let j = k; j < n; j++) A[i * n + j] -= m * A[k * n + j]; for (let j = 0; j < nrhs; j++) X[i * nrhs + j] -= m * X[k * nrhs + j]; } } for (let i = n - 1; i >= 0; i--) for (let j = 0; j < nrhs; j++) { let s = X[i * nrhs + j]; for (let l = i + 1; l < n; l++) s -= A[i * n + l] * X[l * nrhs + j]; X[i * nrhs + j] = s / A[i * n + i]; } return io_done(X); } io_eff(CID(Lapack.sgesv), lapack_sgesv);