PIPS-IPM++ Solver and Tools
a parallel interior-point method for doubly bordered block diagonal linear programs
PIPS-IPM++ C++ API

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 you 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. This structure is common in many areas, including stochastic programming, decomposition methods, and multi-period planning problems.

The general structure of the constraint matrix solved by PIPS-IPM++ looks like this:

\[\begin{pmatrix} A & B_1 & B_2 & \cdots & B_N \\ C_1 & D_1 & & & \\ C_2 & & D_2 & & \\ \vdots & & & \ddots & \\ C_N & & & & D_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 \end{pmatrix} \]

In this matrix:

  • The first block row (A, B_1, ..., B_N) represents linking rows that connect multiple blocks.
  • The first block column containing A, C_1, ..., C_N corresponds to linking columns or "master" variables that are coupled across subproblems.
  • The diagonal matrices D_1, ..., D_N represent the constraints for individual, independent subproblems or blocks.

A primary application that gives rise to this structure is multi-period planning problems, commonly arising in energy systems applications. In this context:

  • The linking variables \(x_0\) represent global strategic decisions that are imposed over many time periods.
  • Each block \(i=1, \dots, N\) corresponds to a time-period or stage.
  • The matrices \(C_i\) (couple the global decisions to the time-period decisions) and \(B_i\) (couple the time-period decision) define the linkage between stages.
  • The matrices \(D_i\) define the constraints within each time-period.

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.

The DistributedTree

The abstract class DistributedTree is the cornerstone for representing the problem structure. You need to implement a derived class of DistributedTree that provides the actual problem data. For convenience, PIPS-IPM++ provides 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 in CSR format).

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, you need to write a set of C-style functions that provide the problem data. These functions must match the FNNZ, FMAT, and FVEC typedefs.

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;
    }

You must provide a complete set of callbacks for all the required parts of your problem (objective, constraints, bounds, etc.).

Step 2: Build the Problem Tree

Once you have defined the callbacks, you need to construct a DistributedInputTree that represents your problem structure.

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

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

// In main()
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)));
}
DistributedInputTree
Represents a node in the distributed input tree, holding problem data or callback functions to retrie...
Definition DistributedInputTree.h:79

The first argument to the 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, you create an instance of PIPSIPMppInterface, passing it the tree and the MPI communicator. Then, simply call the run() method.

PIPSIPMppInterface pipsIpm(root, MPI_COMM_WORLD);
if (rank == 0)
std::cout << "solving...
";
pipsIpm.run();
PIPSIPMppInterface
The main interface class for PIPS-IPM++.
Definition PIPSIPMppInterface.hpp:37

Step 4: Retrieve Results

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

const double objective = pipsIpm.getObjective();
if (rank == 0)
std::cout << "solving finished ... objective value: " << objective << "
";
// You can also get solution vectors
// auto primalSol = pipsIpm.gatherPrimalSolution();
// auto dualSolEq = pipsIpm.gatherDualSolutionEq();

The PIPSIPMppInterface Class

The 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 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 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 TerminationStatus enum indicating whether the solution was successful, or if the problem was found to be infeasible, unbounded, or if another issue occurred.
  • TerminationStatus presolve() This method runs only the presolving and scaling steps. This can be useful for debugging or analyzing the presolved problem. If the STOP_AFTER_PRESOLVE option is set to true, run() will also stop after presolving.
  • void postsolveComputedSolution() If a presolver was used, this method can be called to manually trigger post-solving, which maps the solution of the presolved problem back to the original problem space. run() calls this automatically if necessary.

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>.

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

There are many more gather... methods to retrieve specific components of the primal and dual solutions, as well as slacks and residuals. Refer to 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.

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:

  1. 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.
  2. 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.
  3. 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.