CraftScript iterative solvers and sparse eigenvalues

/** PCG on a sparse SPD system. */
export function pcgSolve(S, b, { tol = 1e-8, maxit = 0, prec = 'auto', x0 = null, dofsPerNode = 1, coords = null, check = () => {} } = {}) {
  if (S.n !== S.m) throw err(t('script.linalg.pcg_square', { size: `${S.n}×${S.m}` }));
  if (b.length !== S.n) throw err(t('script.linalg.pcg_rhs', { b: b.length, n: S.n }));
  const t0 = Date.now(), M = typeof prec === 'object' ? prec : makePrec(S, prec, { dofsPerNode, coords }), t1 = Date.now();
  const r = kPcg(S, b, { x0, tol, maxIter: maxit, prec: M, cancelled: () => { check(); return false; } });
  return { x: r.x, iterations: r.iterations, residual: r.residual, converged: r.converged, prec: M.kind, setupMs: t1 - t0, ms: Date.now() - t0 };
}

// ─────────────────────────────────────────────────────────────── iterative solvers and preconditioners
function jacobiPrec(S) {
  const d = spDiag(S), inv = d.map((v) => (v ? 1 / v : 1));
  return { kind: 'jacobi', apply(r, z) { for (let i = 0; i < r.length; i++) z[i] = r[i] * inv[i]; return z; } };
}

/** ILU(0) of a square CSR with sorted columns → {apply}. */
export function ilu0(S) {
  const { n, rowPtr, col } = S, LU = Float64Array.from(S.val), diag = new Int32Array(n).fill(-1);
  for (let r = 0; r < n; r++) for (let k = rowPtr[r]; k < rowPtr[r + 1]; k++) if (col[k] === r) diag[r] = k;
  for (let r = 0; r < n; r++) if (diag[r] < 0) throw err(t('script.linalg.ilu0_diag'));
  const pos = new Int32Array(n).fill(-1);
  for (let i = 0; i < n; i++) {
    for (let k = rowPtr[i]; k < rowPtr[i + 1]; k++) pos[col[k]] = k;
    for (let k = rowPtr[i]; k < rowPtr[i + 1] && col[k] < i; k++) {
      const j = col[k], piv = LU[diag[j]];
      if (piv === 0) throw err(t('script.linalg.ilu0_pivot'));
      const f = (LU[k] /= piv);
      for (let q = diag[j] + 1; q < rowPtr[j + 1]; q++) { const p = pos[col[q]]; if (p >= 0) LU[p] -= f * LU[q]; }
    }
    for (let k = rowPtr[i]; k < rowPtr[i + 1]; k++) pos[col[k]] = -1;
  }
  return {
    kind: 'ilu0',
    apply(r, z) {
      for (let i = 0; i < n; i++) { let s = r[i]; for (let k = rowPtr[i]; k < diag[i]; k++) s -= LU[k] * z[col[k]]; z[i] = s; }
      for (let i = n - 1; i >= 0; i--) { let s = z[i]; for (let k = diag[i] + 1; k < rowPtr[i + 1]; k++) s -= LU[k] * z[col[k]]; z[i] = s / LU[diag[i]]; }
      return z;
    },
  };
}

/** Right-preconditioned BiCGSTAB (van der Vorst 1992). */
export function bicgstabSolve(S, b, { tol = 1e-8, maxit = 0, prec = 'ilu0', x0 = null, check = () => {} } = {}) {
  const n = S.n, t0 = Date.now();
  const M = prec === 'ilu0' ? ilu0(S) : prec === 'jacobi' ? jacobiPrec(S) : { kind: 'none', apply(r, z) { z.set(r); return z; } };
  const x = x0 ? Float64Array.from(x0) : new Float64Array(n), r = new Float64Array(n), rh = new Float64Array(n), p = new Float64Array(n), v = new Float64Array(n);
  const s = new Float64Array(n), t = new Float64Array(n), ph = new Float64Array(n), sh = new Float64Array(n);
  const dot = (a, c) => { let q = 0; for (let i = 0; i < n; i++) q += a[i] * c[i]; return q; };
  spMulVec(S, x, r); for (let i = 0; i < n; i++) r[i] = b[i] - r[i];
  rh.set(r);
  const bn = Math.sqrt(dot(b, b)) || 1;
  let rho = 1, alpha = 1, omega = 1, it = 0, res = Math.sqrt(dot(r, r)) / bn;
  const cap = maxit || Math.max(1000, 2 * n);
  for (; it < cap && res >= tol; it++) {
    if ((it & 31) === 0) check();
    const rho1 = dot(rh, r);
    if (rho1 === 0) break;
    const beta = (rho1 / rho) * (alpha / omega); rho = rho1;
    for (let i = 0; i < n; i++) p[i] = r[i] + beta * (p[i] - omega * v[i]);
    M.apply(p, ph); spMulVec(S, ph, v);
    alpha = rho / dot(rh, v);
    for (let i = 0; i < n; i++) s[i] = r[i] - alpha * v[i];
    if (Math.sqrt(dot(s, s)) / bn < tol) { for (let i = 0; i < n; i++) x[i] += alpha * ph[i]; res = Math.sqrt(dot(s, s)) / bn; it++; break; }
    M.apply(s, sh); spMulVec(S, sh, t);
    omega = dot(t, s) / dot(t, t);
    for (let i = 0; i < n; i++) { x[i] += alpha * ph[i] + omega * sh[i]; r[i] = s[i] - omega * t[i]; }
    res = Math.sqrt(dot(r, r)) / bn;
    if (omega === 0) break;
  }
  return { x, iterations: it, residual: res, converged: res < tol, prec: M.kind, ms: Date.now() - t0 };
}

