General concepts » Boundary conditions » Essential constraint elimination

Algebraic formulation, sparse storage contracts, and implementation costs.

This reference expands the algebraic step introduced in the boundary conditions guide. It is intended for callers supplying their own sparse operators and for assembly developers.

Math::LinearSystem also accepts an arbitrary caller-supplied square CSC operator. Essential elimination supports nonsymmetric coefficients, asymmetric sparse patterns and missing constrained diagonals. Duplicate triplets must first be assembled into valid sorted, unique-row sparse storage, as with Eigen's setFromTriplets.

Projector formulation

Let $C$ denote the prescribed indices, $F$ the remaining indices, and $g_C$ their prescribed coefficients. Their interpretation is determined by the finite element DOF functionals. The free equations follow by substituting $x_C=g_C$ into the original system.

\[ A_{FF}x_F=b_F-A_{FC}g_C. \]

Define the diagonal constraint projector $Q$ , its complement $P=I-Q$ , and extend $g_C$ by zero on $F$ . Keeping the original system size gives

\[ \widetilde A=PAP+Q,\qquad \widetilde b=P(b-Ag)+g. \]

This leaves $A_{FF}$ unchanged and replaces constrained equations by identities. The formula applies over both real and complex scalar fields; the prescribed coefficients retain their values. The operation preserves symmetry or Hermitian symmetry when present and also applies to nonsymmetric systems.

For example, prescribing $x_0=1$ in

\[ A=\begin{pmatrix}2&3\\0&4\end{pmatrix},\qquad b=\begin{pmatrix}5\\6\end{pmatrix} \]

produces $\widetilde A=\operatorname{diag}(1,4)$ and $\widetilde b=(1,6)^T$ , whether or not the zero at $(1,0)$ is explicitly stored. Discovering constrained rows through their transpose entries would miss the coefficient at $(0,1)$ in this example.

Sparse implementation

The operation validates the system dimensions and shifted constraint indices before mutation, then compresses the operator. It traverses every constrained column of the original matrix to accumulate the load lifting, recording its diagonal position after the traversal. Only after all lifting is complete may an absent diagonal be inserted. Missing diagonals are merged together into the existing CSC arrays by moving entries backwards, avoiding repeated shifts or a full matrix copy. Existing entries and explicit zeros are retained.

A temporary byte mask identifies constrained rows. Consecutive constrained columns are cleared in coalesced storage ranges. Free storage between them is examined in blocks. An empty row mask performs no coefficient stores, while a nonempty mask visits only its set bits. The final scalar tail follows the same zeroing rule. Actual zero assignment or bit clearing is required because multiplication by zero would retain NaNs or create them from infinities. Constraint diagonals and prescribed load entries are restored last.

A single contiguous constraint interval can select sorted row-range pruning. The selection uses sixteen samples at equal stored-entry quantiles, requires at least eight stored entries per column on average, and requires at least half the usable samples to reject the interval. These measured heuristics select the clearing kernel; the column bounds certify each omitted interval. Partial overlaps use binary searches to locate the constrained coefficient range. Fragmented constraints retain the block kernel. No particular FEM boundary numbering is assumed.

For $n$ rows, $m$ stored entries and $k$ constraints, block clearing and batch insertion require $O(n+m+k)$ work and $O(n+k)$ temporary storage. Pruning adds up to two binary searches per overlapping free column. Operator storage grows only for absent constrained diagonals. It remains compressed, with sorted row indices, after the call. The solution vector is untouched, preserving warm starts. Repeating the same constraints is idempotent up to floating-point rounding. An empty constraint map only validates dimensions and compresses storage.

This algorithm is serial and constructs no persistent transpose index or parallel worker pool. The linear scan visits stored entries once while clearing the constrained rows and columns. Tests compare the matrix and load against independent projector constructions, including asymmetric storage, missing diagonals, complex data, offsets and mixed finite-element systems. The opt-in RodinSparseEliminationBenchmarks executable measures the public operation with restoration outside timing.