diff --git a/CMakeLists.txt b/CMakeLists.txt
index 716c45c2..7f5b3677 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -142,7 +142,7 @@ find_library(BQPD bqpd)
 if(NOT BQPD)
    message(WARNING "Optional library BQPD was not found.")
 else()
-   list(APPEND UNO_SOURCE_FILES uno/solvers/BQPD/BQPDSolver.cpp uno/solvers/BQPD/wdotd.f)
+   list(APPEND UNO_SOURCE_FILES uno/solvers/BQPD/BQPDSolver.cpp)
    list(APPEND TESTS_UNO_SOURCE_FILES unotest/functional_tests/BQPDSolverTests.cpp)
    link_to_uno(bqpd ${BQPD})
 endif()
diff --git a/uno/solvers/BQPD/BQPDSolver.cpp b/uno/solvers/BQPD/BQPDSolver.cpp
index 678780ea..c3ab4deb 100644
--- a/uno/solvers/BQPD/BQPDSolver.cpp
+++ b/uno/solvers/BQPD/BQPDSolver.cpp
@@ -14,14 +14,14 @@
 #include "options/Options.hpp"
 #include "fortran_interface.h"
 
-#define WSC FC_GLOBAL(wsc,WSC)
-#define KKTALPHAC FC_GLOBAL(kktalphac,KKTALPHAC)
-#define BQPD FC_GLOBAL(bqpd,BQPD)
+#define WSC FC_GLOBAL(wsc, WSC)
+#define ALPHAC FC_GLOBAL(alphac, ALPHAC)
+#define BQPD FC_GLOBAL(bqpd, BQPD)
+#define hessian_vector_product FC_GLOBAL(gdotx, GDOTX)
 
