Introduction

This guide provides an overview of the PIPS-IPM++ C++ API for solving large-scale optimization problems with a bordered block-diagonal structure. The API allows users to define a problem structure and provide data using a flexible callback mechanism.

The core of the API revolves around two main components:

A tree structure (DistributedTree) to represent the stages and blocks of the optimization problem.
A set of user-defined callback functions to provide the solver with the problem's dimensions, matrix, and vector data on demand.

This document will walk through the main concepts and provide a step-by-step guide to using the API, with references to concrete examples.

Core Concepts

Problem Structure

PIPS-IPM++ is designed to solve large-scale linear programs (LPs) where the constraint matrix has a bordered block-diagonal structure (BlockLP).

In difference to (BlockLP), the structure of the problem that PIPS-IPM++ expects distinguishes equality and inequality constraints, allows arbitrary variable bounds and more variable types, and supports a quadratic objective function. It looks like this:

\begin{align*}\text{minimize}\; & C + \sum_{k=0}^N c_k^T x_k + x_k^T Q_k x_k \\ \text{subject to}\; & \begin{pmatrix} A_0 & & & & \\ A_1 & B_1 & & & \\ A_2 & & B_2 & & \\ \vdots & & & \ddots & \\ A_N & & & & B_N \\ Bl_0 & Bl_1 & Bl_2 & \cdots & Bl_N \\ \end{pmatrix} \begin{pmatrix} x_0 \\ x_1 \\ x_2 \\ \vdots \\ x_N \end{pmatrix} = \begin{pmatrix} b_0 \\ b_1 \\ b_2 \\ \vdots \\ b_N \\ bl \end{pmatrix} \\ & \begin{pmatrix} c_0^{\text{low}} \\ c_1^{\text{low}} \\ c_2^{\text{low}} \\ \vdots \\ c_N^{\text{low}} \\ dl^{\text{low}} \end{pmatrix} \leq \begin{pmatrix} C_0 & & & & \\ C_1 & D_1 & & & \\ C_2 & & D_2 & & \\ \vdots & & & \ddots & \\ C_N & & & & D_N \\ Dl_0 & Dl_1 & Dl_2 & \cdots & Dl_N \\ \end{pmatrix} \begin{pmatrix} x_0 \\ x_1 \\ x_2 \\ \vdots \\ x_N \end{pmatrix} \leq \begin{pmatrix} c_0^{\text{upp}} \\ c_1^{\text{upp}} \\ c_2^{\text{upp}} \\ \vdots \\ c_N^{\text{upp}} \\ dl^{\text{upp}} \end{pmatrix} \\ & x_k^{\text{low}} \leq x_k \leq x_k^{\text{upp}} \quad \forall k \\ & (x_k)_i \in \mathbb{R} \;\text{or}\; \mathbb{Z} \;\text{or}\; \{0,1\} \quad \forall i\,\forall k \end{align*}

where \(C\) is a constant, \(c_k\), \(b_k\), \(c_k^{\text{low}}\), \(c_k^{\text{upp}}\), \(x_k^{\text{low}}\), \(x_k^{\text{upp}}\), \(k=1,\ldots,N\), \(bl\), \(dl^{\text{low}}\), and \(dl^{\text{upp}}\) are vectors, and \(Q_k\), \(A_k\), \(B_k\), \(Bl_k\), \(C_k\), \(D_k\), and \(Dl_k\) are matrices.

In the constraints:

The last row of blocks \((Bl_0, Bl_1, \ldots, Bl_N)\) and \((Dl_0, Dl_1, \ldots, Dl_N)\) represent linking rows that connect multiple blocks.
The first column of blocks containing \(A_0, \ldots, A_N\), \(Bl_0\), \(C_0, \ldots, C_N\), and \(Dl_0\) correspond to linking columns or "master" variables that are coupled across subproblems.
The diagonal blocks \(B_1, \ldots, B_N\), \(D_1, \ldots, D_N\) represent the constraints for individual, independent subproblems or blocks.

The PIPS-IPM++ API represents this structure as a tree. The root of the tree corresponds to the linking variables/constraints, and each child of the root corresponds to one of the \(N\) blocks.

Note: Even though the API supports a quadratic objective function and non-continuous variables, PIPS-IPM++ can currently only solve continuous linear problems. The \(Q_k\) matrices should therefore always set to be empty and all variables should be declared to be continuous.

The DistributedTree

The abstract class pipsipmpp::DistributedTree is the cornerstone for representing the problem structure. One needs to implement a derived class of DistributedTree that provides the actual problem data. For convenience, PIPS-IPM++ provides pipsipmpp::DistributedInputTree, which takes C-style function pointers as callbacks.

Each node in the tree is an object that can provide information about its part of the problem (local variables, constraints, and connections to its parent and children).

