diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/ISQP/ISQP.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/ISQP/ISQP.C
index 9e0eb0d667..f6ab3de548 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/ISQP/ISQP.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/ISQP/ISQP.C
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2021 PCOpt/NTUA
-    Copyright (C) 2020-2021 FOSS GP
+    Copyright (C) 2020-2023 PCOpt/NTUA
+    Copyright (C) 2020-2023 FOSS GP
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,7 +34,7 @@ License
 
 namespace Foam
 {
-    defineTypeNameAndDebug(ISQP, 0);
+    defineTypeNameAndDebug(ISQP, 2);
     addToRunTimeSelectionTable
     (
         updateMethod,
@@ -50,6 +50,15 @@ namespace Foam
 }
 
 
+const Foam::Enum<Foam::ISQP::preconditioners>
+Foam::ISQP::preconditionerNames
+({
+    { preconditioners::diag, "diag" },
+    { preconditioners::invHessian, "invHessian" },
+    { preconditioners::ShermanMorrison, "ShermanMorrison" }
+});
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::ISQP::updateSizes()
@@ -208,6 +217,7 @@ void Foam::ISQP::zeroUpdates()
 
 void Foam::ISQP::solveDeltaPEqn()
 {
+    addProfiling(ISQP, "ISQP::solveDeltaPEqn");
     // Explicit part of the right hand side of the deltaX equation
     scalarField FDx(-resFL());
     if (includeBoundConstraints_)
@@ -233,29 +243,6 @@ void Foam::ISQP::solveDeltaPEqn()
         }
     }
     CGforDeltaP(FDx);
-    /*
-  //Info<< "FDx::" << FDx << endl;
-  //Info<< "invHFL::" << invHFL() << endl;
-
-    // Loop to obtain deltaX. Part of the LHS is treated explicitly, to
-    // avoid solving a (potentially) dense system with the size of the design
-    // variables
-    scalarField rhs(computeRHSForDeltaX(FDx));
-    scalar res(sqrt(globalSum(magSqr(deltaP_ - rhs))));
-    scalar resInit(res);
-    label iter(0);
-    do
-    {
-        deltaP_ = rhs;
-        rhs = computeRHSForDeltaX(FDx);
-        res = sqrt(globalSum(magSqr(deltaP_ - rhs)));
-        DebugInfo
-            << "Solving for deltaX, Initial Residual " << resInit
-            << ", Final Residual " << res << endl;
-        resInit = res;
-    } while (iter++ < maxDxIters_ && res > 1.e-07);
-//  Info<< "deltaX solution " << deltaP_ << endl;
-//  */
 }
 
 
@@ -302,9 +289,10 @@ Foam::tmp<Foam::scalarField> Foam::ISQP::computeRHSForDeltaX
 
 void Foam::ISQP::CGforDeltaP(const scalarField& FDx)
 {
-    scalarField precond(deltaPDiagPreconditioner());
-    scalarField r(FDx - DeltaPMatrixVectorProduct(deltaP_));
-    scalarField z(precond*r);
+    addProfiling(ISQP, "ISQP::CGforDeltaP");
+    autoPtr<scalarField> precon(nullptr);
+    scalarField r(FDx - matrixVectorProduct(deltaP_));
+    scalarField z(preconVectorProduct(r, precon));
     scalarField p(z);
     scalar res(sqrt(globalSum(r*r)));
     scalar resInit(res);
@@ -313,17 +301,17 @@ void Foam::ISQP::CGforDeltaP(const scalarField& FDx)
     label iter(0);
     do
     {
-        scalarField Ap(DeltaPMatrixVectorProduct(p));
+        scalarField Ap(matrixVectorProduct(p));
         scalar a = rz/globalSum(p*Ap);
         deltaP_ += a*p;
         r -= a*Ap;
         res = sqrt(globalSum(r*r));
-        z = precond*r;
+        z = preconVectorProduct(r, precon);
         rz = globalSum(r*z);
         scalar beta = rz/rzOld;
         p = z + beta*p;
         rzOld = rz;
-    } while (iter++ < maxDxIters_ && res > 1.e-09);
+    } while (iter++ < maxDxIters_ && res > relTol_*resInit);
     DebugInfo
         << "CG, Solving for deltaP, Initial Residual " << resInit
         << ", Final Residual " << res
