Files
tensor/integrate/tridiag.go
petrbalvin af4ee19703
Release / gates (push) Successful in 4m38s
Test / test (push) Successful in 5m16s
Release / release (push) Successful in 35s
feat: initial release
Assisted-by: GLM 5.3 Flash
2026-09-03 10:00:00 +02:00

40 lines
1.7 KiB
Go

// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
// SPDX-License-Identifier: MIT
package integrate
// The tridiagonal scratch the PDE step sweeps reuse. The Crank-
// Nicolson and Peaceman-Rachford schemes solve one tridiagonal system
// per line per step against matrices that are constants of the scheme,
// so the elimination scratch belongs to the solve, not to the line:
// the sweeps below allocate it once and reuse it across every line and
// every step, where the library's array-level SolveTridiagonal
// allocated its working vectors per call. The elimination itself is
// internal/base's TriSolve, shared with the public solver, so the
// arithmetic exists once.
// triScratch holds one tridiagonal solve's reusable buffers: the
// right-hand side, the eliminated superdiagonal and right side of the
// Thomas algorithm, and the solution buffer for callers whose output
// is not written straight into their state. aux is a second right-side
// lane for the sweeps that build two neighbouring systems in one pass.
// Every buffer is fully overwritten before the kernel reads it, except
// cp, whose prefix is written and read in the same sweep order the
// fresh buffers saw.
type triScratch struct {
rhs, cp, dp, dst, aux []float64
}
// triSized sizes s for a system of sys unknowns whose right side is a
// line of lineLen entries, reusing storage that is already large
// enough. The solution buffer dst is sized for callers that scatter
// the result elsewhere; a caller writing straight into its state leaves
// it unread.
func triSized(s *triScratch, sys, lineLen int) {
s.rhs = sizedBuf(s.rhs, lineLen)
s.aux = sizedBuf(s.aux, lineLen)
s.cp = sizedBuf(s.cp, sys)
s.dp = sizedBuf(s.dp, sys)
s.dst = sizedBuf(s.dst, sys)
}