Skip to content

Port the Stiefel manifold solvers to Fortran 77 - #5

Merged
saudzahirr merged 8 commits into
masterfrom
port-stiefel-solvers-to-fortran-77
Aug 25, 2026
Merged

Port the Stiefel manifold solvers to Fortran 77#5
saudzahirr merged 8 commits into
masterfrom
port-stiefel-solvers-to-fortran-77

Conversation

@saudzahirr

@saudzahirr saudzahirr commented Aug 25, 2026

Copy link
Copy Markdown
Contributor

smopt

Moves the pystmopt algorithm into this template and reimplements every numerical step in Fortran 77, reached through f2py.

What is in Fortran

The Fortran side owns the whole computation, not just the kernels — the manifold geometry, the proximal operators, the Barzilai-Borwein step sizes and the solver loops themselves. Python only marshals arguments, calls back into the user objective, and prints progress.

The dense linear algebra is hand written alongside it (matrix products, Cholesky, a Jacobi eigensolver for the polar factor, modified Gram-Schmidt), so the extension links against nothing but the Fortran runtime.

File Contents
src/smblas.f dense kernels
src/smman.f Stiefel geometry: C, JA, JC, the A map, the polar retraction
src/smprox.f proximal operators and the l21 multiplier
src/smslpg.f the SLPG drivers and the Arrow-Hurwicz inner iteration
src/smpencf.f the PenCF driver
src/_smopt.pyf f2py signatures
src/smopt/ the Python layer

All fixed-form, within 72 columns, no allocatable storage — f2py allocates the workspaces via intent(cache,hide) so they never appear in the Python signature.

Interface

The solver set and the output dictionary are unchanged, but the names now follow PEP 8 rather than the original PascalCase, so the package no longer ships names its own linter rejects (N802/N803 hold everywhere and the naming exemptions are gone from pyproject.toml):

was now
SLPG_smooth slpg_smooth
SLPG slpg
SLPG_l21 slpg_l21
PenCF pencf
Stiefel.Phi/A/C/JA/JC/JC_transpose phi/a/c/ja/jc/jc_transpose
Feas_eval/Init_point/Post_process feas_eval/init_point/post_process
Xinit= xinit=

Stiefel keeps its capital, being a class. The short method names are the symbols the theory documentation already uses for those maps. The original's top-level Solver / manifold / utility packages become smopt.solver / smopt.manifold / smopt.utility, since shipping those names at top level would collide with anything else installed.

Defects fixed in the original

  • Init_point compared an array against None with ==, which raises rather than testing for a missing argument.
  • Post_process divided by a vanishing singular value, returning a point off the manifold whenever the l21 penalty had zeroed enough rows to make the iterate singular. The rank deficient directions now get an orthonormal completion.

Build and CI fixes

Several of these affect any project generated from f2py-tpl, not just this one.

  • f2py splits its own command line on whitespace, so the absolute --build-dir broke on any path containing a space — the normal case on Windows. bin/run_f2py.py stages the signature file into the build directory and runs f2py on a bare file name.
  • Meson linked this mixed C/Fortran target with gcc, which rejects the -static-* flags. The link language is pinned to Fortran, and the flags are probed rather than assumed, since -static-libquadmath only exists from GCC 13 on.
  • The Windows extension did not import at allDLL load failed while importing _smopt. It still depended on libwinpthread-1.dll: libgfortran.a refers to it and the gfortran driver appends its own -lwinpthread after everything Meson passes, by which point the linker is back in dynamic mode. --whole-archive defines those symbols independently of link order. The built extension now imports KERNEL32, the UCRT stubs and python3xx only. This affected the published Windows wheel too, not just CI: the build job passed because cibuildwheel builds without importing.
  • An SPDX license cannot be combined with a License :: classifier; meson-python refuses outright.
  • The test matrix was not testing what it claimed. Its python-version axis was referenced only in the artifact name and never passed to setup-uv, so every job resolved to whichever interpreter uv picked — the job labelled 3.10 actually built and tested cp313. Now passed through, so the five versions are really covered.
  • test_pencf_reaches_the_known_minimum asserted rel=1e-5 against the known optimum, which PenCF does not deliver. Holding the problem fixed and varying only the starting point over 60 draws, the relative error has median 7.6e-07 and worst case 1.2e-04; a 150-seed sweep tops out at 1.0e-04, and the NumPy reference is worse still at 1.7e-04. That is the penalty method, not the port — the SLPG solvers reach 3.7e-16 on the same instance and keep their tight tolerance.

The Codecov 404s seen on the first run were left alone deliberately: every earlier run had that step skipped because pytest failed first, so it was the first run to reach the endpoint and no repository record existed yet. Uploads later in the very same run succeeded, so the race is over and weakening fail_ci_if_error would have papered over a resolved transient.