-namespace uno {
-#define BIG 1e30
+extern "C" {
+   void hessian_vector_product(int *n, const double x[], const double ws[], const int lws[], double v[]);
 
-   extern "C" {
    // fortran common block used in bqpd/bqpd.f
    extern struct {
       int kk, ll, kkk, lll, mxws, mxlws;
@@ -30,13 +30,16 @@ namespace uno {
    // fortran common for inertia correction in wdotd
    extern struct {
       double alpha;
-   } KKTALPHAC;
+   } ALPHAC;
 
    extern void
    BQPD(const int* n, const int* m, int* k, int* kmax, double* a, int* la, double* x, double* bl, double* bu, double* f, double* fmin, double* g,
          double* r, double* w, double* e, int* ls, double* alp, int* lp, int* mlp, int* peq, double* ws, int* lws, const int* mode, int* ifail,
          int* info, int* iprint, int* nout);
-   }
+}
+
+namespace uno {
+   #define BIG 1e30
 
    // preallocate a bunch of stuff
    BQPDSolver::BQPDSolver(size_t number_variables, size_t number_constraints, size_t number_objective_gradient_nonzeros, size_t number_jacobian_nonzeros,
@@ -104,7 +107,7 @@ namespace uno {
       WSC.ll = static_cast<int>(this->size_hessian_sparsity);
       WSC.mxws = static_cast<int>(this->size_hessian_workspace);
       WSC.mxlws = static_cast<int>(this->size_hessian_sparsity_workspace);
-      KKTALPHAC.alpha = 0; // inertia control
+      ALPHAC.alpha = 0; // inertia control
 
       if (this->print_subproblem) {
          DEBUG << "objective gradient: " << linear_objective;
@@ -311,4 +314,24 @@ namespace uno {
       }
       throw std::invalid_argument("The BQPD ifail is not consistent with the Uno status values");
    }
-} // namespace
\ No newline at end of file
+} // namespace
+
+void hessian_vector_product(int *n, const double x[], const double ws[], const int lws[], double v[]) {
+   assert(n != nullptr && "BQPDSolver::hessian_vector_product: the dimension n passed by pointer is NULL");
+
+   for (size_t i = 0; i < static_cast<size_t>(*n); i++) {
+      v[i] = 0.;
+   }
+
+   int footer_start = lws[0];
+   for (int i = 0; i < *n; i++) {
+      for (int k = lws[footer_start + i]; k < lws[footer_start + i + 1]; k++) {
+         int j = lws[k] - 1;
+         v[i] += ws[k - 1] * x[j];
+         if (j != i) {
+            // off-diagonal term
+            v[j] += ws[k - 1] * x[i];
+         }
+      }
+   }
+}
\ No newline at end of file
diff --git a/uno/solvers/BQPD/wdotd.f b/uno/solvers/BQPD/wdotd.f
deleted file mode 100644
index e7ab5954..00000000
--- a/uno/solvers/BQPD/wdotd.f
+++ /dev/null
@@ -1,118 +0,0 @@
-C Copyright (c) 2018-2024 Sven Leyffer
-C Licensed under the MIT license. See LICENSE file in the project directory for details.
-
-C     hristen this file wdotd.f
-
-      subroutine wdotd (n, x, ws, lws, result)
-
-c     ==========================================================
-c     Computes result = W.x where W is Hessian and x is a vector for AMPL
-c     Assumes v=0 on entry (OK, if called from gdotx, see QPsolve*.f)
-c     ==========================================================
-
-      implicit none
-
-c     ... declaration of passed parameters
-      integer          n, lws(0:*)
-      double precision x(n), result(n), ws(*)
-
-c     ... declaration of internal variables
-      integer i, j, k, footer_start
-
-c     inertia control for diagonal terms
-      double precision alpha
-      common /kktalphac/ alpha
-
-c     ========================  procedure body  =========================
-
-c     ... form result = W.x from sparse, upper triangular Hessian
-      footer_start = lws(0)
-      do i=1,n
-         do k=lws(footer_start + i - 1), lws(footer_start + i)-1
-            j = lws(k)
-            result(i) = result(i) + ws(k)*x(j)
-            if (j.ne.i) then
-c               off-diagonal term
-               result(j) = result(j) + ws(k)*x(i)
-            else
-c               diagonal term
-               result(i) = result(i) + alpha*x(j)
-            endif
-         enddo
-      enddo
-      return
-      end
-
-c     ******************************************************************
-
-      subroutine gdotx (n, x, ws, lws, result)
-
-      implicit none
-
-c     ... declaration of passed parameters
-      integer n, lws(*)
-      double precision    x(n), result(n), ws(*)
-
-c     ... declaration of internal variables
-      integer i
-
-c     ... storage map for hessian and scale_mode
-      integer         scale_mode, phe
-      common /scalec/ scale_mode, phe
-
-c     ========================  procedure body  =========================
-
-c     ... set result = 0
-      do i=1,n
-         result(i) = 0.D0
-      enddo
-
-c     ... allow for scaling of variables 
-      if ((scale_mode.eq.1).or.(scale_mode.eq.3)) then
-         do i=1,n
-            x(i) = x(i) * ws(i)
-         enddo
-      endif
-
-c     ... form v = W.d from sparse, upper triangular Hessian
-      call Wdotd (n, x, ws(phe+1), lws, result)
-
-c     ... allow for scaling of variables 
-      if ((scale_mode.eq.1).or.(scale_mode.eq.3)) then
-         do i=1,n
-            result(i) = result(i) * ws(i)
-            x(i) = x(i) / ws(i)
-         enddo
-      endif
-
-      return
-      end
-c     ******************************************************************
-      subroutine saipy2(s,a,la,i,y,n)
-      implicit double precision (a-h,o-z)
-      dimension a(*),la(0:*),y(*)
-c     ========================  procedure body  =========================
-c     saxpy with column i of A: y + s*A_{i, :}
-      if(s.eq.0.D0) return
-      j_column_start = la(0) + i
-      do j = la(j_column_start), la(j_column_start+1)-1
-         i_variable = la(j)
-         y(i_variable) = y(i_variable) + s*a(j)
-      enddo
-      return
-      end
-
-c     **************************** E N D *********************************
-      function daiscpr2(n,a,la,i,x,b)
-      implicit double precision (a-h,o-z)
-      dimension a(*),la(0:*),x(*)
-      DOUBLE PRECISION daiscpr2
-c     dot product of x and row i of A
-      daiscpr2 = dble(b)
-      j_column_start = la(0) + i
-      do j = la(j_column_start), la(j_column_start+1)-1
-         i_variable = la(j)
-         daiscpr2 = daiscpr2 + dble(x(i_variable))*dble(a(j))
-      enddo
-      return
-      end