The Callback Mechanism

Instead of building the entire problem matrix in memory, PIPS-IPM++ uses a "callback" mechanism. The solver executes user-provided "callbacks" to request the data it needs at a given time. This is a memory-efficient way to handle very large-scale problems.

The primary types of callbacks are:

Size Callbacks: To get the dimensions of variables and constraints (e.g., number of variables n, number of equality constraints my).
Non-zero Count Callbacks: To get the number of non-zero elements in a given matrix. This allows the solver to pre-allocate the correct amount of memory.
Data Callbacks: To get the actual values of vectors (like the objective c or right-hand-side b) and matrices (like A, B, C, D in CSR format).

The callbacks are C-style functions that must match the PipsFNNZ, PipsFMAT, PipsFVEC, or PipsFSCALAR typedefs. When PIPS-IPM++ executes these callbacks, it passes on the block number \(k\) for which it requests data.

The callbacks are passed to the pipsipmpp::DistributedInputTree::DistributedInputNode constructor and have the following correspondence to the above problem representation:

n: dimension of \(x_k\)
objconst: objective offset \(C\)
fc: objective coefficients \(c_k\)
fQ: matrix \(Q_k\) (should always be empty)
fnnzQ: number of nonzeros in matrix \(Q_k\) (should always be 0)
fxlow: lower bounds \(x_k^{\text{low}}\)
fixlow: indicator which entries of \(x_k^{\text{low}}\) are active (i.e., not \(-\infty\))
fxupp: upper bounds \(x_k^{\text{upp}}\)
fixupp: indicator which entries of \(x_k^{\text{upp}}\) are active (i.e., not \(+\infty\))
fixtyp: variable types of \(x_k\) (should always be continuous)
my: number of rows in matrices \(A_k\) and \(B_k\)
fA: matrix \(A_k\)
fnnzA: number of nonzeros in matrix \(A_k\)
fB: matrix \(B_k\)
fnnzB: number of nonzeros in matrix \(B_k\)
fb: right-hand-side vector \(b_k\)
myl: number of rows in matrix \(Bl_k\)
fBl: matrix \(Bl_k\)
fnnzBl: number of nonzeros in matrix \(Bl_k\)
fbl: right-hand-side vector \(bl\)
mz: number of rows in matrices \(C_k\) and \(D_k\)
fC: matrix \(C_k\)
fnnzC: number of nonzeros in matrix \(C_k\)
fD: matrix \(D_k\)
fnnzD: number of nonzeros in matrix \(D_k\)
fclow: left-hand-side vector \(c^{\text{low}}_k\)
ficlow: indicator which entries of \(c^{\text{low}}_k\) are active (i.e., not \(-\infty\))
fcupp: right-hand-side vector \(c_k^{\text{upp}}\)
ficupp: indicator which entries of \(c_k^{\text{upp}}\) are active (i.e., not \(+\infty\))
mzl: number of rows in matrix \(Dl_k\)
fDl: matrix \(Dl_k\)
fnnzDl: number of nonzeros in matrix \(Dl_k\)
fdlow: left-hand-side vector \(dl^{\text{low}}\)
fidlow: indicator which entries of \(dl^{\text{low}}\) are active (i.e., not \(-\infty\))
fdupp: right-hand-side vector \(dl^{\text{upp}}\)
fidupp: indicator which entries of \(dl^{\text{upp}}\) are active (i.e., not \(+\infty\))

The callbacks for dimensions (n, my, mz, myl, mzl) must always be provided. For all other callbacks, it is possible to pass on nullptr, to indicate that the corresponding scalar, vector, or matrix is zero.

Step-by-Step Guide

Here is a summary of the steps to define and solve a problem using the C++ API.

Step 1: Define the Callback Functions

First, one needs to write the callbacks that provide the problem data.

A good example is Drivers/CallbackExample/callbackExample.cpp. Let's look at a few key callbacks from that file.

nSize: Provides the number of variables for a given node.
int nSize(void*, int id, int* nnz) {

if (id == 2)

*nnz = 4; // Node 2 has 4 variables

else

*nnz = 2; // Nodes 0 and 1 have 2 variables

return 0;

}
nnzMatEqStage1: Provides the number of non-zeros in the A matrix (first-stage equality constraints).
int nnzMatEqStage1(void*, int id, int* nnz) {

if (id == 0) // Root node

*nnz = 2;

else // Block nodes

*nnz = 2;

return 0;

}
vecEqRhs: Fills a vector with the right-hand-side values for the equality constraints.
int vecEqRhs(void*, int id, double* vec, int) {

if (id == 0) {

vec[0] = 2.0;

vec[1] = 7.0;

}

// ... other blocks

return 0;

}
matEqStage1: Provides the A matrix data in compressed row storage (CRS) format.
int matEqStage1(void*, int id, int* krowM, int* jcolM, double* M) {

if (id == 0) {

M[0] = 2.0; M[1] = 7.0;

krowM[0] = 0; krowM[1] = 1; krowM[2] = 2;

jcolM[0] = 0; jcolM[1] = 1;

}

// ... other blocks

return 0;

}