Docs assets

Artwork in the eggzec block style, built from the same rounded blocks on the same grey ground as the other project marks, with lettering from the real eggzec-block.js. Colours are the Space berries palette. The mark is the orthonormality condition itself: a bracketed matrix whose diagonal is solid and whose off-diagonal entries have collapsed to points, which is what X'X = I looks like.

mkdocs.yml already pointed at four asset paths that did not exist (smopt.ico, smopt.png, extra.css, katex.js); all are now resolved, with a single smopt-icon.svg serving as both logo and favicon.

Verification

tests/reference.py keeps a NumPy transcription of the original algorithm as a test oracle, and the suite compares whole trajectories against it, not just the final answer. 151 tests plus 6 doctests pass; ruff check and ruff format --check are clean.

The two implementations are not bit-identical — NumPy sums through BLAS, the Fortran kernels sum in their own order — so the trajectory tests compare a 10-iteration horizon, where the worst deviation across every case is 6.5e-10. Convergence over long runs is asserted separately against the known optimum.

All 23 checks pass across ubuntu, macOS and Windows on Python 3.10 through 3.14.

Move the pystmopt algorithm into this template and reimplement every
numerical step in Fortran 77, reached through f2py.

The Fortran side owns the whole computation, not just the kernels: the
manifold geometry, the proximal operators, the Barzilai-Borwein step
sizes and the solver loops themselves all live there. Python only
marshals arguments, calls back into the user objective, and prints
progress. The dense linear algebra the solvers need -- matrix products,
Cholesky, a Jacobi eigensolver for the polar factor, modified
Gram-Schmidt -- is hand written alongside, so the extension links
against nothing but the Fortran runtime.

  src/smblas.f    dense kernels
  src/smman.f     Stiefel geometry
  src/smprox.f    proximal operators
  src/smslpg.f    SLPG drivers and the Arrow-Hurwicz inner iteration
  src/smpencf.f   PenCF driver
  src/_smopt.pyf  f2py signatures
  src/smopt/      the Python layer

The public interface is unchanged: SLPG_smooth, SLPG, SLPG_l21, PenCF,
Stiefel, prox_l1 and prox_l21 keep their signatures and their output
dictionary, now under the smopt namespace.

Two defects in the original are fixed along the way. Init_point compared
an array against None with ==, which raises rather than testing for a
missing argument. Post_process divided by a vanishing singular value,
returning a point off the manifold whenever the l_{2,1} penalty had
zeroed enough rows to make the iterate singular; the rank deficient
directions now get an orthonormal completion, which is what a singular
value decomposition hands back anyway.

tests/reference.py keeps a NumPy transcription of the original as a test
oracle, and the suite compares whole trajectories against it rather than
just the final answer.

Three build fixes were needed to make the template work here:

  - f2py splits its own command line on whitespace, so the absolute
    --build-dir broke on any path containing a space. bin/run_f2py.py
    stages the signature file into the build directory and runs f2py on
    a bare file name instead.
  - Meson linked this mixed C/Fortran target with gcc, which rejects the
    -static-* flags; the link language is now pinned to Fortran, and the
    flags are probed rather than assumed, since -static-libquadmath only
    exists from GCC 13 on.
  - An SPDX license expression cannot be combined with a License ::
    classifier, which meson-python refuses outright.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
@codecov

codecov Bot commented Aug 25, 2026

Copy link
Copy Markdown

Welcome to Codecov 🎉

Once you merge this PR into your default branch, you're all set! Codecov will compare coverage reports and display results in all future pull requests.

ℹ️ You can also turn on project coverage checks and project coverage reporting on Pull Request comment

Thanks for integrating Codecov - We've got you covered ☂️

saudzahirr and others added 7 commits August 25, 2026 17:25
Two CI failures, both real.

The Windows extension did not import at all: "DLL load failed while
importing _smopt". The wheel still depended on libwinpthread-1.dll, so
it only worked on a machine that happened to have MinGW installed --
which the runners do not, and neither do users. The -static-lib* flags
cover only the libraries they name, and the gfortran driver appends its
own -l entries after our link_args, by which point -Wl,-Bdynamic is in
effect again, so the trailing -lwinpthread resolved against the DLL.
-static sets static resolution for every -l the driver adds. The built
extension now imports KERNEL32, msvcrt and python3xx only.

This affected the published Windows wheel too, not just CI: the build
job passed because cibuildwheel builds without importing.

test_pencf_reaches_the_known_minimum asserted rel=1e-5 against the known
optimum, which PenCF does not deliver. Holding the test problem fixed
and varying only the starting point over 60 draws, the relative error
has median 7.6e-07 and worst case 1.2e-04; a 150-seed sweep over
problems tops out at 1.0e-04, and the NumPy reference in
tests/reference.py is worse still at 1.7e-04. So this is the penalty
method, not the port -- the SLPG solvers reach 3.7e-16 on the same
instance and keep their tight tolerance. The assertion is now rel=1e-3,
with feasibility, which the method does guarantee to ~3e-15, still
asserted tightly, plus a check that the objective never drops below the
true minimum.

