Class NetworkLinearSolver

java.lang.Object
neqsim.process.equipment.network.NetworkLinearSolver
public class NetworkLinearSolver extends Object
Sparse and dense linear system solvers for pipeline network equations.

Replaces the O(n³) Gaussian elimination in the Newton-Raphson GGA solver with efficient alternatives from the EJML library (already a NeqSim dependency):

  • Sparse LU: For large networks (>50 nodes), uses compressed sparse column (CSC) format with fill-reducing ordering. Complexity is O(nnz) amortized where nnz is the number of non-zeros in the Schur complement matrix.
  • Dense LU: For small-to-medium networks (≤50 nodes), uses EJML's optimized dense LU with partial pivoting. Faster than hand-coded Gaussian for n > 10.

The Schur complement matrix in the Todini-Pilati GGA is structurally sparse: entry (i,j) is nonzero only if nodes i and j are connected by at least one pipe. For a gas gathering network with n=100 nodes and average degree 3, matrix density is ~6%, making sparse solvers 10–50x faster than dense Gaussian.

Version:
1.0
Author:
Even Solbraa
See Also:

  • Field Summary

    Fields
    Modifier and Type
    Field
    Description
    private static final int
    EJML_THRESHOLD
    Threshold for switching from Gaussian to EJML solvers.
    private static final org.apache.logging.log4j.Logger
    logger
     
    private static final int
    SPARSE_THRESHOLD
    Threshold for switching from dense to sparse EJML solver within the EJML path.

  • Constructor Summary

    Constructors
    Constructor
    Description
    NetworkLinearSolver()
     

  • Method Summary

    Modifier and Type
    Method
    Description
    static double[]
    estimateSparsity(int nodeCount, int pipeCount)
    Estimate the sparsity pattern of the network Schur complement.
    static double[]
    solve(double[][] matA, double[] vecB, int n)
    Solve the linear system Ax = b using the most appropriate method.
    static double[]
    solveDense(double[][] matA, double[] vecB, int n)
    Solve using EJML dense LU decomposition with partial pivoting.
    static double[]
    solveGaussian(double[][] matAOrig, double[] vecBOrig, int n)
    Fallback: Gaussian elimination with partial pivoting.
    static double[]
    solveSparse(double[][] matA, double[] vecB, int n)
    Solve using EJML sparse CSC LU decomposition.

    Methods inherited from class Object

    clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

  • Field Details

    • logger

      private static final org.apache.logging.log4j.Logger logger

    • EJML_THRESHOLD

      private static final int EJML_THRESHOLD
      Threshold for switching from Gaussian to EJML solvers. Networks with more free nodes than this threshold use EJML (dense LU for n ≤ 100, sparse CSC LU for n > 100). Below this threshold, Gaussian elimination is used for backward compatibility with existing solver convergence behavior.
      See Also:

    • SPARSE_THRESHOLD

      private static final int SPARSE_THRESHOLD
      Threshold for switching from dense to sparse EJML solver within the EJML path.
      See Also:

  • Constructor Details

    • NetworkLinearSolver

      public NetworkLinearSolver()

  • Method Details

    • solve

      public static double[] solve(double[][] matA, double[] vecB, int n)
      Solve the linear system Ax = b using the most appropriate method.

      For small systems (n ≤ 30), uses Gaussian elimination with partial pivoting for backward compatibility. For medium systems, uses EJML dense LU. For large systems (n > 100), uses EJML sparse CSC LU. Falls back to Gaussian elimination if EJML solvers fail.

      Parameters:
      matA - coefficient matrix (n x n)
      vecB - right-hand side vector (n)
      n - system size
      Returns:
      solution vector x

    • solveDense

      public static double[] solveDense(double[][] matA, double[] vecB, int n)
      Solve using EJML dense LU decomposition with partial pivoting.

      Uses EJML LinearSolverFactory_DDRM.lu() which provides O(n³/3) factorization with BLAS-optimized operations. For n=50, approximately 5x faster than hand-coded Gaussian elimination.

      Parameters:
      matA - coefficient matrix (n x n)
      vecB - right-hand side vector (n)
      n - system size
      Returns:
      solution vector x

    • solveSparse

      public static double[] solveSparse(double[][] matA, double[] vecB, int n)
      Solve using EJML sparse CSC LU decomposition.

      Converts the dense Schur complement to compressed sparse column (CSC) format, discarding structural zeros. Uses fill-reducing column ordering and sparse LU factorization. For a 100x100 matrix with 6% density, approximately 10-50x faster than dense Gaussian.

      Parameters:
      matA - coefficient matrix (n x n) — may have many zeros
      vecB - right-hand side vector (n)
      n - system size
      Returns:
      solution vector x

    • solveGaussian

      public static double[] solveGaussian(double[][] matAOrig, double[] vecBOrig, int n)
      Fallback: Gaussian elimination with partial pivoting.

      This is the original O(n³) solver, kept as a robust fallback when EJML solvers encounter issues (singular matrices, numerical edge cases). Makes a defensive copy of the input arrays.

      Parameters:
      matAOrig - coefficient matrix (n x n) — not modified
      vecBOrig - right-hand side vector (n) — not modified
      n - system size
      Returns:
      solution vector x

    • estimateSparsity

      public static double[] estimateSparsity(int nodeCount, int pipeCount)
      Estimate the sparsity pattern of the network Schur complement.

      For diagnostics: Returns the estimated density and recommended solver type.

      Parameters:
      nodeCount - number of free nodes
      pipeCount - number of pipe elements
      Returns:
      array [density, nonzeros, recommended_threshold] where recommended_threshold is 0 for dense and 1 for sparse