Linear static analysis (structsolve.static)

structsolve.static(K, fext, silent=False)

Linear static analysis

Solves the linear system K u = fext using solve(). For linear and

non-linear static analyses based on callables that calculate the force vectors and stiffness matrices, see Analysis.

Parameters:

Ksparse_matrix

Stiffness matrix. Should include initial stress stiffness matrix, aerodynamic matrix and so forth when applicable. It may be unsymmetric, e.g. K0 + kCfollower for a follower load on a geometrically linear structure, whose load stiffness is kCfollower = -dF/dc.

fextarray-like

Vector of external loads.

silentbool, optional

A boolean to tell whether the log messages should be printed.

Returns:

incrementslist

[1.], the load factor of the solution.

cslist

List with the solution vector.

structsolve.solve(a, b, silent=False, **kwargs)

Solve a linear system of equations removing null rows and columns

Wrapper for scipy.sparse.linalg.spsolve(). The rows and columns of a without any non-zero term, typically the constrained degrees of freedom of a Ritz or finite element model, are removed before solving the linear system of equations. The corresponding values of the solution x are zero.

Parameters:

andarray or sparse matrix

A square matrix, converted to CSR form in the solution. It may be unsymmetric, e.g. with the load stiffness of follower loads, since scipy.sparse.linalg.spsolve() uses an LU factorization. The degrees of freedom whose row and column are both null are removed.

bndarray

1-D array representing the right-hand side of the system of equations.

silentbool, optional

A boolean to tell whether the log messages should be printed.

kwargskeyword arguments, optional

Other arguments directly passed to scipy.sparse.linalg.spsolve().

Returns:

xndarray

The solution of the system of equations, a 1-D array of size a.shape[1] with the same dtype as b.

Notes

With \(u\) the set of the degrees of freedom whose row or column of a is not null, the reduced system \([a_{uu}]\boldsymbol{x}_u = \boldsymbol{b}_u\) is solved and \(x = 0\) at the other degrees of freedom. The terms of b at the removed degrees of freedom are ignored: x solves the full system only when they are null, otherwise the full system has no solution, since a null row and column make a singular.

Factorizations of the tangent stiffness matrix

structsolve.static._factorize(kT, linear_solver='splu', constraints=None, cache=None)

Factorization of the tangent stiffness matrix

Null rows and columns are removed, the corresponding values of the solution are zero, as in solve(). The factorization is reused by the non-linear solvers while the tangent stiffness matrix is not updated.

Parameters:

kTsparse matrix

Tangent stiffness matrix.

linear_solverstr, optional

'splu' (default) for the sparse LU factorization of SuperLU, or 'pardiso' for the symmetric indefinite factorization of PARDISO, which requires the shared library mkl_rt of Intel oneMKL, see pardiso and LINEAR_SOLVERS.

constraintsarray-like, sparse matrix or None, optional

The matrix \(\boldsymbol{C}\) of the linear constraints \(\boldsymbol{C}^\top \boldsymbol{x} = \boldsymbol{0}\), of shape (n, m) for m constraints, see the Notes. None (default) for none.

cache_FactorizationCache or None, optional

Cache of the work that depends only on the sparsity pattern of kT, reused by the next factorizations of a matrix of the same pattern, see the Notes. None (default) for none.

Returns:

solve_kTcallable

solve_kT(b, d=None) returns the solution x of kT x = b, or of the system bordered by the constraints, whose values are d, of shape (m,), zeros when None, see the Notes; it has the attributes:

• linear_solver: the linear solver used, 'splu' when 'pardiso' fell back to it for an unsymmetric matrix;

• det_sign: the sign of the determinant of the reduced matrix, 1 or -1, on the space of the constraints when given;

• num_negative_eigvals: the number of negative eigenvalues of the reduced matrix with 'pardiso', on the space of the constraints when given, None with 'splu';

• num_constraints: the number m of constraints of the border;

• num_perturbed_pivots: the number of pivots perturbed by PARDISO, None with 'splu';