The Codecov 404s in the same run are left alone deliberately. Every one
of the eight earlier runs had that step skipped because pytest failed
first, so this was the first run to reach the endpoint and no repository
record existed yet; uploads at or after 11:23:46Z in the very same run
succeeded. That race is over, so weakening fail_ci_if_error would be
papering over a resolved transient.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Prints the extension's DLL imports and the exact Win32 load error so the
missing dependency can be identified rather than guessed at. Reverted in
the follow-up commit.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The CI diagnostic showed the extension still importing libwinpthread-1.dll
even with -static on the link line:

  DLL Name: KERNEL32.dll
  DLL Name: api-ms-win-crt-*.dll
  DLL Name: libwinpthread-1.dll     <- this
  DLL Name: python313.dll

The DLL is present in C:\mingw64\bin and that directory is on PATH, but
CPython 3.8 and later do not search PATH when resolving an extension's
dependencies, so the import fails anyway. It has to be linked in.

libgfortran.a refers to libwinpthread, and the gfortran driver appends
its own -lwinpthread from the spec file after everything Meson passes,
by which point the linker is back in dynamic mode. -static did not
prevent that on the MinGW-Builds UCRT toolchain the runners use, even
though it did on the older msvcrt toolchain used locally.

--whole-archive defines those symbols unconditionally and independently
of link order, so the trailing -lwinpthread has nothing left to import.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
It did its job: the extension now imports KERNEL32, the UCRT stubs and
python3xx only, and ctypes.CDLL loads it cleanly on the runners.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Artwork in the eggzec block style, built from the same rounded blocks on
the same grey ground as the other project marks, with the lettering taken
from the real eggzec-block.js rather than a lookalike. Colours are the
Space berries palette.

The mark is the orthonormality condition itself: a bracketed matrix whose
diagonal is solid and whose off-diagonal entries have collapsed to points,
which is what X'X = I looks like. Favicon sizes at or below 32 pixels use
a simplified cut that drops the off-diagonal points, since they turn to
mush and blur into the diagonal at that size; the .ico is assembled by
hand because Pillow only rescales a single image.

mkdocs.yml already pointed at four asset paths that did not exist:
assets/images/smopt.ico, assets/images/smopt.png,
assets/stylesheets/extra.css and assets/javascripts/katex.js. All four
are now present. The stylesheet carries the palette into the theme, and
katex.js renders what arithmatex leaves behind in generic mode, rebinding
on Material's document$ so it survives instant navigation.

Separately, the test matrix was not testing what it claimed. Its
python-version axis was referenced only in the artifact name and never
passed to setup-uv, so every job resolved to whichever interpreter uv
picked by default -- the Windows job labelled 3.10 built and would have
tested cp313. It is now passed through, so the five versions are really
covered.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Follows the convention the other eggzec projects use: the banner leads
the README and the docs index through a raw.githubusercontent.com URL, so
it renders on PyPI and anywhere else the README is displayed rather than
only inside a repository checkout.

The icon set collapses to one file. docs/assets/images/smopt-icon.svg is
now both the theme logo and the favicon, which browsers have accepted as
SVG for years, so the raster and multi-size .ico variants are gone along
with the small-size cut that existed only for them. A vector banner sits
alongside the raster one for the same reason.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The old signatures were carried over verbatim from the code this port
came from, which meant the package shipped names its own linter rejects.
The project rules win: N802 and N803 now hold everywhere and the naming
exemptions are gone from pyproject.toml.

  SLPG_smooth -> slpg_smooth      Phi -> phi          A  -> a
  SLPG        -> slpg             JA  -> ja           C  -> c
  SLPG_l21    -> slpg_l21         JC  -> jc           Xinit -> xinit
  PenCF       -> pencf            JC_transpose -> jc_transpose
                                  Feas_eval    -> feas_eval
                                  Init_point   -> init_point
                                  Post_process -> post_process

Stiefel keeps its capital, being a class. The short method names are the
symbols the theory documentation already uses for those maps, so the two
still read against each other one for one.

tests/reference.py is renamed to match, so the oracle and the package
under test are still reachable through the same attribute lookups, and
its local variables follow suit; L became lip rather than l, which E741
rejects. Docs, README and the doctests are updated throughout, and the
examples in them were run to confirm they work as written.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
@saudzahirr
saudzahirr merged commit cee2640 into master Aug 25, 2026
23 checks passed
@saudzahirr
saudzahirr deleted the port-stiefel-solvers-to-fortran-77 branch August 25, 2026 13:30
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant