diff --git a/sc/Library_solver.mk b/sc/Library_solver.mk
index 8bca81033ffd..cbc2d2bf988d 100644
--- a/sc/Library_solver.mk
+++ b/sc/Library_solver.mk
@@ -52,6 +52,7 @@ $(eval $(call gb_Library_add_exception_objects,solver,\
 	sc/source/solver/SwarmSolver \
 	sc/source/solver/SolverComponent \
 	$(if $(ENABLE_COINMP), sc/source/solver/CoinMPSolver) \
+	$(if $(ENABLE_COINMP), sc/source/solver/CoinMPModel) \
 ))
 
 # vim: set noet sw=4 ts=4:
diff --git a/sc/source/solver/CoinMPModel.cxx b/sc/source/solver/CoinMPModel.cxx
new file mode 100644
index 000000000000..f20ce4adb644
--- /dev/null
+++ b/sc/source/solver/CoinMPModel.cxx
@@ -0,0 +1,207 @@
+/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+/*
+ * This file is part of the LibreOffice project.
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ *
+ */
+
+#include <CoinMP.h>
+#include <CoinError.hpp>
+
+#include "CoinMPModel.hxx"
+
+#include <memory>
+#include <stdexcept>
+
+using namespace com::sun::star;
+
+CoinMpSolveResult coinmpSolveLinearModel(const uno::Sequence<table::CellAddress>& rVariables,
+                                         const uno::Sequence<sheet::SolverConstraint>& rConstraints,
+                                         const table::CellAddress& rObjective,
+                                         const ScSolverCellHashMap& rCells,
+                                         const std::vector<double>& rLowerBounds,
+                                         const std::vector<double>& rUpperBounds,
+                                         const std::vector<char>& rColumnType, bool bMaximize,
+                                         double fTimeoutSeconds)
+{
+    CoinMpSolveResult aResult;
+
+    const size_t nVariables = rVariables.size();
+    size_t nVar = 0;
+
+    // objective: slope of each variable, plus the constant term at the origin
+    const std::vector<double>& rObjectCoefficients = rCells.at(rObjective);
+    std::unique_ptr<double[]> pObjectCoefficients(new double[nVariables]);
+    for (nVar = 0; nVar < nVariables; nVar++)
+        pObjectCoefficients[nVar] = rObjectCoefficients[nVar + 1];
+    double fObjectConstant = rObjectCoefficients[0];
+
+    // rows: collect all coefficients row-wise first
+    size_t nRows = rConstraints.getLength();
+    size_t nDenseSize = nVariables * nRows;
+    std::unique_ptr<double[]> pDenseMatrix(new double[nDenseSize]);
+    for (size_t i = 0; i < nDenseSize; i++)
+        pDenseMatrix[i] = 0.0;
+
+    std::unique_ptr<double[]> pRHS(new double[nRows]);
+    std::unique_ptr<char[]> pRowType(new char[nRows]);
+    for (size_t i = 0; i < nRows; i++)
+    {
+        pRHS[i] = 0.0;
+        pRowType[i] = 'N';
+    }
+
+    for (sal_Int32 nConstraintPosition = 0; nConstraintPosition < rConstraints.getLength();
+         ++nConstraintPosition)
+    {
+        // integer and binary constraints are reflected in the column types
+        sheet::SolverConstraintOperator eOperator = rConstraints[nConstraintPosition].Operator;
+        if (eOperator == sheet::SolverConstraintOperator_LESS_EQUAL
+            || eOperator == sheet::SolverConstraintOperator_GREATER_EQUAL
+            || eOperator == sheet::SolverConstraintOperator_EQUAL)
+        {
+            double fDirectValue = 0.0;
+            bool bRightCell = false;
+            table::CellAddress aRightAddress;
+            const uno::Any& rRightAny = rConstraints[nConstraintPosition].Right;
+            if (rRightAny >>= aRightAddress)
+                bRightCell = true;
+            else
+                rRightAny >>= fDirectValue;
+
+            table::CellAddress aLeftAddress = rConstraints[nConstraintPosition].Left;
+
+            const std::vector<double>& rLeftCoefficients = rCells.at(aLeftAddress);
+            double* pValues = &pDenseMatrix[nConstraintPosition * nVariables];
+            for (nVar = 0; nVar < nVariables; nVar++)
+                pValues[nVar] = rLeftCoefficients[nVar + 1];
+
+            // if the left hand cell has a constant term, move it to the right
+            double fRightValue = -rLeftCoefficients[0];
+
+            if (bRightCell)
+            {
+                const std::vector<double>& rRightCoefficients = rCells.at(aRightAddress);
+                for (nVar = 0; nVar < nVariables; nVar++)
+                    pValues[nVar] -= rRightCoefficients[nVar + 1];
+
+                fRightValue += rRightCoefficients[0];
+            }
+            else
+                fRightValue += fDirectValue;
+
+            switch (eOperator)
+            {
+                case sheet::SolverConstraintOperator_LESS_EQUAL:
+                    pRowType[nConstraintPosition] = 'L';
+                    break;
+                case sheet::SolverConstraintOperator_GREATER_EQUAL:
+                    pRowType[nConstraintPosition] = 'G';
+                    break;
+                case sheet::SolverConstraintOperator_EQUAL:
+                    pRowType[nConstraintPosition] = 'E';
+                    break;
+                default:
+                    break;
+            }
+            pRHS[nConstraintPosition] = fRightValue;
+        }
+    }
+
+    // turn the row-wise matrix into the column-wise form CoinMP wants
+    std::unique_ptr<int[]> pMatrixBegin(new int[nVariables + 1]);
+    std::unique_ptr<int[]> pMatrixCount(new int[nVariables]);
+    std::unique_ptr<double[]> pMatrix(new double[nDenseSize]);
+    std::unique_ptr<int[]> pMatrixIndex(new int[nDenseSize]);
+    int nMatrixPosition = 0;
+    for (nVar = 0; nVar < nVariables; nVar++)
+    {
+        int nBegin = nMatrixPosition;
+        for (size_t nRow = 0; nRow < nRows; nRow++)
+        {
+            double fCoefficient = pDenseMatrix[nRow * nVariables + nVar];
+            if (fCoefficient != 0.0)
+            {
+                pMatrix[nMatrixPosition] = fCoefficient;
+                pMatrixIndex[nMatrixPosition] = nRow;
+                ++nMatrixPosition;
+            }
+        }
+        pMatrixBegin[nVar] = nBegin;
+        pMatrixCount[nVar] = nMatrixPosition - nBegin;
+    }
+    pMatrixBegin[nVariables] = nMatrixPosition;
+    pDenseMatrix.reset();
+
+    std::unique_ptr<double[]> pLowerBounds(new double[nVariables]);
+    std::unique_ptr<double[]> pUpperBounds(new double[nVariables]);
+    std::unique_ptr<char[]> pColumnType(new char[nVariables]);
+    for (nVar = 0; nVar < nVariables; nVar++)
+    {
+        pLowerBounds[nVar] = rLowerBounds[nVar];
+        pUpperBounds[nVar] = rUpperBounds[nVar];
+        pColumnType[nVar] = rColumnType[nVar];
+    }
+
+    int nObjectSense = bMaximize ? SOLV_OBJSENS_MAX : SOLV_OBJSENS_MIN;
+
+    HPROB hProblem = CoinCreateProblem("");
+    int nResult = CoinLoadProblem(hProblem, nVariables, nRows, nMatrixPosition, 0, nObjectSense,
+                                  fObjectConstant, pObjectCoefficients.get(), pLowerBounds.get(),
+                                  pUpperBounds.get(), pRowType.get(), pRHS.get(), nullptr,
+                                  pMatrixBegin.get(), pMatrixCount.get(), pMatrixIndex.get(),
+                                  pMatrix.get(), nullptr, nullptr, nullptr);
+    if (nResult == SOLV_CALL_SUCCESS)
+        nResult = CoinLoadInteger(hProblem, pColumnType.get());
+
+    pColumnType.reset();
+    pMatrixIndex.reset();
+    pMatrix.reset();
+    pMatrixCount.reset();
+    pMatrixBegin.reset();
+    pUpperBounds.reset();
+    pLowerBounds.reset();
+    pRowType.reset();
+    pRHS.reset();
+    pObjectCoefficients.reset();
+
+    CoinSetRealOption(hProblem, COIN_REAL_MAXSECONDS, fTimeoutSeconds);
+    CoinSetRealOption(hProblem, COIN_REAL_MIPMAXSEC, fTimeoutSeconds);
+
+    if (nResult == SOLV_CALL_SUCCESS)
+        nResult = CoinCheckProblem(hProblem);
+
+    if (nResult == SOLV_CALL_SUCCESS)
+    {
+        try
+        {
+            nResult = CoinOptimizeProblem(hProblem, 0);
+        }
+        catch (const CoinError& rError)
+        {
+            CoinUnloadProblem(hProblem);
+            throw std::runtime_error(rError.message());
+        }
+    }
+
+    aResult.bSuccess = (nResult == SOLV_CALL_SUCCESS);
+    if (aResult.bSuccess)
+    {
+        aResult.aSolution.resize(nVariables);
+        CoinGetSolutionValues(hProblem, aResult.aSolution.data(), nullptr, nullptr, nullptr);
+        aResult.fObjective = CoinGetObjectValue(hProblem);
+    }
+    else
+    {
+        aResult.nStatus = CoinGetSolutionStatus(hProblem);
+    }
+
+    CoinUnloadProblem(hProblem);
+
+    return aResult;
+}
+
+/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
diff --git a/sc/source/solver/CoinMPModel.hxx b/sc/source/solver/CoinMPModel.hxx
new file mode 100644
index 000000000000..09928268c767
--- /dev/null
+++ b/sc/source/solver/CoinMPModel.hxx
@@ -0,0 +1,46 @@
+/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+/*
+ * This file is part of the LibreOffice project.
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ *
+ */
+
+#pragma once
+
+#include <com/sun/star/sheet/SolverConstraint.hpp>
+#include <com/sun/star/table/CellAddress.hpp>
+#include <com/sun/star/uno/Sequence.hxx>
+
+#include "SolverComponent.hxx"
+
+#include <vector>
+
+// Result of solving one linear problem with CoinMP.
+struct CoinMpSolveResult
+{
+    bool bSuccess = false;
+    // CoinGetSolutionStatus when the solve did not succeed: 1 infeasible,
+    // 2 unbounded, other values for the remaining cases.
+    int nStatus = 0;
+    std::vector<double> aSolution;
+    double fObjective = 0.0;
+};
+
+// Build a linear problem from already sampled coefficients and solve it with
+// CoinMP. For the objective and for each constraint cell, rCells holds a vector
+// whose first entry is the value of an affine model of that cell at the origin
+// (all variables zero) and whose remaining entries are the slope of each
+// variable. The caller supplies the variable bounds and column types (which
+// encode integer and binary variables), so the same routine serves both a
+// one-shot linear solve and a linearized step of an outer nonlinear loop.
+CoinMpSolveResult coinmpSolveLinearModel(
+    const css::uno::Sequence<css::table::CellAddress>& rVariables,
+    const css::uno::Sequence<css::sheet::SolverConstraint>& rConstraints,
+    const css::table::CellAddress& rObjective, const ScSolverCellHashMap& rCells,
+    const std::vector<double>& rLowerBounds, const std::vector<double>& rUpperBounds,
+    const std::vector<char>& rColumnType, bool bMaximize, double fTimeoutSeconds);
+
+/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
diff --git a/sc/source/solver/CoinMPSolver.cxx b/sc/source/solver/CoinMPSolver.cxx
index 5f8969f197ac..de66301f3f0c 100644
--- a/sc/source/solver/CoinMPSolver.cxx
+++ b/sc/source/solver/CoinMPSolver.cxx
@@ -17,9 +17,7 @@
  *   the License at http://www.apache.org/licenses/LICENSE-2.0 .
  */
 