@@ -331,42 +319,13 @@ void Foam::ISQP::CGforDeltaP(const scalarField& FDx)
 }
 
 
-Foam::tmp<Foam::scalarField> Foam::ISQP::deltaPDiagPreconditioner()
-{
-    tmp<scalarField> tpreconditioner(HessianDiag());
-  //tmp<scalarField> tpreconditioner(SR1HessianDiag());
-    scalarField& preconditioner = tpreconditioner.ref();
-
-    // Part related to the constraints
-    forAll(constraintDerivatives_, cI)
-    {
-        scalarField cDerivs(constraintDerivatives_[cI], activeDesignVars_);
-        scalar mult(gs_[cI]/lamdas_[cI]);
-        if (includeExtraVars_)
-        {
-            mult += extraVars_()[cI]/z_()[cI];
-        }
-        preconditioner += sqr(cDerivs)/mult;
-    }
-
-    if (includeBoundConstraints_)
-    {
-        preconditioner += lTilda_()/ls_() + uTilda_()/us_();
-    }
-
-    preconditioner = 1./preconditioner;
-
-    return tpreconditioner;
-}
-
-
-Foam::tmp<Foam::scalarField> Foam::ISQP::DeltaPMatrixVectorProduct
+Foam::tmp<Foam::scalarField> Foam::ISQP::matrixVectorProduct
 (
     const scalarField& vector
 )
 {
+    addProfiling(ISQP, "ISQP::MatrixVectorProduct");
     tmp<scalarField> tAp(HessianVectorProduct(vector));
-  //tmp<scalarField> tAp(SR1HessianVectorProduct(vector));
     scalarField& Ap = tAp.ref();
     scalarField GsLAv(cValues_.size(), Zero);
     forAll(constraintDerivatives_, cI)
@@ -403,6 +362,134 @@ Foam::tmp<Foam::scalarField> Foam::ISQP::DeltaPMatrixVectorProduct
 }
 
 
