Add docstrings

src/buffer.jl

@@ -81,13 +81,42 @@ function _add_vec!(I, J, V, j, offset, A::SparseArrays.SparseMatrixCSC)
     return offset + length(Ai)
 end
 
-# `SparseMatrixCSC` is stored with an offset by column.
-# This means that getting view `view(A, :, I)` can be handles efficently,
-# these give `SparseMatrixCSCView` (if `I` is a `UnitRange`) and
-# `SparseMatrixCSCColumnSubset` otherwise.
-# In `schur.jl`, we therefore get a `SparseMatrixCSCColumnSubset`.
-# Since we want to use subsets of constraint indices, we use the columns
-# of `A` for constraint indices and the rows of `A` for matrix indices.
+"""
+    buffer_for_jprod(model::Model, i::MatrixIndex)
+
+Return the sparse matrix collecting the vectorization of every constraint
+matrix of the `i`th PSD block:
+```
+𝐀ᵢ = [vec(Aᵢ₁) vec(Aᵢ₂) … vec(Aᵢₙ)] ∈ ℝ^(mᵢ² × n)
+```
+where `mᵢ` is the side dimension of the block and `n` is the number of
+constraints. This is the matrix `𝒜` of [HKS24, Section 3.1], which is
+stated there for a single PSD block, here built once per block.
+
+It represents the linear equality-constraint operator
+`Xᵢ ↦ (⟨Aᵢⱼ, Xᵢ⟩)ⱼ` as a single matrix, so that `add_jprod!` is one
+sparse matrix-vector product. It is computed once at problem setup and
+reused at every interior-point iteration, both for the Jacobian products
+and for the dense columns of the Schur complement assembled in
+`schur.jl`.
+
+`SparseMatrixCSC` is stored with an offset by column.
+This means that getting view `view(A, :, I)` can be handles efficently,
+these give `SparseMatrixCSCView` (if `I` is a `UnitRange`) and
+`SparseMatrixCSCColumnSubset` otherwise.
+In `schur.jl`, we therefore get a `SparseMatrixCSCColumnSubset`.
+Since we want to use subsets of constraint indices, we use the columns
+of `A` for constraint indices and the rows of `A` for matrix indices.
+The subsets used in `schur.jl` are suffixes of the constraints sorted by
+decreasing number of nonzeros, following [FKN97].
+
+[HKS24] Habibi, Soodeh, Michal Kočvara, and Michael Stingl. "Loraine -- an
+interior-point solver for low-rank semidefinite programming."
+Optimization Methods and Software 39.6 (2024): 1185-1215.
+[FKN97] Fujisawa, Katsuki, Masakazu Kojima, and Kazuhide Nakata.
+"Exploiting sparsity in primal-dual interior-point methods for
+semidefinite programming." Mathematical Programming 79 (1997): 235-253.
+"""
 function buffer_for_jprod(model::Model{T}, i::MatrixIndex) where {T}
     nnz = sum(1:(model.meta.ncon); init = 0) do j
         return _nnz(model.A[i.value, j])
@@ -124,6 +153,21 @@ function _add_jprod!(V, Jv::AbstractArray, A)
     return _add_mul!(Jv, A', _vec(V), true)
 end
 
+"""
+    add_sub_jprod!(
+        model::BufferedModelForSchur,
+        i::MatrixIndex,
+        V::AbstractMatrix,
+        Jv::AbstractVector,
+        I,
+    )
+
+Same as `add_jprod!` but restricted to the constraints of `I`:
+adds `⟨Aᵢⱼ, V⟩` to `Jv[k]` for the `k`th entry `j` of `I`. This is the
+product with the column subset `𝐀ᵢ[:, I]` of the buffer described in
+`buffer_for_jprod`, which is why the buffer is stored with
+constraints as columns.
+"""
 function add_sub_jprod!(
     model::BufferedModelForSchur,
     i::MatrixIndex,
@@ -136,6 +180,20 @@ function add_sub_jprod!(
     return _add_jprod!(V, Jv, view(A, :, I))
 end
 
+"""
+    add_jprod!(
+        model::BufferedModelForSchur,
+        V::AbstractMatrix,
+        Jv::AbstractVector,
+        i::MatrixIndex,
+    )
+
+Add the contribution of the `i`th PSD block to the product between the
+Jacobian of the equality constraints and `V`, that is, add `𝐀ᵢᵀ vec(V)`
+to `Jv`, whose `j`th entry is `⟨Aᵢⱼ, V⟩`. Here `𝐀ᵢ` is the buffer built
+by `buffer_for_jprod` so this is a single sparse matrix-vector
+product.
+"""
 function add_jprod!(
     model::BufferedModelForSchur,
     V::AbstractMatrix,

src/schur.jl

@@ -1,6 +1,7 @@
 # This code computes the Schur complement using the ideas detailed in [FKN97, Section 3]
 # This is useful to compute search direction in primal-dual interior-point methods for semidefinite programs [FKN97]
 # [FKN97] Fujisawa, Katsuki, Masakazu Kojima, and Kazuhide Nakata. "Exploiting sparsity in primal-dual interior-point methods for semidefinite programming." Mathematical Programming 79 (1997): 235-253.
+# [HKS24] Habibi, Soodeh, Michal Kočvara, and Michael Stingl. "Loraine -- an interior-point solver for low-rank semidefinite programming." Optimization Methods and Software 39.6 (2024): 1185-1215.
 # It was adapted from dapted from Michal Kocvara's code in
 # https://github.com/kocvara/Loraine.jl/blob/bd2821ba830786a78f04081d7e8f5cac25e56cac/src/makeBBBB.jl
 
@@ -81,6 +82,16 @@ function add_schur_complement!(
     return H
 end
 
+# Adds the contribution `𝐀ᵢᵀ (W ⊗ W) 𝐀ᵢ` of the PSD block `mat_idx` to `H`,
+# where `𝐀ᵢ` is the buffer documented in `buffer_for_jprod`.
+# `H` is assembled column by column. For a constraint `i` in the dense
+# regime, that is, one of the `last_dense[ilmi]` first entries of the
+# permutation `σ` sorting the constraints by decreasing number of nonzeros,
+# the dense matrix `W Aᵢ W` is formed once and all the remaining inner
+# products `⟨Aⱼ, W Aᵢ W⟩`, `j ∈ I`, are then obtained with the single
+# sparse matrix-vector product `𝐀ᵢ[:, I]ᵀ vec(W Aᵢ W)` of
+# `add_sub_jprod!`. The constraints with at most one nonzero use the
+# dedicated low-nnz paths below, also following [FKN97].
 # /!\ W needs to be symmetric
 function add_schur_complement!(
     model::BufferedModelForSchur,
@@ -162,6 +173,9 @@ end
 # [HKS24, (5b)]
 # Returns the matrix equal to the sum, for each equation, of
 # ⟨A_i, WA_jW⟩
+# [HKS24, (5b)] is stated for a single PSD block; summing over the blocks
+# and adding the contribution of the scalar block gives
+# H = ∑ᵢ 𝐀ᵢᵀ (Wᵢ ⊗ Wᵢ) 𝐀ᵢ + C_lin * Diagonal(w) * C_linᵀ
 function schur_complement!(model::BufferedModelForSchur, W::AbstractVector, H)
     fill!(H, zero(eltype(H)))
     if num_matrices(model) > 0