-#include <CoinMP.h>
-#include <CoinError.hpp>
-
+#include "CoinMPModel.hxx"
 #include "SolverComponent.hxx"
 #include "strings.hrc"
 
@@ -27,7 +25,7 @@
 #include <com/sun/star/table/CellAddress.hpp>
 
 #include <rtl/math.hxx>
-#include <stdexcept>
+#include <algorithm>
 #include <vector>
 #include <float.h>
 
@@ -136,128 +134,12 @@ void SAL_CALL CoinMPSolver::solve()
     if ( !maStatus.isEmpty() )
         return;
 
-    //
-    // build parameter arrays for CoinMP
-    //
-
-    // set objective function
-
-    const std::vector<double>& rObjCoeff = aCellsHash[maObjective];
-    std::unique_ptr<double[]> pObjectCoeffs(new double[nVariables]);
-    for (nVar=0; nVar<nVariables; nVar++)
-        pObjectCoeffs[nVar] = rObjCoeff[nVar+1];
-    double nObjectConst = rObjCoeff[0];             // constant term of objective
-
-    // add rows
-
-    size_t nRows = maConstraints.getLength();
-    size_t nCompSize = nVariables * nRows;
-    std::unique_ptr<double[]> pCompMatrix(new double[nCompSize]);    // first collect all coefficients, row-wise
-    for (size_t i=0; i<nCompSize; i++)
-        pCompMatrix[i] = 0.0;
-
-    std::unique_ptr<double[]> pRHS(new double[nRows]);
-    std::unique_ptr<char[]> pRowType(new char[nRows]);
-    for (size_t i=0; i<nRows; i++)
-    {
-        pRHS[i] = 0.0;
-        pRowType[i] = 'N';
-    }
-
-    for (sal_Int32 nConstrPos = 0; nConstrPos < maConstraints.getLength(); ++nConstrPos)
-    {
-        // integer constraints are set later
-        sheet::SolverConstraintOperator eOp = maConstraints[nConstrPos].Operator;
-        if ( eOp == sheet::SolverConstraintOperator_LESS_EQUAL ||
-             eOp == sheet::SolverConstraintOperator_GREATER_EQUAL ||
-             eOp == sheet::SolverConstraintOperator_EQUAL )
-        {
-            double fDirectValue = 0.0;
-            bool bRightCell = false;
-            table::CellAddress aRightAddr;
-            const uno::Any& rRightAny = maConstraints[nConstrPos].Right;
-            if ( rRightAny >>= aRightAddr )
-                bRightCell = true;                  // cell specified as right-hand side
-            else
-                rRightAny >>= fDirectValue;         // constant value
-
-            table::CellAddress aLeftAddr = maConstraints[nConstrPos].Left;
-
-            const std::vector<double>& rLeftCoeff = aCellsHash[aLeftAddr];
-            double* pValues = &pCompMatrix[nConstrPos * nVariables];
-            for (nVar=0; nVar<nVariables; nVar++)
-                pValues[nVar] = rLeftCoeff[nVar+1];
-
-            // if left hand cell has a constant term, put into rhs value
-            double fRightValue = -rLeftCoeff[0];
-
-            if ( bRightCell )
-            {
-                const std::vector<double>& rRightCoeff = aCellsHash[aRightAddr];
-                // modify pValues with rhs coefficients
-                for (nVar=0; nVar<nVariables; nVar++)
-                    pValues[nVar] -= rRightCoeff[nVar+1];
-
-                fRightValue += rRightCoeff[0];      // constant term
-            }
-            else
-                fRightValue += fDirectValue;
-
-            switch ( eOp )
-            {
-                case sheet::SolverConstraintOperator_LESS_EQUAL:    pRowType[nConstrPos] = 'L'; break;
-                case sheet::SolverConstraintOperator_GREATER_EQUAL: pRowType[nConstrPos] = 'G'; break;
-                case sheet::SolverConstraintOperator_EQUAL:         pRowType[nConstrPos] = 'E'; break;
-                default:
-                    OSL_ENSURE( false, "unexpected enum type" );
-            }
-            pRHS[nConstrPos] = fRightValue;
-        }
-    }
-
-    // Find non-zero coefficients, column-wise
-
-    std::unique_ptr<int[]> pMatrixBegin(new int[nVariables+1]);
-    std::unique_ptr<int[]> pMatrixCount(new int[nVariables]);
-    std::unique_ptr<double[]> pMatrix(new double[nCompSize]);    // not always completely used
-    std::unique_ptr<int[]> pMatrixIndex(new int[nCompSize]);
-    int nMatrixPos = 0;
-    for (nVar=0; nVar<nVariables; nVar++)
-    {
-        int nBegin = nMatrixPos;
-        for (size_t nRow=0; nRow<nRows; nRow++)
-        {
-            double fCoeff = pCompMatrix[ nRow * nVariables + nVar ];    // row-wise
-            if ( fCoeff != 0.0 )
-            {
-                pMatrix[nMatrixPos] = fCoeff;
-                pMatrixIndex[nMatrixPos] = nRow;
-                ++nMatrixPos;
-            }
-        }
-        pMatrixBegin[nVar] = nBegin;
-        pMatrixCount[nVar] = nMatrixPos - nBegin;
-    }
-    pMatrixBegin[nVariables] = nMatrixPos;
-    pCompMatrix.reset();
-
-    // apply settings to all variables
-
-    std::unique_ptr<double[]> pLowerBounds(new double[nVariables]);
-    std::unique_ptr<double[]> pUpperBounds(new double[nVariables]);
-    for (nVar=0; nVar<nVariables; nVar++)
-    {
-        pLowerBounds[nVar] = mbNonNegative ? 0.0 : -DBL_MAX;
-        pUpperBounds[nVar] = DBL_MAX;
-
-        // bounds could possibly be further restricted from single-cell constraints
-    }
-
-    std::unique_ptr<char[]> pColType(new char[nVariables]);
-    for (nVar=0; nVar<nVariables; nVar++)
-        pColType[nVar] = mbInteger ? 'I' : 'C';
-
-    // apply single-var integer constraints
+    // variable bounds and column types. The default range is unbounded unless
+    // non-negativity was requested. Integer and binary single-variable
+    // constraints refine these.
+    std::vector<double> aLowerBounds(nVariables, mbNonNegative ? 0.0 : -DBL_MAX);
+    std::vector<double> aUpperBounds(nVariables, DBL_MAX);
+    std::vector<char> aColumnType(nVariables, mbInteger ? 'I' : 'C');
 
     for (const auto& rConstr : maConstraints)
     {
@@ -266,92 +148,44 @@ void SAL_CALL CoinMPSolver::solve()
              eOp == sheet::SolverConstraintOperator_BINARY )
         {
             table::CellAddress aLeftAddr = rConstr.Left;
-            // find variable index for cell
             for (nVar=0; nVar<nVariables; nVar++)
                 if ( AddressEqual( aVariableCells[nVar], aLeftAddr ) )
                 {
                     if ( eOp == sheet::SolverConstraintOperator_INTEGER )
-                        pColType[nVar] = 'I';
+                        aColumnType[nVar] = 'I';
                     else
                     {
-                        pColType[nVar] = 'B';
-                        pLowerBounds[nVar] = 0.0;
-                        pUpperBounds[nVar] = 1.0;
+                        aColumnType[nVar] = 'B';
+                        aLowerBounds[nVar] = 0.0;
+                        aUpperBounds[nVar] = 1.0;
                     }
                 }
         }
     }
 
