This specifies the functions and operators that comprise an ODE problem with the form

$\dfrac{dy}{dt} = Ly - B_1 z + r_1(y,t)$

$B_2 y + C z = r_2(x,t)$

where $y$ is the state, $z$ is a constraint force, and $x$ is an auxiliary state describing the constraints.

The optional linear operator L defaults to zeros. The B1 and B2 functions must be of the respective in-place forms B1(dy,z,x,p) (to compute the action of B1 on z) and B2(dz,y,x,p) (to compute the action of B2 on y). The function r1 must of the in-place form r1(dy,y,x,p,t), and r2 must be in the in-place form r2(dz,x,p,t). The C function can be omitted, but if it is included, then it must be of the form C(dz,z,x,p) (to compute the action of C on z). Alternatively, one can supply out-of-place forms, respectively, as B1(z,x,p), B2(y,x,p), C(z,x,p), r1(y,x,p,t) and r2(x,p,t).

An optional keyword argument param_update_func can be used to set a function that updates problem parameters with the current solution. This function must take the in-place form f(q,u,p,t) or out of place form f(u,p,t) to create some q based on u, where y = state(u), z = constraint(u) and x = aux_state(u). (Note that q might enter the function simply as p, to be mutated.) This function can be used to update B1, B2, and C, for example.

We can also include another (unconstrained) set of equations to the set above in order to update x:

$\dfrac{dx}{dt} = r_{1x}(u,p,t)$

In this case, the right-hand side has access to the entire u vector. We would pass the pair of r1 functions as an ArrayPartition.


Construct a saddle-point system operator from the constituent operator blocks. The resulting object can be used with * and \ to multiply and solve. The saddle-point problem has the form

$\begin{bmatrix}A & B_1^T \\ B_2 & C \end{bmatrix} \begin{pmatrix} u \\ f \end{pmatrix} = \begin{pmatrix} r_1 \\ r_2 \end{pmatrix}$


SaddleSystem(A::AbstractMatrix,B₂::AbstractMatrix,B₁ᵀ::AbstractMatrix,C::AbstractMatrix[,eltype=Float64]). Blocks are given as matrices. Must have consistent sizes to stack appropriately. If this is called with SaddleSystem(A,B₂,B₁ᵀ), it sets C to zero automatically.

SaddleSystem(A,B₂,B₁ᵀ,C,u,f[,eltype=Float64]). Operators A, B₂, B₁ᵀ, C are given in various forms, including matrices, functions, and function-like objects. u and f are examples of the data types in the corresponding solution and right-hand side vectors. Guidelines:

  • The entries A and B₂ must be able to act upon u (either by multiplication or as a function) and B₁ᵀ and C must be able to act on f (also, either by multiplication or as a function).
  • A and B₁ᵀ should return data of type u, and B₂ and C should return data of type f.
  • A must be invertible and be outfitted with operators `andldiv!`.
  • Both u and f must be subtypes of AbstractArray: they must be equipped with size and vec functions and with a constructor of the form T(data) where T is the data type of u or f and data is the wrapped data array.

If called as SaddleSystem(A,B₂,B₁ᵀ,u,f), the C block is omitted and assumed to be zero.

If called with SaddleSystem(A,u), this is equivalent to calling SaddleSystem(A,nothing,nothing,u,[]), then this reverts to the unconstrained system described by operator A.

The list of vectors u and f in any of these constructors can be bundled together as a SaddleVector, e.g. SaddleSystem(A,B₂,B₁ᵀ,SaddleVector(u,f)).

An optional keyword argument solver= can be used to specify the type of solution for the Schur complement system. By default, this is set to Direct, and the Schur complement matrix is formed, factorized, and stored. This can be changed to a variety of iterative solvers, e.g. BiCGStabl, CG, GMRES, in which case an iterative solver from IterativeSolvers.jl is used.


Build a vector of the r1 functions for the state ODEs and auxiliary state ODEs.


Build a solution vector for a constrained system. This takes three optional keyword arguments: state, constraint, and aux_state. If only a state is supplied, then the constraint is set to an empty vector and the system is assumed to correspond to an unconstrained system. (aux_state is ignored in this situation.)