Port the Stiefel manifold solvers to Fortran 77 - #5
Merged
Conversation
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>
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 ☂️ |
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>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Moves the
pystmoptalgorithm 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.
src/smblas.fsrc/smman.fC,JA,JC, theAmap, the polar retractionsrc/smprox.fsrc/smslpg.fsrc/smpencf.fsrc/_smopt.pyfsrc/smopt/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/N803hold everywhere and the naming exemptions are gone frompyproject.toml):SLPG_smoothslpg_smoothSLPGslpgSLPG_l21slpg_l21PenCFpencfStiefel.Phi/A/C/JA/JC/JC_transposephi/a/c/ja/jc/jc_transposeFeas_eval/Init_point/Post_processfeas_eval/init_point/post_processXinit=xinit=Stiefelkeeps its capital, being a class. The short method names are the symbols the theory documentation already uses for those maps. The original's top-levelSolver/manifold/utilitypackages becomesmopt.solver/smopt.manifold/smopt.utility, since shipping those names at top level would collide with anything else installed.Defects fixed in the original
Init_pointcompared an array againstNonewith==, which raises rather than testing for a missing argument.Post_processdivided 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.--build-dirbroke on any path containing a space — the normal case on Windows.bin/run_f2py.pystages the signature file into the build directory and runs f2py on a bare file name.-static-*flags. The link language is pinned to Fortran, and the flags are probed rather than assumed, since-static-libquadmathonly exists from GCC 13 on.DLL load failed while importing _smopt. It still depended onlibwinpthread-1.dll:libgfortran.arefers to it and the gfortran driver appends its own-lwinpthreadafter everything Meson passes, by which point the linker is back in dynamic mode.--whole-archivedefines those symbols independently of link order. The built extension now importsKERNEL32, the UCRT stubs andpython3xxonly. This affected the published Windows wheel too, not just CI: the build job passed because cibuildwheel builds without importing.licensecannot be combined with aLicense ::classifier; meson-python refuses outright.python-versionaxis was referenced only in the artifact name and never passed tosetup-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_minimumassertedrel=1e-5against 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_errorwould 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 whatX'X = Ilooks like.mkdocs.ymlalready pointed at four asset paths that did not exist (smopt.ico,smopt.png,extra.css,katex.js); all are now resolved, with a singlesmopt-icon.svgserving as both logo and favicon.Verification
tests/reference.pykeeps 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 checkandruff format --checkare 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.