// Banded symmetric positive-definite linear algebra for the static condensation: a reverse
// Cuthill-McKee ordering of the node graph to make the bandwidth small, a banded Cholesky
// factorisation, and forward/back substitution for many right-hand sides at once.

/** Reverse Cuthill-McKee ordering of an undirected graph given as CSR adjacency (no self loops
 * required). Returns perm with perm[newIndex] = oldIndex. */
export function reverseCuthillMcKee(n: number, ptr: Uint32Array, idx: Uint32Array): Uint32Array {
  const deg = new Uint32Array(n);
  for (let i = 0; i < n; i++) deg[i] = ptr[i + 1] - ptr[i];
  const visited = new Uint8Array(n);
  const order: number[] = [];
  const byDegree = Array.from({ length: n }, (_, i) => i).sort((a, b) => deg[a] - deg[b]);
  for (const start of byDegree) {
    if (visited[start]) continue;
    // pseudo-peripheral start: walk to the last node of a BFS level structure twice
    let s = start;
    for (let rep = 0; rep < 2; rep++) {
      const seen = new Uint8Array(n);
      let frontier = [s];
      seen[s] = 1;
      let last = s;
      while (frontier.length) {
        const next: number[] = [];
        for (const u of frontier) {
          for (let k = ptr[u]; k < ptr[u + 1]; k++) {
            const v = idx[k];
            if (!seen[v] && !visited[v]) { seen[v] = 1; next.push(v); }
          }
        }
        if (next.length) { next.sort((a, b) => deg[a] - deg[b]); last = next[0]; }
        frontier = next;
      }
      s = last;
    }
    const queue: number[] = [s];
    visited[s] = 1;
    let head = 0;
    while (head < queue.length) {
      const u = queue[head++];
      order.push(u);
      const nb: number[] = [];
      for (let k = ptr[u]; k < ptr[u + 1]; k++) {
        const v = idx[k];
        if (!visited[v]) { visited[v] = 1; nb.push(v); }
      }
      nb.sort((a, b) => deg[a] - deg[b]);
      for (const v of nb) queue.push(v);
    }
  }
  order.reverse();
  return Uint32Array.from(order);
}

/** Symmetric band matrix, lower part stored row-wise: entry (i, j) with i - b <= j <= i lives at
 * data[i * (b + 1) + (j - i + b)]. */
export interface BandMatrix {
  n: number;
  b: number;
  data: Float64Array;
}

export function bandMatrix(n: number, b: number): BandMatrix {
  return { n, b, data: new Float64Array(n * (b + 1)) };
}

export function bandAdd(A: BandMatrix, i: number, j: number, v: number): void {
  if (j > i) { const t = i; i = j; j = t; }
  const off = j - i + A.b;
  if (off < 0) throw new Error(`entry (${i}, ${j}) outside the band ${A.b}`);
  A.data[i * (A.b + 1) + off] += v;
}

/** In-place Cholesky A = L L^T of a banded SPD matrix; returns false if a pivot is not positive. */
export function bandedCholesky(A: BandMatrix, onProgress?: (frac: number) => void): boolean {
  const { n, b, data } = A;
  const w = b + 1;
  for (let i = 0; i < n; i++) {
    const ri = i * w;
    const j0 = Math.max(0, i - b);
    for (let j = j0; j <= i; j++) {
      const rj = j * w;
      let s = data[ri + (j - i + b)];
      const k0 = Math.max(j0, j - b);
      // s -= sum_k L[i,k] L[j,k] for k in [k0, j)
      let pi = ri + (k0 - i + b), pj = rj + (k0 - j + b);
      for (let k = k0; k < j; k++) s -= data[pi++] * data[pj++];
      if (j === i) {
        if (!(s > 0)) return false;
        data[ri + b] = Math.sqrt(s);
      } else {
        data[ri + (j - i + b)] = s / data[rj + b];
      }
    }
    if (onProgress && (i & 255) === 0) onProgress(i / n);
  }
  return true;
}

/** Solve L L^T X = R for X, with R given row-major (n x m) and overwritten by X. */
export function bandedSolve(L: BandMatrix, R: Float64Array, m: number, onProgress?: (frac: number) => void): void {
  const { n, b, data } = L;
  const w = b + 1;
  // forward: L y = r
  for (let i = 0; i < n; i++) {
    const ri = i * w, rowI = i * m;
    const k0 = Math.max(0, i - b);
    for (let k = k0; k < i; k++) {
      const l = data[ri + (k - i + b)];
      if (l === 0) continue;
      const rowK = k * m;
      for (let c = 0; c < m; c++) R[rowI + c] -= l * R[rowK + c];
    }
    const d = 1 / data[ri + b];
    for (let c = 0; c < m; c++) R[rowI + c] *= d;
    if (onProgress && (i & 255) === 0) onProgress(0.5 * i / n);
  }
  // backward: L^T x = y  (column i of L^T is row i of L)
  for (let i = n - 1; i >= 0; i--) {
    const ri = i * w, rowI = i * m;
    const d = 1 / data[ri + b];
    for (let c = 0; c < m; c++) R[rowI + c] *= d;
    const k0 = Math.max(0, i - b);
    for (let k = k0; k < i; k++) {
      const l = data[ri + (k - i + b)];
      if (l === 0) continue;
      const rowK = k * m;
      for (let c = 0; c < m; c++) R[rowK + c] -= l * R[rowI + c];
    }
    if (onProgress && (i & 255) === 0) onProgress(0.5 + 0.5 * (n - i) / n);
  }
}