-    int nObjectSense = mbMaximize ? SOLV_OBJSENS_MAX : SOLV_OBJSENS_MIN;
-
-    HPROB hProb = CoinCreateProblem("");
-    int nResult = CoinLoadProblem( hProb, nVariables, nRows, nMatrixPos, 0,
-                    nObjectSense, nObjectConst, pObjectCoeffs.get(),
-                    pLowerBounds.get(), pUpperBounds.get(), pRowType.get(), pRHS.get(), nullptr,
-                    pMatrixBegin.get(), pMatrixCount.get(), pMatrixIndex.get(), pMatrix.get(),
-                    nullptr, nullptr, nullptr );
-    if (nResult == SOLV_CALL_SUCCESS)
-    {
-        nResult = CoinLoadInteger( hProb, pColType.get() );
-    }
-
-    pColType.reset();
-    pMatrixIndex.reset();
-    pMatrix.reset();
-    pMatrixCount.reset();
-    pMatrixBegin.reset();
-    pUpperBounds.reset();
-    pLowerBounds.reset();
-    pRowType.reset();
-    pRHS.reset();
-    pObjectCoeffs.reset();
-
-    CoinSetRealOption( hProb, COIN_REAL_MAXSECONDS, mnTimeout );
-    CoinSetRealOption( hProb, COIN_REAL_MIPMAXSEC, mnTimeout );
-
     // TODO: handle (or remove) settings: epsilon, B&B depth
 