+Foam::tmp<Foam::scalarField> Foam::ISQP::preconVectorProduct
+(
+    const scalarField& vector,
+    autoPtr<scalarField>& precon
+)
+{
+    addProfiling(ISQP, "ISQP::preconVectorProduct");
+    if (preconType_ == preconditioners::diag)
+    {
+        if (!precon)
+        {
+            precon.reset(diagPreconditioner().ptr());
+        }
+        return precon()*vector;
+    }
+    else if (preconType_ == preconditioners::invHessian)
+    {
+        return invHessianVectorProduct(vector);
+    }
+    else if (preconType_ == preconditioners::ShermanMorrison)
+    {
+        return ShermanMorrisonPrecon(vector);
+    }
+    return nullptr;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::ISQP::diagPreconditioner()
+{
+    addProfiling(ISQP, "ISQP::deltaPDiagPreconditioner");
+    tmp<scalarField> tpreconditioner(HessianDiag());
+    scalarField& preconditioner = tpreconditioner.ref();
+
+    // Part related to the constraints
+    forAll(constraintDerivatives_, cI)
+    {
+        scalarField cDerivs(constraintDerivatives_[cI], activeDesignVars_);
+        scalar mult(gs_[cI]/lamdas_[cI]);
+        if (includeExtraVars_)
+        {
+            mult += extraVars_()[cI]/z_()[cI];
+        }
+        preconditioner += sqr(cDerivs)/mult;
+    }
+
+    if (includeBoundConstraints_)
+    {
+        preconditioner += lTilda_()/ls_() + uTilda_()/us_();
+    }
+
+    preconditioner = 1./preconditioner;
+
+    return tpreconditioner;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::ISQP::ShermanMorrisonPrecon
+(
+    const scalarField& vector
+)
+{
+    // Recursively update the inv(LHS)*vector since the LHS consists of the
+    // L-BFGS-based Hessian, computed with rank-2 updates, and the part related
+    // to flow constraints, computed as rank-1 updates. In the inversion
+    // process, the diagonal matrix related to bound constraints is treated as
+    // the initial matrix of the L-BFGS update.
+
+    // Constribution from bound constraints, treated as the seed of the
+    // L-BFGS inverse
+    tmp<scalarField> tdiag(nullptr);
+    if (includeBoundConstraints_)
+    {
+        tdiag.reset(lTilda_()/ls_() + uTilda_()/us_());
+    }
+
+    // List of vectors to be used in the rank-1 updates related to the flow
+    // constraitns
+    PtrList<scalarField> r1Updates(cValues_.size());
+
+    forAll(constraintDerivatives_, cI)
+    {
+        const scalarField& cDerivsI = constraintDerivatives_[cI];
+        r1Updates.set(cI, new scalarField(cDerivsI, activeDesignVars_));
+    }
+    // Multipliers of the rank-1 updates
+    scalarField mult(gs_/lamdas_);
+    if (includeExtraVars_)
+    {
+        mult += extraVars_()/z_();
+    }
+
+    return
+        ShermanMorrisonRank1Update(r1Updates, vector, tdiag, mult, mult.size());
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::ISQP::ShermanMorrisonRank1Update
+(
+    const PtrList<scalarField>& r1Updates,
+    const scalarField& p,
+    const tmp<scalarField>& diag,
+    const scalarField& mult,
+    label n
+)
+{
+    auto tAp(tmp<scalarField>::New(activeDesignVars_.size(), Zero));
+    scalarField& Ap = tAp.ref();
+
+    if (n == 0)
+    {
+        Ap = invHessianVectorProduct(p, counter_, diag);
+        return tAp;
+    }
+
+    do
+    {
+        --n;
+        Ap = ShermanMorrisonRank1Update(r1Updates, p, diag, mult, n);
+        tmp<scalarField> tAv =
+            ShermanMorrisonRank1Update(r1Updates, r1Updates[n], diag, mult, n);
+        scalarField& Av = tAv.ref();
+        scalar yHs = globalSum(r1Updates[n]*Av)/mult[n];
+        Ap -= Av*globalSum(r1Updates[n]*Ap)/(1 + yHs)/mult[n];
+    } while (n > 0);
+    return tAp;
+}
+
+
 void Foam::ISQP::computeNewtonDirection()
 {
     // Zero the updates computed in the previous optimisation cycle
@@ -410,6 +497,7 @@ void Foam::ISQP::computeNewtonDirection()
 
     // Solve equation for deltaP_. The expensive part of the step. Everything
     // else can be computed based on this
+    addProfiling(ISQP, "ISQP::computeNewtonDirection");
     solveDeltaPEqn();
 
     // deltaLamda
@@ -698,7 +786,10 @@ Foam::tmp<Foam::scalarField> Foam::ISQP::invHFL()
 
 Foam::tmp<Foam::scalarField> Foam::ISQP::resFGs()
 {
-    tmp<scalarField> tFGs(tmp<scalarField>::New(gs_ + cValues_));
+    tmp<scalarField> tFGs
+    (
+        tmp<scalarField>::New(gs_ + cValues_ - max((1 - cRed_)*cValues_, Zero))
+    );
     scalarField& FGs = tFGs.ref();
 
     forAll(constraintDerivatives_, cI)
@@ -795,6 +886,7 @@ Foam::tmp<Foam::scalarField> Foam::ISQP::resFz()
 
 void Foam::ISQP::solveSubproblem()
 {
+    addProfiling(ISQP, "ISQP::solveSubproblem");
     zeroUpdates();
     if (includeBoundConstraints_ || !cValues_.empty())
     {
@@ -825,13 +917,14 @@ void Foam::ISQP::solveSubproblem()
         (
             iter++ < maxNewtonIters_ && (eps_ > epsMin_ || resMax > 0.9*eps_)
         );
+        Info<< "Finished solving the QP subproblem" << nl << endl;
         if (iter == maxNewtonIters_)
         {
             WarningInFunction
                 << "Iterative solution of the QP problem did not converge"
                 << endl;
         }
-        if (debug)
+        if (debug > 2)
         {
             scalarField vars(designVars_().getVars(), activeDesignVars_);
             scalarField newVars(vars + p_);
@@ -961,8 +1054,21 @@ Foam::ISQP::ISQP
         coeffsDict(type).getOrDefault<label>("maxLineSearchIters", 10)
     ),
     maxDxIters_(coeffsDict(type).getOrDefault<label>("maxDpIters", 1000)),
+    relTol_(coeffsDict(type).getOrDefault<scalar>("relTol", 0.01)),
+    preconType_
+    (
+        preconditionerNames.getOrDefault
+        (
+            "preconditioner", coeffsDict(type), preconditioners::diag
+        )
+    ),
+    cRed_
+        (coeffsDict(type).getOrDefault<scalar>("targetConstraintReduction", 1)),
     meritFunctionFile_(nullptr)
 {
+    Info<< "Preconditioner type of the SQP subproblem is ::"
+        << preconditionerNames.names()[preconType_]
+        << endl;
     // Always apply damping of s in ISQP
     useYDamping_ = true;
     useSDamping_ = false;
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/ISQP/ISQP.H b/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/ISQP/ISQP.H
index 56336c41cb..d94f444ded 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/ISQP/ISQP.H
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/ISQP/ISQP.H
@@ -34,14 +34,17 @@ Description
     solved using the interior point method (hence the I in the classe name,
     which is not standard in the terminology used in the literature). The
     (potentially dense) linear system formulated by the interior point method
-    is solved using Conjugate Gradient with a diagonal preconditioner, using
-    matrix-vector products to avoid storing the LHS matrix.
+    is solved using Conjugate Gradient and a choice from three preconditioners
+    (diagonal, inverse Hessian, Sherman-Morrison), using matrix-vector products
+    to avoid storing the LHS matrix.
 
     Bound constraints on the design variables will also be included, if set by
     the designVariables known by the updateMethod. If the QP problem is
     infeasible, the algorithm can still be used by setting includeExtraVars_
     to true, to allow a computation of an update, despite not being able to
-    satisfy all the constraints of the QP problem.
+    satisfy all the constraints of the QP problem.  Alternatively,
+    targetConstraintReduction can be set to a number lower than 1 to help with
+    the feasibility of the dual problem.
 
 
 SourceFiles
@@ -69,6 +72,18 @@ class ISQP
     public LBFGS,
     public SQPBase
 {
+public:
+
+        //- Enumeration of preconditioners for the subproblem
+        enum preconditioners
+        {
+            diag, invHessian, ShermanMorrison
+        };
+
+        //- Names of preconditioners for the subproblem
+        static const Enum<preconditioners> preconditionerNames;
+
+
 protected:
 
     // Protected data
@@ -92,10 +107,6 @@ protected:
 
         // Lagrange multipliers and slack variables
 
-            //- Lagrange multipliers for the inequality constraints
-            //  Inheritated from SQPBase
-            //scalarField lamdas_;
-
             //- Slack variables for the inequality constraints
             scalarField gs_;
 
@@ -154,6 +165,16 @@ protected:
         //- Maxmimum number of iterations for the deltaX equation
         label maxDxIters_;
 
+        //- Relative tolerance of the CG solver, solving for deltaP_
+        scalar relTol_;
+
+        //- Which preconditioner to use for the solution of the subproblem
+        label preconType_;
+
+        //- Percentage reduction for the constraints.
+        //  Can be used to relax QP problems with no feasible points
+        scalar cRed_;
+
         //- File including the l1 merit function
         autoPtr<OFstream> meritFunctionFile_;
 
@@ -193,20 +214,45 @@ protected:
             void solveDeltaPEqn();
 
             //- Compute the RHS for the deltaX equation
+            //  Currently not in use
             tmp<scalarField> computeRHSForDeltaX(const scalarField& FDx);
 
             //- CG algorithm for the solution of the deltaP eqn
             void CGforDeltaP(const scalarField& FDx);
 
-            //- Diagonal preconditioner of the deltaP eqn
-            tmp<scalarField> deltaPDiagPreconditioner();
-
             //- Procudt of the LHS of the deltaP eqn with a vector
-            tmp<scalarField> DeltaPMatrixVectorProduct
+            tmp<scalarField> matrixVectorProduct
             (
                 const scalarField& vector
             );
 
+            //- Preconditioner-vector product for CG
+            //  In case a diagonal preconditioner is used, it is stored in
+            //  precon for all CG iterations
+            tmp<scalarField> preconVectorProduct
+            (
+                const scalarField& vector,
+                autoPtr<scalarField>& precon
+            );
+
+            //- Diagonal preconditioner of the deltaP eqn
+            tmp<scalarField> diagPreconditioner();
+
+            //- Use the Sherman-Morrison formula to compute the matrix-free
+            //- product of the approximate inverse of the LHS with a vector
+            tmp<scalarField> ShermanMorrisonPrecon(const scalarField& vector);
+
+            //- Recursive (naive) implementation of the rank-1 update
+            //  WIP, efficient for up to 2 flow-related constraints
+            tmp<scalarField> ShermanMorrisonRank1Update
+            (
+                const PtrList<scalarField>& r1Updates,
+                const scalarField& p,
+                const tmp<scalarField>& diag,
+                const scalarField& mult,
+                label n
+            );
+
             //- Perform line search and return max residual corresponding to
             //- the updated solution
             scalar lineSearch();
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/LBFGS/LBFGS.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/LBFGS/LBFGS.C
index 6d82c97299..3617635957 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/LBFGS/LBFGS.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/LBFGS/LBFGS.C
@@ -174,29 +174,36 @@ Foam::tmp<Foam::scalarField>
 Foam::LBFGS::invHessianVectorProduct
 (
     const scalarField& vector,
-    const label counter
+    const label counter,
+    tmp<scalarField> diag
 )
 {
     // Sanity checks
     tmp<scalarField> tq(tmp<scalarField>::New(activeDesignVars_.size(), Zero));
     scalarField& q = tq.ref();
-    if (vector.size() == designVars_().getVars().size())
+    label nv = designVars_().getVars().size();
+    label nav = activeDesignVars_.size();
+    if (vector.size() == nv)
     {
         q.map(vector, activeDesignVars_);
     }
-    else if (vector.size() == activeDesignVars_.size())
+    else if (vector.size() == nav)
     {
         q = vector;
     }
     else
     {
         FatalErrorInFunction
-            << "Size of input vector is equal to neither the number of "
-            << " design variabes nor that of the active design variables"
+            << "Size of input vector "
+            << "(" << vector.size() << ") "
+            << "is equal to neither the number of design variabes "
+            << "(" << nv << ")"
+            << " nor that of the active design variables "
+            << "(" << nav << ")"
             << exit(FatalError);
     }
 
-    if (counter != 0)
+    if (counter)
     {
         // L-BFGS two loop recursion
         //~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -206,16 +213,21 @@ Foam::LBFGS::invHessianVectorProduct
         scalarField r(nSteps, 0.);
         for (label i = nLast; i > -1; --i)
         {
-          //Info << "Y " << y_[i] << endl;
-          //Info << "S " << s_[i] << endl;
             r[i] = 1./globalSum(y_[i]*s_[i]);
             a[i] = r[i]*globalSum(s_[i]*q);
             q -= a[i]*y_[i];
         }
 
         scalar gamma =
-            globalSum(y_[nLast]*s_[nLast])/globalSum(y_[nLast]*y_[nLast]);
-        q *= gamma;
+            globalSum(y_[nLast]*y_[nLast])/globalSum(y_[nLast]*s_[nLast]);
+        if (diag)
+        {
+            q /= (gamma + diag());
+        }
+        else
+        {
+            q /= gamma;
+        }
 
         scalarField b(activeDesignVars_.size(), Zero);
         for (label i = 0; i < nSteps; ++i)
@@ -224,6 +236,10 @@ Foam::LBFGS::invHessianVectorProduct
             q += s_[i]*(a[i] - b);
         }
     }
+    else if (diag)
+    {
+        q /= (1 + diag());
+    }
 
     return tq;
 }
@@ -243,6 +259,7 @@ Foam::LBFGS::HessianVectorProduct
     const label counter
 )
 {
+    addProfiling(LBFGS, "LBFGS::HessianVectorProduct");
     // Sanity checks
     tmp<scalarField> tq(tmp<scalarField>::New(activeDesignVars_.size(), Zero));
     scalarField& q = tq.ref();
@@ -273,7 +290,7 @@ Foam::LBFGS::HessianVectorProduct
 
         // Product of the last matrix on the right with the input vector
         scalarField SKsource(2*nSteps, Zero);
-        for(label i = 0; i < nSteps; ++i)
+        for (label i = 0; i < nSteps; ++i)
         {
             SKsource[i] = delta*globalSum(s_[i]*source);
             SKsource[i + nSteps] = globalSum(y_[i]*source);
@@ -374,11 +391,11 @@ Foam::tmp<Foam::scalarField> Foam::LBFGS::HessianDiag()
 
         // Product of the inverse of the middle matrix with the right vector
         List<scalarField> MR(2*nSteps, scalarField(n, Zero));
-        for(label k = 0; k < n; ++k)
+        for (label k = 0; k < n; ++k)
         {
-            for(label i = 0; i < 2*nSteps; ++i)
+            for (label i = 0; i < 2*nSteps; ++i)
             {
-                for(label j = 0; j < nSteps; ++j)
+                for (label j = 0; j < nSteps; ++j)
                 {
                     MR[i][k] +=
                         invM[i][j]*delta*s_[j][k]
@@ -390,9 +407,9 @@ Foam::tmp<Foam::scalarField> Foam::LBFGS::HessianDiag()
         // Part of the Hessian diagonal computed by the multiplication
         // of the above matrix with the left matrix of the recursive Hessian
         // reconstruction
-        for(label k = 0; k < n; ++k)
+        for (label k = 0; k < n; ++k)
         {
-            for(label j = 0; j < nSteps; ++j)
+            for (label j = 0; j < nSteps; ++j)
             {
                 diag[k] -=
                     delta*s_[j][k]*MR[j][k] + y_[j][k]*MR[j + nSteps][k];
@@ -448,7 +465,7 @@ Foam::LBFGS::SR1HessianVectorProduct
 
         // Product of the last matrix on the right with the input vector
         scalarField YBSsource(nSteps, Zero);
-        for(label i = 0; i < nSteps; ++i)
+        for (label i = 0; i < nSteps; ++i)
         {
             YBSsource[i] = globalSum((y_[i] - delta*s_[i])*source);
         }
@@ -544,11 +561,11 @@ Foam::tmp<Foam::scalarField> Foam::LBFGS::SR1HessianDiag()
 
         // Product of the inverse of the middle matrix with the right vector
         List<scalarField> MR(nSteps, scalarField(n, Zero));
-        for(label k = 0; k < n; ++k)
+        for (label k = 0; k < n; ++k)
         {
-            for(label i = 0; i < nSteps; ++i)
+            for (label i = 0; i < nSteps; ++i)
             {
-                for(label j = 0; j < nSteps; ++j)
+                for (label j = 0; j < nSteps; ++j)
                 {
                     MR[i][k] += invM[i][j]*(y_[j][k] - delta*s_[j][k]);
                 }
@@ -558,9 +575,9 @@ Foam::tmp<Foam::scalarField> Foam::LBFGS::SR1HessianDiag()
         // Part of the Hessian diagonal computed by the multiplication
         // of the above matrix with the left matrix of the recursive Hessian
         // reconstruction
-        for(label k = 0; k < n; ++k)
+        for (label k = 0; k < n; ++k)
         {
-            for(label j = 0; j < nSteps; ++j)
+            for (label j = 0; j < nSteps; ++j)
             {
                 diag[k] += (y_[j][k] - delta*s_[j][k])*MR[j][k];
             }
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/LBFGS/LBFGS.H b/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/LBFGS/LBFGS.H
index 2f6db81888..e785e7997a 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/LBFGS/LBFGS.H
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/LBFGS/LBFGS.H
@@ -117,9 +117,10 @@ protected:
         //- Update design variables
         virtual void update();
 
-        //- Compute the inverse Hessian - vector product
+        //- Compute the inverse Hessian - vector product.
         //  Input should have the size of all design variables or the active
-        //  ones, output is the size of the active design variables
+        //  ones, output is the size of the active design variables.
+        //  Can optionally received a seed diagonal matrix
         virtual tmp<scalarField> invHessianVectorProduct
         (
             const scalarField& vector
@@ -129,7 +130,8 @@ protected:
         tmp<scalarField> invHessianVectorProduct
         (
             const scalarField& vector,
-            const label counter
+            const label counter,
+            tmp<scalarField> diag = nullptr
         );
 
         //- Compute the Hessian - vector product
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/updateMethod/updateMethod.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/updateMethod/updateMethod.C
index 84a7ad4f7f..4178195c9b 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/updateMethod/updateMethod.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/updateMethod/updateMethod.C
@@ -174,16 +174,11 @@ Foam::updateMethod::inv(SquareMatrix<scalar> A)
 
 Foam::scalar Foam::updateMethod::globalSum(const scalarField& field)
 {
-    scalar value(0);
     if (globalSum_)
     {
-        value = gSum(field);
+        return gSum(field);
     }
-    else
-    {
-        value = sum(field);
-    }
-    return value;
+    return sum(field);
 }
 
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/parameterization/NURBS/NURBS3DSurface/NURBS3DSurface.C b/src/optimisation/adjointOptimisation/adjoint/parameterization/NURBS/NURBS3DSurface/NURBS3DSurface.C
index 5e602bd3ab..c0cbf085bc 100644
--- a/src/optimisation/adjointOptimisation/adjoint/parameterization/NURBS/NURBS3DSurface/NURBS3DSurface.C
+++ b/src/optimisation/adjointOptimisation/adjoint/parameterization/NURBS/NURBS3DSurface/NURBS3DSurface.C
@@ -1479,7 +1479,7 @@ scalar NURBS3DSurface::lengthU
     }
 
     // Integrate.
-    for(label uI = 0; uI<(uLenSize - 1); uI++)
+    for (label uI = 0; uI<(uLenSize - 1); uI++)
     {
         const label ptI((uIStart+uI)*nVPts_ + vIConst);
 
@@ -1512,7 +1512,7 @@ scalar NURBS3DSurface::lengthU
     }
 
     // Integrate.
-    for(label uI = 0; uI<(nPts - 1); uI++)
+    for (label uI = 0; uI<(nPts - 1); uI++)
     {
         uLength +=
             0.5*(mag(dxdu[uI + 1]) + mag(dxdu[uI]))*(localU[uI + 1]-localU[uI]);
@@ -1556,7 +1556,7 @@ scalar NURBS3DSurface::lengthV
     }
 
     // Integrate.
-    for(label vI = 0; vI<(vLenSize - 1); vI++)
+    for (label vI = 0; vI<(vLenSize - 1); vI++)
     {
         const label ptI((uIConst)*nVPts_ + (vIStart + vI));
 
@@ -1589,7 +1589,7 @@ scalar NURBS3DSurface::lengthV
     }
 
     // Integrate.
-    for(label vI = 0; vI<(nPts - 1); vI++)
+    for (label vI = 0; vI<(nPts - 1); vI++)
     {
         vLength +=
             0.5*(mag(dxdv[vI + 1]) + mag(dxdv[vI]))*(localV[vI + 1]-localV[vI]);
@@ -1979,7 +1979,7 @@ scalar NURBS3DSurface::lengthDerivativeU
     }
 
     // Integrate.
-    for(label uI=0; uI<(nPts-1); uI++)
+    for (label uI=0; uI<(nPts-1); uI++)
     {
         ulDerivative +=
             0.5
@@ -2017,7 +2017,7 @@ scalar NURBS3DSurface::lengthDerivativeV
     }
 
     // Integrate.
-    for(label vI=0; vI<(nPts-1); vI++)
+    for (label vI=0; vI<(nPts-1); vI++)
     {
         vlDerivative +=
             0.5
@@ -2046,9 +2046,9 @@ const labelList& NURBS3DSurface::getBoundaryCPIDs()
 
         // v-constant cps
         label bID(0);
-        for(label vI=0; vI<vNCPs; vI+=vNCPs-1)
+        for (label vI=0; vI<vNCPs; vI+=vNCPs-1)
         {
-            for(label uI=0; uI<uNCPs; uI++)
+            for (label uI=0; uI<uNCPs; uI++)
             {
                 const label CPI(vI*uNCPs + uI);
                 whichBoundaryCPID_()[CPI] = bID;
@@ -2056,10 +2056,10 @@ const labelList& NURBS3DSurface::getBoundaryCPIDs()
             }
         }
         // u-constant cps
-        for(label uI=0; uI<uNCPs; uI+=uNCPs-1)
+        for (label uI=0; uI<uNCPs; uI+=uNCPs-1)
         {
             // corner CPS already accounted for
-            for(label vI=1; vI<vNCPs-1; vI++)
+            for (label vI=1; vI<vNCPs-1; vI++)
             {
                 const label CPI(vI*uNCPs + uI);
                 whichBoundaryCPID_()[CPI] = bID;