• smallest_eigvals, with 'splu' only: smallest_eigvals(k, s=None) returns the k eigenvalues of smallest magnitude of the reduced matrix, scaled by diag(s) on both sides when given, on the space of the constraints, by shift-invert with ARPACK and this factorization, see _smallest_eigvals().

Raises:

RuntimeError

When the matrix is singular, for 'splu', or when PARDISO fails.

ImportError

When linear_solver='pardiso' and oneMKL is not found.

ValueError

For an invalid linear_solver.

Notes

With 'splu', the reduced matrix is factorized once with SuperLU, scipy.sparse.linalg.splu() with its default options, i.e. a sparse LU factorization \(\boldsymbol{P}_r \boldsymbol{A} \boldsymbol{P}_c = \boldsymbol{L} \boldsymbol{U}\) with row and column permutations, valid for unsymmetric matrices. Each call of solve_kT is then only a forward and a backward substitution, which is why the factorization is reused while the tangent stiffness matrix is not updated, e.g. in the modified Newton-Raphson method and in the two solutions of the bordered system of the arc-length methods. \(\boldsymbol{L}\) having a unit diagonal, the sign of the determinant is the product of the signs of the diagonal of \(\boldsymbol{U}\) and of the two permutations; for a symmetric matrix it is \((-1)^{n_{neg}}\), with \(n_{neg}\) its number of negative eigenvalues, which LU does not give.

With 'pardiso', the reduced matrix \(\boldsymbol{A}\) is scaled to a unit diagonal, \(\boldsymbol{A}_s = \boldsymbol{D} \boldsymbol{A} \boldsymbol{D}\), \(\boldsymbol{D} = \mathrm{diag}(|a_{ii}|^{-1/2})\), a null diagonal term being replaced by 1, and the upper triangle of \(\boldsymbol{A}_s\) is factorized by PARDISO as a real symmetric indefinite matrix (mtype=-2), \(\boldsymbol{L} \boldsymbol{D} \boldsymbol{L}^\top\) with Bunch-Kaufman pivoting. By Sylvester’s law of inertia, the number of negative pivots, iparm(23) of PARDISO, is the number of negative eigenvalues of \(\boldsymbol{A}_s\), hence of \(\boldsymbol{A}\), the congruence by \(\boldsymbol{D}\) preserving the inertia: zero on the stable parts of an equilibrium path, it changes at the limit points and at the bifurcation points crossed. PARDISO perturbs the pivots that are too small, iparm(14) of them, which makes the count uncertain when it is not zero. When the first factorization perturbed pivots, the matrix is factorized again with the parameters PARDISO_FALLBACK_IPARM, iparm(24) = 1, the two-level factorization algorithm, kept, solve_kT.pardiso_fallback = True, only when no pivot is perturbed and the relative residual of its solution of a random right-hand side, without refinement, is not larger than PARDISO_ACCURACY_RTOL = 1e-6; otherwise the matrix is factorized by SuperLU, as 'splu', with a RuntimeWarning, and num_negative_eigvals is None: a factorization of iparm(24) = 1 without perturbed pivots had residuals up to 1e3 and a wrong inertia for some Ritz models of the tests (MKL 2025.0). With cache, the next factorizations of the same pattern use the fallback directly, with the phase 12 and the same check, or SuperLU directly. The default factorizations without perturbed pivots are not checked, at no cost. The tangent stiffness matrices of a Ritz model of a cylinder (panels 0.15.0, cylshell_clpt_donnell, 8 x 8 terms, 114 degrees of freedom, a dense matrix) and dense random indefinite matrices of 114 and 200 degrees of freedom had 18 to 40 perturbed pivots with the default parameters (MKL 2026.1, iparm(10) = 8, iparm(11) = iparm(13) = 0, iparm(21) = 1, iparm(24) = 0), solutions of relative residual 0.4 to 1, recovered by the GMRES below, and a number of negative pivots short by one, e.g. 0 from the load factor 0.827 of the critical point, where a dense eigenvalue analysis has one negative eigenvalue, up to 0.903; neither iparm(10) = 13 or 20, nor the scaling and the weighted matching, iparm(11) = iparm(13) = 1, nor iparm(21) = 3 changed it, and iparm(24) = 1 gave no perturbed pivot, the counts of the dense eigenvalue analysis along the whole path, and residuals of 1e-13. The tangent of the DOE09 cylinder below, sparse, had no perturbed pivot. The diagonal of the tangent of a model with derivatives of the displacements among its unknowns may span many orders of magnitude, which the scaling compensates. Every solution is checked against its relative residual \(\|\boldsymbol{A} \boldsymbol{x} - \boldsymbol{b}\|/\|\boldsymbol{b}\|\) and, when it is larger than PARDISO_REFINE_TOL, improved by GMRES with the factorization as the preconditioner, to the relative tolerance PARDISO_GMRES_RTOL; a RuntimeWarning is issued when the residual remains larger than PARDISO_WARN_TOL. These are the settings of the Riks analyses of a tow-steered cylinder for the DCSE website (2026), where the unsymmetric LU of PARDISO was wrong and its Cholesky factorization failed beyond the limit point. The symmetric factorization would be wrong for an unsymmetric matrix, e.g. with the load stiffness of a follower load: 'splu' is then used, with a RuntimeWarning. With constraints, the bordered system