-    // solve model
+    CoinMpSolveResult aSolveResult
+        = coinmpSolveLinearModel(maVariables, maConstraints, maObjective, aCellsHash, aLowerBounds,
+                                 aUpperBounds, aColumnType, mbMaximize, mnTimeout);
 
-    if (nResult == SOLV_CALL_SUCCESS)
-    {
-        nResult = CoinCheckProblem( hProb );
-    }
-
-    if (nResult == SOLV_CALL_SUCCESS)
-    {
-        try
-        {
-            nResult = CoinOptimizeProblem( hProb, 0 );
-        }
-        catch (const CoinError& e)
-        {
-            CoinUnloadProblem(hProb);
-            throw std::runtime_error(e.message());
-        }
-    }
-
-    mbSuccess = ( nResult == SOLV_CALL_SUCCESS );
+    mbSuccess = aSolveResult.bSuccess;
     if ( mbSuccess )
     {
-        // get solution
-
         maSolution.realloc( nVariables );
-        CoinGetSolutionValues( hProb, maSolution.getArray(), nullptr, nullptr, nullptr );
-        mfResultValue = CoinGetObjectValue( hProb );
+        std::copy(aSolveResult.aSolution.begin(), aSolveResult.aSolution.end(),
+                  maSolution.getArray());
+        mfResultValue = aSolveResult.fObjective;
     }
     else
     {
-        int nSolutionStatus = CoinGetSolutionStatus( hProb );
-        if ( nSolutionStatus == 1 )
+        if (aSolveResult.nStatus == 1)
             maStatus = SolverComponent::GetResourceString( RID_ERROR_INFEASIBLE );
-        else if ( nSolutionStatus == 2 )
+        else if (aSolveResult.nStatus == 2)
             maStatus = SolverComponent::GetResourceString( RID_ERROR_UNBOUNDED );
         // TODO: detect timeout condition and report as RID_ERROR_TIMEOUT
         // (currently reported as infeasible)
     }
-
-    CoinUnloadProblem( hProb );
 }
 
 extern "C" SAL_DLLPUBLIC_EXPORT css::uno::XInterface *