// Symmetric banded LDL^T factorization of K - sigma M, its solve, and its
// inertia (the number of negative pivots, which by Sylvester's law equals
// the number of eigenvalues of the pencil (K, M) below sigma). The band
// layout is that of assemble.ts: row i holds K[i, i - d] at offset d.

export interface Factorization {
  n: number;
  bw: number;
  /** Unit lower band factor L, same layout as the input (the diagonal offset is unused). */
  L: Float64Array;
  /** Pivots D. */
  D: Float64Array;
  /** Number of negative pivots: eigenvalues of (K, M) below the shift. */
  negative: number;
  shift: number;
}

/** LDL^T of K - shift * diag(mass). Throws on an exactly zero pivot (a shift equal to an eigenvalue). */
export function factorize(band: Float64Array, mass: Float64Array, n: number, bw: number, shift: number): Factorization {
  const s = bw + 1;
  const L = new Float64Array(band);
  for (let i = 0; i < n; i++) L[i * s] -= shift * mass[i];
  const D = new Float64Array(n);
  // W[j][k] = D[k] L[j][k], kept alongside L so that the inner loop of the
  // factorization multiplies two loads rather than three.
  const W = new Float64Array(n * s);
  let negative = 0;
  for (let i = 0; i < n; i++) {
    const rowI = i * s;
    const kmin = Math.max(0, i - bw);
    for (let j = kmin; j < i; j++) {
      const rowJ = j * s;
      let acc = L[rowI + (i - j)];
      const k0 = Math.max(kmin, j - bw);
      // acc -= sum_k L[i][k] W[j][k]; L[i][k] sits at rowI + i - k, W[j][k] at rowJ + j - k.
      const offI = rowI + i, offJ = rowJ + j;
      for (let k = k0; k < j; k++) acc -= L[offI - k] * W[offJ - k];
      L[rowI + (i - j)] = acc / D[j];
    }
    let d = L[rowI];
    for (let k = kmin; k < i; k++) {
      const l = L[rowI + (i - k)];
      const w = l * D[k];
      W[rowI + (i - k)] = w;
      d -= l * w;
    }
    if (d === 0 || !isFinite(d)) throw new Error(`zero pivot at row ${i} for shift ${shift}`);
    D[i] = d;
    if (d < 0) negative++;
  }
  return { n, bw, L, D, negative, shift };
}

/** Solves (K - shift M) x = b in place (x may be b). */
export function solve(f: Factorization, b: Float64Array, x: Float64Array): void {
  const { n, bw, L, D } = f;
  const s = bw + 1;
  if (x !== b) x.set(b);
  // Forward: L y = b
  for (let i = 0; i < n; i++) {
    const rowI = i * s;
    const kmin = Math.max(0, i - bw);
    let acc = x[i];
    for (let k = kmin; k < i; k++) acc -= L[rowI + (i - k)] * x[k];
    x[i] = acc;
  }
  for (let i = 0; i < n; i++) x[i] /= D[i];
  // Backward: L^T x = y
  for (let i = n - 1; i >= 0; i--) {
    const xi = x[i];
    const rowI = i * s;
    const kmin = Math.max(0, i - bw);
    for (let k = kmin; k < i; k++) x[k] -= L[rowI + (i - k)] * xi;
  }
}

/** Number of eigenvalues of (K, M) below the shift. */
export function inertia(band: Float64Array, mass: Float64Array, n: number, bw: number, shift: number): number {
  return factorize(band, mass, n, bw, shift).negative;
}