\[\begin{split}\begin{bmatrix} \boldsymbol{0} & \boldsymbol{C}^\top \\ \boldsymbol{C} & \boldsymbol{K}_T \end{bmatrix} \begin{Bmatrix} \boldsymbol{\mu} \\ \boldsymbol{x} \end{Bmatrix} = \begin{Bmatrix} \boldsymbol{d} \\ \boldsymbol{b} \end{Bmatrix}\end{split}\]

is factorized instead, the border first, such that no row of the upper triangle factorized by PARDISO is empty, and solve_kT returns \(\boldsymbol{x}\), which satisfies \(\boldsymbol{C}^\top \boldsymbol{x} = \boldsymbol{d}\), \(\boldsymbol{\mu}\) being the Lagrange multipliers, the reactions of the constraints, \(\boldsymbol{d} = \boldsymbol{0}\) unless d is given, e.g. a prescribed displacement. The rows of \(\boldsymbol{C}\) of the null rows and columns of kT are removed, and the constraints left without any non-zero term, satisfied with \(d_j = 0\), are dropped; a non-zero \(d_j\) of such a constraint raises a ValueError, the degrees of freedom of the null rows and columns being zero. The bordered matrix is indefinite even when \(\boldsymbol{K}_T\) is positive definite on the space of the constraints: for \(\boldsymbol{C}\) of full rank its inertia is that of \(\boldsymbol{Z}^\top \boldsymbol{K}_T \boldsymbol{Z}\) plus \(m\) positive and \(m\) negative eigenvalues, \(\boldsymbol{Z}\) a basis of the null space of \(\boldsymbol{C}^\top\) (Gould, N.I.M., On practical conditions for the existence and uniqueness of solutions to the general equality quadratic programming problem, Mathematical Programming, 32, 90-99, 1985), hence num_negative_eigvals is the number of negative pivots minus \(m\), and det_sign is multiplied by \((-1)^m\). This is the border of the inertia relief of the Riks analyses of a tow-steered cylinder for the DCSE website (2026), with a single constraint.

PARDISO is called through PardisoSolver, whose memory is freed when solve_kT is deleted.

With cache, created by Analysis.static() for each non-linear analysis, along which the pattern of the tangent stiffness matrix usually does not change, the work that depends only on the pattern is done once per pattern: the indices of the reduced and bordered matrix, of its CSC form, of the transposed terms of the symmetry check, of its diagonal and of the scaled upper triangle, from which each factorization only gathers the values of kT; and the analysis (reordering and symbolic factorization) of PARDISO, its phase 11, the next factorizations being its phase 22 instead of 12 with a solver that keeps the analysis of the pattern; up to CACHE_MAX_PARDISO_SOLVERS such solvers are kept, since the factorizations of the last converged state and of the current iteration coexist. The pattern is the shape, indptr and indices of kT in CSR form, compared at every call, the degrees of freedom whose row or column is not null, which depend on the values, since explicit zeros do not count, computed at every call, and constraints, compared with its values; any change rebuilds the cache. The values enter only through the gathered terms, and the symmetry check, the scaling and the null rows and columns are evaluated at every call, such that the factorization is bitwise the one without the cache, verified by the tests, also the inertia and the perturbed pivots of PARDISO (with MKL_NUM_THREADS=1; the analysis of PARDISO depends only on the pattern, with its default iparm(11) = iparm(13) = 0 for symmetric indefinite matrices). A matrix that is not float64 or whose CSR form has duplicated or unsorted terms is factorized without the cache. On a symmetric matrix of 57000 degrees of freedom with 4.5 million terms, the factorization took 1.65 s without and 0.40 s with the cache with 'pardiso', with or without one constraint, and 1.1 s with 'splu', whose preprocessing went from 0.11 s to 0.05 s, SuperLU computing its ordering at every factorization.

