# Modernizing \`quadprog\` for convex quadratic programming **URL:** \ **Category:** Uncategorized\ **Created:** [March 15, 2025, 1:22pm UTC](https://fortran-lang.discourse.group/t/modernizing-quadprog-for-convex-quadratic-programming/9355 "2025-03-15T13:22:14Z")\ **Posts on this page:** 15\ **Page:** 1
**Author:** ![loiseaujc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/loiseaujc/32/4582_2.png) [@loiseaujc](https://fortran-lang.discourse.group/u/loiseaujc)\ **Post date:** [March 15, 2025, 1:22pm UTC](https://fortran-lang.discourse.group/t/modernizing-quadprog-for-convex-quadratic-programming/9355/1 "2025-03-15T13:22:14Z")
Working on control theory problems, I’ve recently needed to be able to solve convex quadratic programs (QP) of the form \begin{aligned} \textrm{minimize} & \quad \dfrac12 \mathbf{x}^T \mathbf{Px} - \mathbf{x}^T \mathbf{q} \\ \textrm{subject to} & \quad \mathbf{Ax} = \mathbf{b} \\ & \quad \mathbf{Cx} \geq \mathbf{d}, \end{aligned} where \mathbf{P} is a dense symmetric positive definite matrix of moderate size (say 512 x 512). The matrices \mathbf{A} and \mathbf{C} are also dense. One of the best algorithms for such small-to-moderate scale QP is the Godfarb & Idnani algorithm. A reference implementation is the one by Berwin Turlach [(here)](https://cran.r-project.org/web/packages/quadprog/) which is used under the hood by `quadprog` in R and other languages. The original implementation being written in F77, I’ve converted it to F90 and made very small changes here and there to the code to make it completely stand-alone (it does not even need `blas`, `linpack` or `lapack`). You can find the code on this [Github repo](https://github.com/loiseaujc/QuadProg) (still work-in-progress). The legacy solver has a pretty verbose interface that reads like so ```fortran interface pure module subroutine qpgen2(dmat, dvec, fddmat, n, sol, lagr, crval, amat, bvec, fdamat, q, & meq, iact, nact, iter, work, ierr) integer, intent(in) :: fddmat, n !! Dimensions of the symmetric positive definite matrix Dmat. integer, intent(in) :: fdamat, q !! Dimensions of the constraint matrix Amat integer, intent(in) :: meq !! Number of equality constraints. integer, intent(out) :: iact(*), nact !! Indices and number of active constraints at the optimum. integer, intent(out) :: iter(*) !! Number of iterations. integer, intent(inout) :: ierr !! Information flag. real(dp), intent(inout) :: dmat(fddmat, *), dvec(*) !! Sym. pos. def. matrix and vector defining the quadratic cost. real(dp), intent(out) :: lagr(*), sol(*) !! Lagrange multipliers and solution vector. real(dp), intent(inout) :: amat(fdamat, *), bvec(*) !! Matrix and vector defining the (in-)equality constraints. real(dp), intent(inout) :: work(*) !! Workspace. real(dp), intent(out) :: crval !! Cost function at the optimum. end subroutine end interface ``` In the process of modernizing the code and making it more user-friendly, I wrapped this `qpgen2` subroutine into another function with the following interface ````fortran interface !! #### Description !! !! ... !! !! #### Syntax !! !! - To solve an unconstrained quadratic program: !! !! ```fortran !! call solve_qp(P, q, x) !! ``` !! !! #### Arguments !! module function solve_qp(P, q, Aeq, beq, C, d, y, obj, & factorized, overwrite_p, info) result(x) real(dp), intent(inout), contiguous, target :: P(:, :) !! n x n Symmetric positive-definite matrix defining the quadratic form. real(dp), intent(in) :: q(:) !! Vector defining the linear term of the quadratic cost. real(dp), optional, intent(in) :: Aeq(:, :) !! Matrix defining the linear equality constraints Aeq @ x = beq. real(dp), optional, intent(in) :: beq(:) !! Right-hand side vector of the linear equality constraints Aeq @ x = beq. real(dp), optional, intent(in) :: C(:, :) !! Matrix defining the linear inequality constraints C @ x >= d. real(dp), optional, intent(in) :: d(:) !! Right-hand side vector of the linear inequality constraints C @ x >= d. real(dp), optional, allocatable, target, intent(out) :: y(:) !! Vector of Lagrange multipliers at the optimum. real(dp), optional, intent(out), target :: obj !! Cost function at the optimum. logical, optional, intent(inout) :: factorized !! Whether P has already been factorized using ??? (default .false.) logical, optional, intent(in) :: overwrite_p !! Whether P can be overwritten (default .false.) integer, intent(out) :: info !! Information flag. real(dp), allocatable :: x(:) !! Vector corresponding to the minimizer. end function end interface ```` which is much more flexible on the user-side and also quite a bit more natural. In my set of unit tests, everything works as expected if I use directly `qpgen2`, no matter whether I use `gfortran` or `ifx` to compile the code with `fpm`. However, if I use `solve_qp`, this is where the sh\*t hits the fan and I may need help. Solving an unconstrained QP with the following call ```fortran x = solve_qp(P, q, obj=obj, info=info) ``` causes not problem whatsoever, no matter the compiler. However, if I want to solve a constrained problem, and using the following calling sequence ```fortran x = solve_qp(P=P, q=q, Aeq=C, beq=b, y=y, obj=obj, info=info) ``` only `ifx` can compile the code and produce the expected output. If I use `gfortran`, it encounters an internal compile error ```auto 33%] Compiling... test/problems.f90:158:67: 158 | x = solve_qp(P=P, q=q, Aeq=C, beq=b, y=y, obj=obj, info=info) | 1 internal compiler error: Segmentation fault 0x7fc99922551f ??? ./signal/../sysdeps/unix/sysv/linux/x86_64/libc_sigaction.c:0 0x7fc99920cd8f __libc_start_call_main ../sysdeps/nptl/libc_start_call_main.h:58 0x7fc99920ce3f __libc_start_main_impl ../csu/libc-start.c:392 Please submit a full bug report, with preprocessed source (by using -freport-bug). Please include the complete backtrace with any bug report. ``` I’ve never really had to deal with internal compiler errors so far and so I was wondering if anyone could give me some insights on how to generate a proper bug report for `gfortran` given that I ain’t sure how to reproduce this issue with a MWE (since I don’t know exactly where it comes from). ---
**Author:** ![hkvzjal](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/hkvzjal/32/3055_2.png) [@hkvzjal](https://fortran-lang.discourse.group/u/hkvzjal)\ **Post date:** [March 15, 2025, 1:45pm UTC](https://fortran-lang.discourse.group/t/modernizing-quadprog-for-convex-quadratic-programming/9355/2 "2025-03-15T13:45:09Z")
I suspect the problem coming from the fact that `info` is a non-optional argument but you placed it at the end of the procedure signature. My advise, try again placing all non-optional first and then all optional arguments. Non-related to your post: `forall` is a deprecated feature of the language, I know it is neat for one-liners, but I would advise against its inclusion in new code. ---
**Author:** ![loiseaujc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/loiseaujc/32/4582_2.png) [@loiseaujc](https://fortran-lang.discourse.group/u/loiseaujc)\ **Post date:** [March 15, 2025, 1:49pm UTC](https://fortran-lang.discourse.group/t/modernizing-quadprog-for-convex-quadratic-programming/9355/3 "2025-03-15T13:49:44Z")
Nope, didn’t change anything. That was my first idea as well. As for `forall`, yep, it was just out of convenience. I eventually plan to make the code as much F2018 compliant as possible with the exception of some `go to` in `qpgen2` which I don’t want to mess with. ---
**Author:** ![hkvzjal](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/hkvzjal/32/3055_2.png) [@hkvzjal](https://fortran-lang.discourse.group/u/hkvzjal)\ **Post date:** [March 15, 2025, 2:17pm UTC](https://fortran-lang.discourse.group/t/modernizing-quadprog-for-convex-quadratic-programming/9355/4 "2025-03-15T14:17:06Z")
Humm, I tweaked it a bit here [Compiler Explorer](https://godbolt.org/z/KPdKdhdE6) it seems to be coming from your choice of splitting the procedure signature within the interface declaration and having a `module procedure` for the actual implementation. Moving your whole signature below, the ICE goes away. ---
**Author:** ![loiseaujc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/loiseaujc/32/4582_2.png) [@loiseaujc](https://fortran-lang.discourse.group/u/loiseaujc)\ **Post date:** [March 15, 2025, 2:24pm UTC](https://fortran-lang.discourse.group/t/modernizing-quadprog-for-convex-quadratic-programming/9355/5 "2025-03-15T14:24:34Z")
Awesome ! I’ll make the changes on Monday as I get back to the lab’. I’ll also add some derived types to have an even simpler interface. On a different note, at some point, I think we, as a community, would need to start a discussion about the future of `stdlib`, what is beyond its scope, etc. Basically, if we think of `stdlib` as the Fortran equivalent of `numpy`, what about `scipy` then? Given this topic, I’m obviously thinking about where an optimization module would fit. ---
**Author:** ![hkvzjal](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/hkvzjal/32/3055_2.png) [@hkvzjal](https://fortran-lang.discourse.group/u/hkvzjal)\ **Post date:** [March 15, 2025, 2:31pm UTC](https://fortran-lang.discourse.group/t/modernizing-quadprog-for-convex-quadratic-programming/9355/6 "2025-03-15T14:31:23Z")
> [@loiseaujc](#): > > On a different note, at some point, I think we, as a community, would need to start a discussion about the future of `stdlib`, what is beyond its scope, etc. Basically, if we think of `stdlib` as the Fortran equivalent of `numpy`, what about `scipy` then. Given this topic, I’m obviously thinking about where an optimization module would fit. Yep, I’m personally convinced that `stdlib` ought to be a complete scientific toolbox for the community. It should cover all those nice tools and algorithms one can find in `numpy`, `scipy`, `scikit` and more. I wanted to suggest, if you were willing to do this modernization in the scope of stdlib but refrained doing so because it should be your decision 🙂 . If you share the same vision, then that’s awesome!! ---
**Author:** ![fxm](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/fxm/32/5726_2.png) [@fxm](https://fortran-lang.discourse.group/u/fxm)\ **Post date:** [March 15, 2025, 2:50pm UTC](https://fortran-lang.discourse.group/t/modernizing-quadprog-for-convex-quadratic-programming/9355/7 "2025-03-15T14:50:12Z")
I’m all for stdlib becoming a complete scientific toolbox, given that Fortran is so heavily used in science. I’m currently stretched a little thin, but I’ve been meaning to (and still plan to) dive a little deeper into stdlib, to understand how it works, and to see if/how I could contribute. Would love nothing more than to have a single lib as a toolbox for all regular analyses (i.e. everything I don’t need specialised models for). It would also make it less complicated to re-introduce it to university classrooms for a number of reasons. ---
**Author:** ![hkvzjal](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/hkvzjal/32/3055_2.png) [@hkvzjal](https://fortran-lang.discourse.group/u/hkvzjal)\ **Post date:** [March 15, 2025, 2:54pm UTC](https://fortran-lang.discourse.group/t/modernizing-quadprog-for-convex-quadratic-programming/9355/8 "2025-03-15T14:54:02Z")
This warms my heart 😁 ---
**Author:** ![rwmsu](https://avatars.discourse-cdn.com/v4/letter/r/48db29/32.png) [@rwmsu](https://fortran-lang.discourse.group/u/rwmsu)\ **Post date:** [March 15, 2025, 5:37pm UTC](https://fortran-lang.discourse.group/t/modernizing-quadprog-for-convex-quadratic-programming/9355/9 "2025-03-15T17:37:40Z")
> [@loiseaujc](#): > > eventually plan to make the code as much F2018 compliant as possible with the exception of some `go to` in `qpgen2` which I don’t want to mess with. Removing the go to’s is fairly straightforward but takes a little thought. I made the mods and they seem to work with ifx. The check program said everything worked. I’ll PM my updated solve.QP file. ---
**Author:** ![loiseaujc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/loiseaujc/32/4582_2.png) [@loiseaujc](https://fortran-lang.discourse.group/u/loiseaujc)\ **Post date:** [March 16, 2025, 3:20pm UTC](https://fortran-lang.discourse.group/t/modernizing-quadprog-for-convex-quadratic-programming/9355/10 "2025-03-16T15:20:07Z")
> [@hkvzjal](#): > > I wanted to suggest, if you were willing to do this modernization in the scope of stdlib but refrained doing so because it should be your decision 🙂 . If you share the same vision, then that’s awesome!! Definitely ! For the moment I needed to have rapidly a slightly modernized stand-alone version of this algorithm for a project I’m working one but in the long term (e.g. couple of months) that is definitely something I had in mind. One caveat though is the licence of the original code (LGPL if I recall correctly) which might be incompatible with that of `stdlib`. I’ll be in the lab in Paris tomorrow with no students, no meeting, no nothing. I’ll take some time to write a longer post reflecting on the past year as an active `stdlib` user and, from this (admittedly biased) point of view, how I’d like to see the project evolve in the near and longer future. ---
**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\ **Post date:** [March 16, 2025, 6:20pm UTC](https://fortran-lang.discourse.group/t/modernizing-quadprog-for-convex-quadratic-programming/9355/11 "2025-03-16T18:20:48Z")
Yes, stdlib should be the full toolbox, kind of like scipy. As to compiler errors, create an MRE (minimal reproducible example) and submit it as a bug report. ---
**Author:** ![loiseaujc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/loiseaujc/32/4582_2.png) [@loiseaujc](https://fortran-lang.discourse.group/u/loiseaujc)\ **Post date:** [April 7, 2025, 12:19pm UTC](https://fortran-lang.discourse.group/t/modernizing-quadprog-for-convex-quadratic-programming/9355/13 "2025-04-07T12:19:35Z")
Quick update for those interested: I’ve wrapped up the first iteration on the original code. It exposes a couple of new derived-types and interfaces for solving (strictly) convex quadratic programs. Here is a simple example ```fortran program test_program use QuadProg implicit none integer, parameter :: n = 3 ! Size of the problem. type(qp_problem) :: problem type(OptimizeResult) :: result real(dp) :: P(n, n), q(n), C(n, n), d(n) ! Qyadratic cost. P = 0.0_dp; forall (i=1:n) P(i, i) = 1.0_dp q = [0.0_dp, 5.0_dp, 0.0_dp] ! Inequality constraints C(1, :) = [-4, -3, 0] C(2, :) = [2, 1, 0] C(3, :) = [0, -2, 1] d = [-8, 2, 0] problem = qp_problem(P, q, C=C, d=d) result = solve(problem) print *, "Success ?", result%success print *, "Solution :", result%x print *, "Lagrange multipliers :", result%y print *, "Objective function :", result%obj end program ``` You can also find a rudimentary model predictive control (MPC) example using a double integrator with bounded actuation signal in the [example](https://github.com/loiseaujc/QuadProg/tree/main/example/mpc) folder. ---
**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\ **Post date:** [April 10, 2025, 12:14am UTC](https://fortran-lang.discourse.group/t/modernizing-quadprog-for-convex-quadratic-programming/9355/14 "2025-04-10T00:14:57Z")
Thanks for your project. I tried to build it on Windows with gfortran using fpm but got ```auto c:\fortran\public_domain\github\QuadProg>fpm build + mkdir build\dependencies Initialized empty Git repository in C:/fortran/public_domain/github/QuadProg/build/dependencies/test-drive/.git/ From https://github.com/nekStab/test-drive * branch HEAD -> FETCH_HEAD [0%] constants.f90 [14%] constants.f90 done. [14%] aind.f90 [28%] aind.f90 done. [28%] QuadProg.f90 [42%] QuadProg.f90 done. [42%] solvers.f90 [57%] solvers.f90 done. [57%] main.f90 [71%] main.f90 done. [71%] libQuadProg.a [85%] libQuadProg.a done. [85%] mpc.exe [100%] mpc.exe done. C:/equation/bin/../lib/gcc/x86_64-w64-mingw32/15.0.0/../../../../x86_64-w64-mingw32/bin/ld.exe: cannot find -lblas: No such file or directory C:/equation/bin/../lib/gcc/x86_64-w64-mingw32/15.0.0/../../../../x86_64-w64-mingw32/bin/ld.exe: cannot find -llapack: No such file or directory collect2.exe: error: ld returned 1 exit status Compilation failed for object " mpc.exe " stopping due to failed compilation STOP 1 ``` Shortening the above, the key errors are ```auto ld.exe: cannot find -lblas: No such file or directory ld.exe: cannot find -llapack: No such file or directory ``` I was able to compile by [manually linking to librefblas.a](https://github.com/Beliavsky/QuadProg-gfortran-windows) but would like to be able to build with fpm. ---
**Author:** ![loiseaujc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/loiseaujc/32/4582_2.png) [@loiseaujc](https://fortran-lang.discourse.group/u/loiseaujc)\ **Post date:** [April 10, 2025, 6:28am UTC](https://fortran-lang.discourse.group/t/modernizing-quadprog-for-convex-quadratic-programming/9355/15 "2025-04-10T06:28:27Z")
I haven’t used Windows in 20 years or so. My understanding of `fpm` was that, provided you have `blas` and `lapack` already available on your system, having ```toml [build] link = ['blas', 'lapack'] ``` in the `fpm.toml` file should be enough for `fpm` to compile the code correctly. It certainly does when using `unbuntu` or `macOS`. On a similar note: I haven’t been able to use `Windows` in my CI with Github Actions. Here is my `ci.yml` file: ```auto name: ci on: [push, pull_request] jobs: test: runs-on: ${{ matrix.os }} strategy: fail-fast: false matrix: include: - os: ubuntu-latest toolchain: {compiler: gcc, version: 14, flags: ['-cpp -O3 -march=native -mtune=native']} # - os: windows-latest # toolchain: {compiler: gcc, version: 13, flags: ['-cpp -O3 -march=native -mtune=native']} - os: macos-13 toolchain: {compiler: gcc, version: 13, flags: ['-cpp -O3 -march=native -mtune=native']} steps: - name: Checkout code uses: actions/checkout@v1 - uses: fortran-lang/setup-fortran@main id: setup-fortran with: compiler: ${{ matrix.toolchain.compiler }} version: ${{ matrix.toolchain.version }} - name: Setup Fortran Package Manager uses: fortran-lang/setup-fpm@v7 with: github-token: ${{ secrets.GITHUB_TOKEN }} - name: Install BLAS and LAPACK if: contains( matrix.os, 'ubuntu') run: | sudo apt-get install libblas-dev liblapack-dev - run: | fpm test --flag "${{ join(matrix.toolchain.flags, ' ') }}" ``` It works perfectly fine for `unbuntu`, installing correctly `blas` and `lapack`. Whenever I uncomment the `Windows` runner I get error about missing `blas` and `lapack`. Any one knows how I should modify this yml file to get them installed when using the `Windows` runner? ---
**Author:** ![HugoMVale](https://avatars.discourse-cdn.com/v4/letter/h/e19adc/32.png) [@HugoMVale](https://fortran-lang.discourse.group/u/HugoMVale)\ **Post date:** [April 10, 2025, 7:45am UTC](https://fortran-lang.discourse.group/t/modernizing-quadprog-for-convex-quadratic-programming/9355/16 "2025-04-10T07:45:01Z")
It can be a bit tricky to get BLAs running across all OSs… you may have a look at my odrpack95 and odrpack-python projects… > **[GitHub - HugoMVale/odrpack95: Software for Weighted Orthogonal Distance...](https://github.com/HugoMVale/odrpack95)** > > Software for Weighted Orthogonal Distance Regression with Bound Constraints > **[GitHub - HugoMVale/odrpack-python: Python bindings for the modernized version of...](https://github.com/HugoMVale/odrpack-python)** > > Python bindings for the modernized version of odrpack95.