/** Near-null space for AMG: constants per dof direction, plus rotations when node coordinates are given. */
function nearNull(n, dpn, coords) {
  const B = [];
  for (let d = 0; d < dpn; d++) { const v = new Float64Array(n); for (let i = d; i < n; i += dpn) v[i] = 1; B.push(v); }
  if (coords && (dpn === 2 || dpn === 3)) {
    const nn = n / dpn, X = coords;
    if (X.length < nn * dpn) throw err(t('script.linalg.pcg_coords', { dpn, nn }));
    let c = [0, 0, 0]; for (let i = 0; i < nn; i++) for (let a = 0; a < dpn; a++) c[a] += X[i * dpn + a] / nn;
    if (dpn === 2) { const v = new Float64Array(n); for (let i = 0; i < nn; i++) { v[2 * i] = -(X[2 * i + 1] - c[1]); v[2 * i + 1] = X[2 * i] - c[0]; } B.push(v); }
    else for (const [a, b] of [[0, 1], [1, 2], [2, 0]]) { const v = new Float64Array(n); for (let i = 0; i < nn; i++) { v[3 * i + a] = -(X[3 * i + b] - c[b]); v[3 * i + b] = X[3 * i + a] - c[a]; } B.push(v); }
  }
  return B;
}

// ─────────────────────────────────────────────────────────────── eigs (small symmetric generalized)
export function eigsSolve(K, M, k, { sigma = 0, tol = 1e-10, maxit = 300, check = () => {} } = {}) {
  const n = isSparse(K) ? K.n : size2(K)[0];
  if (!Number.isInteger(k) || k < 1 || k > n) throw err(t('script.linalg.eigs_k', { n }));
  if (M && (isSparse(M) ? M.n : size2(M)[0]) !== n) throw err(t('script.linalg.eigs_size'));
  if (n <= DENSE_EIG) {
    const Kd = isSparse(K) ? toDense(K) : asTensor(K), Md = M ? (isSparse(M) ? toDense(M) : asTensor(M)) : null;
    const { d, V } = Md ? symGenEigen(Kd, Md) : symEigen(Kd);
    return { d: d.slice(0, k), V: V.slice(0, n * k), iterations: 0, method: 'dense' };
  }
  const Ks = isSparse(K) ? K : fromDense(asTensor(K)), Ms = M ? (isSparse(M) ? M : fromDense(asTensor(M))) : null;
  const shifted = sigma ? spAdd(Ks, Ms || speyeS(n), 1, -sigma) : Ks;
  let F = envelopeCholesky(shifted, { check });
  const solveK = F ? (b, x) => F.solve(b, x) : (() => { const P = makePrec(shifted, 'auto'); return (b, x) => { const r = pcgSolve(shifted, b, { tol: 1e-12, prec: P, check }); x.set(r.x); return x; }; })();
  const Mmul = Ms ? (x, y) => spMulVec(Ms, x, y) : (x, y) => { y.set(x); return y; };
  const p = Math.min(n, Math.max(2 * k, k + 8));
  let X = new Float64Array(n * p);
  let s = 12345; for (let i = 0; i < X.length; i++) { s = (s * 16807) % 2147483647; X[i] = s / 2147483647 - 0.5; }
  let lam = new Float64Array(k).fill(Infinity), it = 0;
  const Y = new Float64Array(n * p), tw = new Float64Array(n), u = new Float64Array(n);
  for (; it < maxit; it++) {
    check();
    for (let j = 0; j < p; j++) { Mmul(X.subarray(j * n, (j + 1) * n), tw); solveK(tw, Y.subarray(j * n, (j + 1) * n)); }
    // Rayleigh–Ritz: Kr = Yᵀ(K−σM)Y, Mr = YᵀMY
    const Kr = new Float64Array(p * p), Mr = new Float64Array(p * p), KY = new Float64Array(n * p), MY = new Float64Array(n * p);
    for (let j = 0; j < p; j++) { spMulVec(shifted, Y.subarray(j * n, (j + 1) * n), u); KY.set(u, j * n); Mmul(Y.subarray(j * n, (j + 1) * n), u); MY.set(u, j * n); }
    for (let a = 0; a < p; a++) for (let b = a; b < p; b++) { let s1 = 0, s2 = 0; for (let i = 0; i < n; i++) { const ya = Y[a * n + i]; s1 += ya * KY[b * n + i]; s2 += ya * MY[b * n + i]; } Kr[a + b * p] = Kr[b + a * p] = s1; Mr[a + b * p] = Mr[b + a * p] = s2; }
    const { d, V: Q } = symGenEigen(mat(p, p, Kr), mat(p, p, Mr));
    X = new Float64Array(n * p);
    for (let j = 0; j < p; j++) for (let a = 0; a < p; a++) { const q = Q[a + j * p]; if (q) for (let i = 0; i < n; i++) X[j * n + i] += Y[a * n + i] * q; }
    let done = true;
    for (let j = 0; j < k; j++) { if (!(Math.abs(d[j] - lam[j]) <= tol * Math.max(1e-300, Math.abs(d[j])))) done = false; lam[j] = d[j]; }
    if (done) break;
  }
  return { d: Float64Array.from(lam, (v) => v + sigma), V: X.slice(0, n * k), iterations: it + 1, method: F ? 'subspace iteration (envelope Cholesky)' : 'subspace iteration (pcg)' };   // i18n-ignore: method names scripts read
}