structsolve.static._smallest_eigvals(lu, m, k, s=None)

The k eigenvalues of smallest magnitude of the matrix of the SuperLU factorization lu, on the space of its m constraints

The matrix is the reduced matrix \(\boldsymbol{K}\) of _factorize(), bordered first by m constraints \(\boldsymbol{C}\). Shift-invert at \(\sigma = 0\) with ARPACK, scipy.sparse.linalg.eigs(), on the operator \(\boldsymbol{T} \boldsymbol{b} = \boldsymbol{x}\), the solution of the bordered system \(\boldsymbol{K} \boldsymbol{x} + \boldsymbol{C} \boldsymbol{\mu} = \boldsymbol{b}\), \(\boldsymbol{C}^\top \boldsymbol{x} = \boldsymbol{0}\), by the forward and backward substitutions of lu, which is not factorized again: an eigenvalue \(\nu \ne 0\) of \(\boldsymbol{T}\) has its eigenvector \(\boldsymbol{b}\) in the null space of \(\boldsymbol{C}^\top\), and \(\boldsymbol{P} \boldsymbol{K} \boldsymbol{b} = \boldsymbol{b}/\nu\), \(\boldsymbol{P}\) the orthogonal projection on that space, i.e. \(\lambda = 1/\nu\) is an eigenvalue of \(\boldsymbol{Z}^\top \boldsymbol{K} \boldsymbol{Z}\), \(\boldsymbol{Z}\) an orthonormal basis of the null space of \(\boldsymbol{C}^\top\); the range of \(\boldsymbol{C}\), of \(\nu = 0\), is not among the largest \(|\nu|\). With k larger than the size minus 2, the limit of ARPACK, all the eigenvalues are computed with numpy.linalg.eigvals() from \(\boldsymbol{T}\) built column by column.

With the scaling s, of the reduced degrees of freedom, the eigenvalues are those of \(\boldsymbol{S} \boldsymbol{K} \boldsymbol{S}\), \(\boldsymbol{S} = \mathrm{diag}(s)\), with the constraints \(\boldsymbol{S} \boldsymbol{C}\), from the operator \(\boldsymbol{S}^{-1} \boldsymbol{T} \boldsymbol{S}^{-1}\): those of \(\boldsymbol{K} \boldsymbol{v} = \lambda \boldsymbol{S}^{-2} \boldsymbol{v}\), with the positive diagonal reference matrix \(\boldsymbol{S}^{-2}\), see newton_raphson._eigvals_inertia(). smallest_eigvals(k, s) of _factorize() takes s of all the degrees of freedom.

Returns:

lamndarray or None

The eigenvalues, complex, sorted by increasing magnitude, None when ARPACK did not converge.

structsolve.static._pardiso_factorize_fallback(A, iparm)

Factorization of phase 12 of A by a new PardisoSolver with the parameters iparm (the defaults filled by a previous factorization), iparm(1) = 1, and those of PARDISO_FALLBACK_IPARM; None when it fails

structsolve.static.PARDISO_FALLBACK_IPARM = {24: 1}

'pardiso': parameters, {i: iparm(i)}, of the second factorization of a matrix whose first factorization, with the default parameters, perturbed pivots: iparm(24) = 1, the two-level factorization algorithm instead of the classic one, see _factorize()