Working on control theory problems, I’ve recently needed to be able to solve convex quadratic programs (QP) of the form
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) 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 (still work-in-progress). The legacy solver has a pretty verbose interface that reads like so
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
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
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
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
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).
. If you share the same vision, then that’s awesome!!