Repository navigation
Conversation
Add src/util/nlopt-blas.h, through which algorithms call the dense
linear algebra routines they need as NLOPT_BLAS(name)(...) with the
Fortran calling convention, either from an external BLAS/LAPACK or from
built-in versions.
New CMake option NLOPT_BLAS (case-insensitive): AUTO (default) uses a
BLAS/LAPACK found by FindLAPACK if it provides all the routines in
nlopt-blas.h, ON requires one, OFF never uses one, and any other value
is passed as the library list (names, paths or flags) to link against.
The Fortran symbol suffix ("_64_", "_64", "_" or none) is found by
test-linking all the routines, without running anything, so this works
when cross-compiling; a "64" suffix means 64-bit integers. Libraries
that provide both interfaces get the one matching the pointer size.
The libraries are added to Libs.private in nlopt.pc.
Without a BLAS/LAPACK, nlopt-blas-builtin.h provides plain C ports of
the reference BLAS/LAPACK 3.12 routines (pointer-sized integers),
restricted to the arguments used (e.g. DGELSY without column pivoting
and with the basic instead of the minimum-norm solution when rank
deficient). They are static inline, so that each algorithm can inline
them. Compared with OpenBLAS on random inputs they agree to rounding
and make the same rank decisions.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…C rewrite Replace the f2c translation of Kraft's SLSQP (main iteration, LSQ, LSEI, LSI, LDP, NNLS, LDL' update and the hand-written BLAS-1 loops) with a port of SciPy's C rewrite of the same routines (scipy/scipy#22524), which follows the same iterates but delegates the dense linear algebra to BLAS/LAPACK (dgeqr2/dorm2r, dgerq2/dormr2, dgelsy, dtrsm, dtpmv, ...) through nlopt-blas.h. This cuts the time per iteration by 6-10x for a few hundred variables with an external BLAS/LAPACK. The NLopt wrapper and its stopping criteria are kept; the NLopt changes of the old translation are carried over (bound clamping, backing off from non-finite trial points, no iteration limit, gradient at the first line-search trial point). Unbounded variables are passed as NaN bounds as in SciPy, and the DGELSY call marks all columns as fixed explicitly instead of relying on stale workspace contents (same result as SciPy after the first iteration, and about twice as fast as pivoted QR). New algorithm parameters: linesearch_eager_grad (default 1, current behavior): evaluate the gradient at the first line-search trial point; 0 restores Kraft's function-only trial points. kraft_acc (default 0, current behavior): enables Kraft's own convergence tests with the given accuracy (SciPy's ftol). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Without inequality constraints, LSEI solves its least-squares problem with DGELSY, and without equality constraints the matrix it passes is the upper triangular LSQ matrix built by LSQ. Its QR factorization leaves it unchanged (every reflector is the identity), but LAPACK's blocked DGEQRF still multiplies the zeros: O(n^3) per iteration, 40% of the run time on a dense unconstrained problem with n = 800. lsq_triangular does the remaining steps of DGELSY directly: the same incremental condition estimate (DLAIC1) for the rank and the same DTRSM solve, so the result is bit-identical. DGELSY is still called whenever it would do something else (rescaling of A or b, or rank deficiency for the minimum-norm solution). Triangularity follows from the problem structure (no equality constraints); debug builds assert it. On a dense unconstrained quadratic with n = 800 this cuts the time per evaluation with OpenBLAS from 8.5 ms to 1.7 ms; all evaluated points are bit-identical. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
In each outer iteration NNLS computed the dual vector with one DDOT per column of the zero set and applied the new Householder transformation with one single-column DLARF per remaining column (as in SciPy). With the LDP of a problem with many inequality constraints, that is millions of tiny BLAS calls. Compute the dual vector with one DGEMV over all columns and apply the transformation with DLARF to the two contiguous column blocks on either side of the new column; the columns of the solution set are exactly zero in the rows involved, so they are unchanged. DLARF needs a work array as long as the number of inequality rows, which is added to the buffer size. On a beam-spectrum problem with 200-320 variables and 400-3600 constraints this cuts the time with OpenBLAS by 14% at equal evaluation counts. The summation order changes with an external BLAS, so results can differ from SciPy's in the last bits (same optima); the built-in routines keep the previous loop order and results. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
yuyichao
force-pushed
the
slsqp-lapack
branch
from
October 30, 2026 14:34
c7b2751 to
f72953b
Compare
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.
The implementation itself is completely generated by Claude since I'm not very familiar with this algorithm. I've guided and checked the blas/lapack detection though...
By default this should not increase dependency of nlopt thanks to the bundled blas/lapack implementation but we can get better implementation with a more optimized blas implementation.
During the implementation, claude noticed another difference in the old implementation regarding termination condition which was the actual reason for the slowness compared to the old scipy implementation when too small of a tolerance was specified. This and the line search change was added as solver parameters for the user to choose.
This has been tested extensively by the AI but I'm not able to fully review this.
@ilayn Please have a look at the third commit which implements an optimization for unconstrained problems. It may be applicable for scipy as well. You can use the unconstrained problem in #649 as a benchmark. The forth commit may also be applicable to scipy by batching blas calls.
Fix #649
Close #600