Step 2: Build the Problem Tree

Once the callbacks have been defined, a pipsipmpp::DistributedInputTree that represents the problem structure needs to be constructed.

Create the root node, passing pointers to all callback functions. The id for the root node is typically 0.
Create child nodes for each block, also passing the callback functions.
Add the child nodes to the root.

Here's how it's done in Drivers/CallbackExample/callbackExample.cpp:

// In main()
using namespace pipsipmpp;
 
ProbData probData(nBlocks);
 
//build the problem tree
std::unique_ptr<DistributedInputTree::DistributedInputNode> data_root = std::make_unique<DistributedInputTree::DistributedInputNode>(&probData, 0, ...);
 
auto* root = new DistributedInputTree(std::move(data_root));
 
for (int id = 1; id <= nBlocks; id++) {
   std::unique_ptr<DistributedInputTree::DistributedInputNode> data_child = std::make_unique<DistributedInputTree::DistributedInputNode>(&probData, id, ...);
   root->add_child(std::make_unique<DistributedInputTree>(std::move(data_child)));
}

The first argument to the pipsipmpp::DistributedInputTree::DistributedInputNode constructor is a void* for user data, which is then passed to every callback. This is useful for passing problem parameters or other data.

Step 3: Create and Run the Solver

After building the tree, an instance of pipsipmpp::PIPSIPMppInterface is created, passing it the tree and the MPI communicator. In addition, parameters can be set by using methods of the pipsipmpp::options namespace. Then, the run() method is called.

options::set_parameter("PRESOLVE", false);
options::set_parameter("SCALER", "geometricmean");
 
PIPSIPMppInterface pipsIpm(root, MPI_COMM_WORLD);
 
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
if (rank == 0)
   std::cout << "solving..." << std::endl;
 
TerminationStatus status = pipsIpm.run();

Since the same code runs in every MPI process, a check has been added here to make sure that solving... is printed only in the first process.

The pipsipmpp::TerminationStatus indicates whether the problem could be solved to optimality, unboundedness or infeasibility was detected, or the solve stopped due to a working limit or numerical problems.

Step 4: Retrieve Results

After run() completes, one can query the solver for the results.

// print optimal objective value
const double objective = pipsIpm.getObjective();
if (rank == 0)
   std::cout << "solving finished... objective value: " << objective << std::endl;
 
// retrieve primal solution (variable values and activity of constraints)
std::vector<double> primalSolVec = pipsIpm.gatherPrimalSolution();
std::vector<double> eqValues = pipsIpm.gatherEqualityConsValues();
std::vector<double> ineqValues = pipsIpm.gatherInequalityConsValues();
 
// retieve dual solution (corresponding to variable bounds and constraints)
std::vector<double> dualSolVarBounds = pipsIpm.gatherDualSolutionVarBounds();
std::vector<double> dualSolEqVec = pipsIpm.gatherDualSolutionEq();
std::vector<double> dualSolIneqVec = pipsIpm.gatherDualSolutionIneq();

The PIPSIPMppInterface Class

The pipsipmpp::PIPSIPMppInterface class is the primary user-facing class for interacting with the PIPS-IPM++ solver. It encapsulates the entire solution process, including reading the problem structure from the pipsipmpp::DistributedTree, presolving, solving with the core interior-point algorithm, and postsolving to return the solution in terms of the original problem.

Key Public Methods

Below are some of the most important methods available in pipsipmpp::PIPSIPMppInterface.

Setup and Execution

PIPSIPMppInterface(DistributedInputTree* tree, MPI_Comm comm, const std::string& settings): The constructor takes the problem representation (tree), an MPI communicator, and an optional path to a settings file (PIPSIPMpp.opt by default).
TerminationStatus run(): This is the main method to start the solution process. It executes the full pipeline: presolve, solve, and postsolve. It returns a pipsipmpp::TerminationStatus enum indicating whether the solution was successful, or if the problem was found to be infeasible, unbounded, or if another issue occurred. Options PRESOLVE, STOP_AFTER_PRESOLVE, and POSTSOLVE can be set to do only presolve, skip presolve, or skip postsolve.

Retrieving Solver Status and Results

TerminationStatus termination_status() const: Returns the final status of the solver, which is the same value returned by run().
int n_iterations() const: Returns the number of interior-point iterations taken by the solver.
double getObjective(): Returns the objective function value for the computed solution. If postsolve has been run, this is the objective of the original problem.

Retrieving Solution Vectors

PIPSIPMppInterface provides a comprehensive set of methods to gather different parts of the solution from all MPI processes to the root process (rank 0). These methods typically return a std::vector<double>.

pipsipmpp::PIPSIPMppInterface::gatherPrimalSolution : Returns the full primal solution vector \(x\).
pipsipmpp::PIPSIPMppInterface::gatherDualSolutionEq : Returns the dual variables for the equality constraints.
pipsipmpp::PIPSIPMppInterface::gatherDualSolutionIneq : Returns the dual variables for the inequality constraints.
pipsipmpp::PIPSIPMppInterface::gatherDualSolutionVarBounds : Returns the dual variables for the variable bounds.
pipsipmpp::PIPSIPMppInterface::gatherEqualityConsValues : Returns the values of the equality constraints.
pipsipmpp::PIPSIPMppInterface::gatherInequalityConsValues : Returns the values of the inequality constraints.

There are many more gather... methods to retrieve specific components of the primal and dual solutions, as well as slacks and residuals. Refer to pipsipmpp::PIPSIPMppInterface for a complete list.

Examples

Basic Example: callbackExample.cpp

The file Drivers/CallbackExample/callbackExample.cpp provides a minimal, self-contained example demonstrating the C-style callback API for a simple two-stage problem. This is the best place to start to understand the basic mechanics. The problem data is hard-coded directly inside the callback functions.

Advanced Example: gmspips_reader.cpp

The file Drivers/gams/gmspips/gmspips_reader.cpp shows a more complex, real-world example of how to interface an external modeling system, in this case GAMS (the General Algebraic Modeling System), with PIPS-IPM++.

In this setup, the optimization problem is not defined in C++ code but is instead read from GDX (GAMS Data Exchange) files. GDX files are binary files that store GAMS model data (parameters, variables, equations, etc.). This approach allows modelers to use the high-level GAMS language to define their large-scale LPs with bordered block diagonal structure, with PIPS-IPM++ acting as the solver engine. See Running on GDX input for more information on this workflow.

The gmspips_reader serves as a bridge, reading the GDX files and translating the data into a format that PIPS-IPM++ can understand, via the callback interface.

Key Mechanisms in gmspips_reader.cpp

Instead of hard-coding data, the gmspips_reader loads problem data on-demand. Here’s a more detailed look at how it works:

Data Loading with GMSPIPSIO: The reader uses a dedicated helper library (GMSPIPSIO) which is responsible for the low-level details of opening GDX files and reading the structured data. The data for each block is loaded into a C-style struct called GMSPIPSBlockData_t. This struct holds all the vectors and matrices for that specific block.
Callback Generation with Macros: To avoid writing repetitive code for each callback, the implementation uses a set of C macros (nCB, nnzCB, vecCB, matCB) to generate the required functions. For example, nCB(ni) generates a function called fsizeni that returns the number of variables, while nnzCB(A) generates fnonzeroA which returns the number of non-zeros in matrix A. This makes the code more concise and less error-prone.
Lazy Loading: A clever checkAndAlloc macro is used to implement lazy loading. The data for a particular block is only read from its GDX file the first time a callback is invoked for that block's ID. This is a crucial optimization when dealing with a very large number of blocks, as data is only loaded into memory when the solver actually needs it for computation.

Example Workflow

Let's trace how the data for a matrix is provided to the solver:

// 1. The solver calls the 'fmatA' callback to get matrix A for a specific 'id'.
// This function was generated by the matCB macro.
int fmatA(void* user_data, int id, int* krowM, int* jcolM, double* M)
{
    // 2. The user_data is cast to the array of block data pointers.
    GMSPIPSBlockData_t** blocks = (GMSPIPSBlockData_t**) user_data;
 
    // 3. 'checkAndAlloc' is called. If blocks[id] is null, it reads the
    //    corresponding GDX file and populates the struct.
    checkAndAlloc(id);
 
    // 4. At this point, the data is guaranteed to be in memory.
    GMSPIPSBlockData_t* blk = blocks[id];
 
    // 5. The matrix data from the 'blk' struct is copied into the
    //    pointers provided by the solver.
    for( int i = 0; i <= blk->mA; i++ ) {
       krowM[i] = blk->rmA[i];
    }
    for( int k = 0; k < blk->nnzA; k++ ) {
       jcolM[k] = blk->ciA[k];
       M[k] = blk->valA[k];
    }
    return 0;
}

This example illustrates how the callback API decouples the solver from the data source. The solver does not need to know where the data comes from; it only needs the callbacks to provide it when requested. This makes the API highly flexible for integrating with different modeling languages and file formats.