ipx/src/model.cc

// Copyright (c) 2018 ERGO-Code. See license.txt for license.

#include "model.h"
#include <algorithm>
#include <cassert>
#include <cmath>
#include "utils.h"

namespace ipx {

Int Model::Load(const Control& control, Int num_constr, Int num_var,
                const Int* Ap, const Int* Ai, const double* Ax,
                const double* rhs, const char* constr_type, const double* obj,
                const double* lbuser, const double* ubuser) {
    clear();
    Int errflag = CopyInput(num_constr, num_var, Ap, Ai, Ax, rhs, constr_type,
                            obj, lbuser, ubuser);
    if (errflag)
        return errflag;
    control.Log()
        << "Input\n"
        << Textline("Number of variables:") << num_var_ << '\n'
        << Textline("Number of free variables:") << num_free_var_ << '\n'
        << Textline("Number of constraints:") << num_constr_ << '\n'
        << Textline("Number of equality constraints:") << num_eqconstr_ << '\n'
        << Textline("Number of matrix entries:") << num_entries_ << '\n';
    PrintCoefficientRange(control);
    ScaleModel(control);

    // Make an automatic decision for dualization if not specified by user.
    Int dualize = control.dualize();
    if (dualize < 0)
        dualize = num_constr > 2*num_var;
    if (dualize)
        LoadDual();
    else
        LoadPrimal();

    A_.clear();
    AIt_ = Transpose(AI_);
    assert(AI_.begin(num_cols_ + num_rows_) == AIt_.begin(num_rows_));
    FindDenseColumns();
    norm_c_ = Infnorm(c_);
    norm_bounds_ = Infnorm(b_);
    for (double x : lb_)
        if (std::isfinite(x))
            norm_bounds_ = std::max(norm_bounds_, std::abs(x));
    for (double x : ub_)
        if (std::isfinite(x))
            norm_bounds_ = std::max(norm_bounds_, std::abs(x));
    PrintPreprocessingLog(control);
    return 0;
}

void Model::GetInfo(Info *info) const {
    info->num_var = num_var_;
    info->num_constr = num_constr_;
    info->num_entries = num_entries_;
    info->num_rows_solver = num_rows_;
    info->num_cols_solver = num_cols_ + num_rows_; // including slack columns
    info->num_entries_solver = AI_.entries();
    info->dualized = dualized_;
    info->dense_cols = num_dense_cols();
}

void Model::clear() {
    // clear computational form model
    dualized_ = false;
    num_rows_ = 0;
    num_cols_ = 0;
    num_dense_cols_ = 0;
    nz_dense_ = 0;
    AI_.clear();
    AIt_.clear();
    b_.resize(0);
    c_.resize(0);
    lb_.resize(0);
    ub_.resize(0);
    norm_bounds_ = 0.0;
    norm_c_ = 0.0;

    // clear user model
    num_constr_ = 0;
    num_eqconstr_ = 0;
    num_var_ = 0;
    num_free_var_ = 0;
    num_entries_ = 0;
    boxed_vars_.clear();
    constr_type_.clear();
    norm_obj_ = 0.0;
    norm_rhs_ = 0.0;
    scaled_obj_.resize(0);
    scaled_rhs_.resize(0);
    scaled_lbuser_.resize(0);
    scaled_ubuser_.resize(0);
    A_.clear();

    flipped_vars_.clear();
    colscale_.resize(0);
    rowscale_.resize(0);
}

void Model::PresolveStartingPoint(const double* x_user,
                                  const double* slack_user,
                                  const double* y_user,
                                  const double* z_user,
                                  Vector& x_solver,
                                  Vector& y_solver,
                                  Vector& z_solver) const {
    const Int m = rows();
    const Int n = cols();
    assert(x_solver.size() == n+m);
    assert(y_solver.size() == m);
    assert(z_solver.size() == n+m);

    Vector x_temp(num_var_);
    Vector slack_temp(num_constr_);
    Vector y_temp(num_constr_);
    Vector z_temp(num_var_);
    if (x_user)
        std::copy_n(x_user, num_var_, std::begin(x_temp));
    if (slack_user)
        std::copy_n(slack_user, num_constr_, std::begin(slack_temp));
    if (y_user)
        std::copy_n(y_user, num_constr_, std::begin(y_temp));
    if (z_user)
        std::copy_n(z_user, num_var_, std::begin(z_temp));
    ScalePoint(x_temp, slack_temp, y_temp, z_temp);
    DualizeBasicSolution(x_temp, slack_temp, y_temp, z_temp,
                         x_solver, y_solver, z_solver);
}

Int Model::PresolveIPMStartingPoint(const double* x_user,
                                    const double* xl_user,
                                    const double* xu_user,
                                    const double* slack_user,
                                    const double* y_user,
                                    const double* zl_user,
                                    const double* zu_user,
                                    Vector& x_solver,
                                    Vector& xl_solver,
                                    Vector& xu_solver,
                                    Vector& y_solver,
                                    Vector& zl_solver,
                                    Vector& zu_solver) const {
    if (!x_user || !xl_user || !xu_user || !slack_user ||
        !y_user || !zl_user || !zu_user)
        return IPX_ERROR_argument_null;
    if (dualized_)
        return IPX_ERROR_not_implemented;

    // Copy user point into workspace and apply model scaling.
    Vector x_temp(x_user, num_var_);
    Vector xl_temp(xl_user, num_var_);
    Vector xu_temp(xu_user, num_var_);
    Vector slack_temp(slack_user, num_constr_);
    Vector y_temp(y_user, num_constr_);
    Vector zl_temp(zl_user, num_var_);
    Vector zu_temp(zu_user, num_var_);
    ScalePoint(x_temp, xl_temp, xu_temp, slack_temp,  y_temp, zl_temp, zu_temp);

    // Check that point is compatible with bounds.
    for (Int j = 0; j < num_var_; ++j) {
        if (!std::isfinite(x_temp[j]))
            return IPX_ERROR_invalid_vector;
    }
    for (Int j = 0; j < num_var_; ++j) {
        if (!(xl_temp[j] >= 0.0) ||
            scaled_lbuser_[j] == -INFINITY && xl_temp[j] != INFINITY ||
            scaled_lbuser_[j] != -INFINITY && xl_temp[j] == INFINITY)
            return IPX_ERROR_invalid_vector;
    }
    for (Int j = 0; j < num_var_; ++j) {
        if (!(xu_temp[j] >= 0.0) ||
            scaled_ubuser_[j] == INFINITY && xu_temp[j] != INFINITY ||
            scaled_ubuser_[j] != INFINITY && xu_temp[j] == INFINITY)
            return IPX_ERROR_invalid_vector;
    }
    for (Int i = 0; i < num_constr_; ++i) {
        if (!std::isfinite(slack_temp[i]) ||
            constr_type_[i] == '=' && !(slack_temp[i] == 0.0) ||
            constr_type_[i] == '<' && !(slack_temp[i] >= 0.0) ||
            constr_type_[i] == '>' && !(slack_temp[i] <= 0.0) )
            return IPX_ERROR_invalid_vector;
    }
    for (Int i = 0; i < num_constr_; ++i) {
        if (!std::isfinite(y_temp[i]) ||
            constr_type_[i] == '<' && !(y_temp[i] <= 0.0) ||
            constr_type_[i] == '>' && !(y_temp[i] >= 0.0) )
            return IPX_ERROR_invalid_vector;
    }
    for (Int j = 0; j < num_var_; ++j) {
        if (!(zl_temp[j] >= 0.0 && zl_temp[j] < INFINITY) ||
            scaled_lbuser_[j] == -INFINITY && zl_temp[j] != 0.0)
            return IPX_ERROR_invalid_vector;
    }
    for (Int j = 0; j < num_var_; ++j) {
        if (!(zu_temp[j] >= 0.0 && zu_temp[j] < INFINITY) ||
            scaled_ubuser_[j] == INFINITY && zu_temp[j] != 0.0)
            return IPX_ERROR_invalid_vector;
    }

    DualizeIPMStartingPoint(x_temp, xl_temp, xu_temp, slack_temp,
                            y_temp, zl_temp, zu_temp,
                            x_solver, xl_solver, xu_solver,
                            y_solver, zl_solver, zu_solver);
    return 0;
}


void Model::PostsolveInteriorSolution(const Vector& x_solver,
                                      const Vector& xl_solver,
                                      const Vector& xu_solver,
                                      const Vector& y_solver,
                                      const Vector& zl_solver,
                                      const Vector& zu_solver,
                                      double* x_user,
                                      double* xl_user, double* xu_user,
                                      double* slack_user,
                                      double* y_user,
                                      double* zl_user, double* zu_user) const {
    const Int m = rows();
    const Int n = cols();
    assert(x_solver.size() == n+m);
    assert(xl_solver.size() == n+m);
    assert(xu_solver.size() == n+m);
    assert(y_solver.size() == m);
    assert(zl_solver.size() == n+m);
    assert(zu_solver.size() == n+m);

    Vector x_temp(num_var_);
    Vector xl_temp(num_var_);
    Vector xu_temp(num_var_);
    Vector slack_temp(num_constr_);
    Vector y_temp(num_constr_);
    Vector zl_temp(num_var_);
    Vector zu_temp(num_var_);
    DualizeBackInteriorSolution(x_solver, xl_solver, xu_solver, y_solver,
                                zl_solver, zu_solver, x_temp, xl_temp, xu_temp,
                                slack_temp, y_temp, zl_temp, zu_temp);
    ScaleBackInteriorSolution(x_temp, xl_temp, xu_temp, slack_temp, y_temp,
                              zl_temp, zu_temp);
    if (x_user)
        std::copy(std::begin(x_temp), std::end(x_temp), x_user);
    if (xl_user)
        std::copy(std::begin(xl_temp), std::end(xl_temp), xl_user);
    if (xu_user)
        std::copy(std::begin(xu_temp), std::end(xu_temp), xu_user);
    if (slack_user)
        std::copy(std::begin(slack_temp), std::end(slack_temp), slack_user);
    if (y_user)
        std::copy(std::begin(y_temp), std::end(y_temp), y_user);
    if (zl_user)
        std::copy(std::begin(zl_temp), std::end(zl_temp), zl_user);
    if (zu_user)
        std::copy(std::begin(zu_temp), std::end(zu_temp), zu_user);
}

void Model::EvaluateInteriorSolution(const Vector& x_solver,
                                     const Vector& xl_solver,
                                     const Vector& xu_solver,
                                     const Vector& y_solver,
                                     const Vector& zl_solver,
                                     const Vector& zu_solver,
                                     Info* info) const {
    const Int m = rows();
    const Int n = cols();
    assert(x_solver.size() == n+m);
    assert(xl_solver.size() == n+m);
    assert(xu_solver.size() == n+m);
    assert(y_solver.size() == m);
    assert(zl_solver.size() == n+m);
    assert(zu_solver.size() == n+m);

    // Build solution to scaled user model.
    Vector x(num_var_);
    Vector xl(num_var_);
    Vector xu(num_var_);
    Vector slack(num_constr_);
    Vector y(num_constr_);
    Vector zl(num_var_);
    Vector zu(num_var_);
    DualizeBackInteriorSolution(x_solver, xl_solver, xu_solver, y_solver,
                                zl_solver, zu_solver, x, xl, xu, slack, y, zl,
                                zu);

    // Build residuals for scaled model.
    // rl = lb-x+xl
    Vector rl(num_var_);
    for (Int j = 0; j < num_var_; j++) {
        if (std::isfinite(scaled_lbuser_[j]))
            rl[j] = scaled_lbuser_[j] - x[j] + xl[j];
        else
            assert(xl[j] == INFINITY);
    }
    // ru = ub-x-xu
    Vector ru(num_var_);
    for (Int j = 0; j < num_var_; j++) {
        if (std::isfinite(scaled_ubuser_[j]))
            ru[j] = scaled_ubuser_[j] - x[j] - xu[j];
        else
            assert(xu[j] == INFINITY);
    }
    // rb = rhs-slack-A*x
    Vector rb = scaled_rhs_ - slack;
    assert(scaled_rhs_.size() == num_constr_);
    assert(rb.size() == num_constr_);
    MultiplyWithScaledMatrix(x, -1.0, rb, 'N');
    // rc = obj-zl+zu-A'y
    Vector rc = scaled_obj_ - zl + zu;
    MultiplyWithScaledMatrix(y, -1.0, rc, 'T');

    ScaleBackResiduals(rb, rc, rl, ru);
    double presidual = Infnorm(rb);
    presidual = std::max(presidual, Infnorm(rl));
    presidual = std::max(presidual, Infnorm(ru));
    double dresidual = Infnorm(rc);

    double pobjective = Dot(scaled_obj_, x);
    double dobjective = Dot(scaled_rhs_, y);
    for (Int j = 0; j < num_var_; j++) {
        if (std::isfinite(scaled_lbuser_[j]))
            dobjective += scaled_lbuser_[j] * zl[j];
        if (std::isfinite(scaled_ubuser_[j]))
            dobjective -= scaled_ubuser_[j] * zu[j];
    }
    assert(std::isfinite(dobjective));
    double objective_gap = (pobjective-dobjective) /
        (1.0 + 0.5*std::abs(pobjective+dobjective));

    double complementarity = 0.0;
    for (Int j = 0; j < num_var_; j++) {
        if (std::isfinite(scaled_lbuser_[j]))
            complementarity += xl[j]*zl[j];
        if (std::isfinite(scaled_ubuser_[j]))
            complementarity += xu[j]*zu[j];
    }
    for (Int i = 0; i < num_constr_; i++)
        complementarity -= y[i]*slack[i];

    // For computing the norms of the user variables, we have to scale back.
    ScaleBackInteriorSolution(x, xl, xu, slack, y, zl, zu);

    info->abs_presidual = presidual;
    info->abs_dresidual = dresidual;
    info->rel_presidual = presidual/(1.0+norm_rhs_);
    info->rel_dresidual = dresidual/(1.0+norm_obj_);
    info->pobjval = pobjective;
    info->dobjval = dobjective;
    info->rel_objgap = objective_gap;
    info->complementarity = complementarity;
    info->normx = Infnorm(x);
    info->normy = Infnorm(y);
    info->normz = std::max(Infnorm(zl), Infnorm(zu));
}

void Model::PostsolveBasicSolution(const Vector& x_solver,
                                   const Vector& y_solver,
                                   const Vector& z_solver,
                                   const std::vector<Int>& basic_status_solver,
                                   double* x_user, double* slack_user,
                                   double* y_user, double* z_user) const {
    const Int m = rows();
    const Int n = cols();
    assert(x_solver.size() == n+m);
    assert(y_solver.size() == m);
    assert(z_solver.size() == n+m);
    assert(basic_status_solver.size() == n+m);

    Vector x_temp(num_var_);
    Vector slack_temp(num_constr_);
    Vector y_temp(num_constr_);
    Vector z_temp(num_var_);
    std::vector<Int> cbasis_temp(num_constr_);
    std::vector<Int> vbasis_temp(num_var_);
    DualizeBackBasicSolution(x_solver, y_solver, z_solver, x_temp, slack_temp,
                             y_temp, z_temp);
    DualizeBackBasis(basic_status_solver, cbasis_temp, vbasis_temp);
    CorrectScaledBasicSolution(x_temp, slack_temp, y_temp, z_temp, cbasis_temp,
                               vbasis_temp);
    ScaleBackBasicSolution(x_temp, slack_temp, y_temp, z_temp);
    if (x_user)
        std::copy(std::begin(x_temp), std::end(x_temp), x_user);
    if (slack_user)
        std::copy(std::begin(slack_temp), std::end(slack_temp), slack_user);
    if (y_user)
        std::copy(std::begin(y_temp), std::end(y_temp), y_user);
    if (z_user)
        std::copy(std::begin(z_temp), std::end(z_temp), z_user);
}

void Model::EvaluateBasicSolution(const Vector& x_solver,
                                  const Vector& y_solver,
                                  const Vector& z_solver,
                                  const std::vector<Int>& basic_status_solver,
                                  Info* info) const {
    const Int m = rows();
    const Int n = cols();
    assert(x_solver.size() == n+m);
    assert(y_solver.size() == m);
    assert(z_solver.size() == n+m);
    assert(basic_status_solver.size() == n+m);

    // Build basic solution to scaled user model.
    Vector x(num_var_);
    Vector slack(num_constr_);
    Vector y(num_constr_);
    Vector z(num_var_);
    std::vector<Int> cbasis(num_constr_);
    std::vector<Int> vbasis(num_var_);
    DualizeBackBasicSolution(x_solver, y_solver, z_solver, x, slack, y, z);
    DualizeBackBasis(basic_status_solver, cbasis, vbasis);
    CorrectScaledBasicSolution(x, slack, y, z, cbasis, vbasis);
    double pobj = Dot(scaled_obj_, x);

    // Build infeasibilities in scaled user model.
    Vector xinfeas(num_var_);
    Vector sinfeas(num_constr_);
    Vector yinfeas(num_constr_);
    Vector zinfeas(num_var_);
    for (Int j = 0; j < num_var_; j++) {
        if (x[j] < scaled_lbuser_[j])
            xinfeas[j] = x[j]-scaled_lbuser_[j];
        if (x[j] > scaled_ubuser_[j])
            xinfeas[j] = x[j]-scaled_ubuser_[j];
        if (vbasis[j] != IPX_nonbasic_lb && z[j] > 0.0)
            zinfeas[j] = z[j];
        if (vbasis[j] != IPX_nonbasic_ub && z[j] < 0.0)
            zinfeas[j] = z[j];
    }
    for (Int i = 0; i < num_constr_; i++) {
        if (constr_type_[i] == '<') {
            if (slack[i] < 0.0)
                sinfeas[i] = slack[i];
            if (y[i] > 0.0)
                yinfeas[i] = y[i];
        }
        if (constr_type_[i] == '>') {
            if (slack[i] > 0.0)
                sinfeas[i] = slack[i];
            if (y[i] < 0.0)
                yinfeas[i] = y[i];
        }
    }

    // Scale back basic solution and infeasibilities.
    ScaleBackBasicSolution(x, slack, y, z);
    ScaleBackBasicSolution(xinfeas, sinfeas, yinfeas, zinfeas);

    info->primal_infeas = std::max(Infnorm(xinfeas), Infnorm(sinfeas));
    info->dual_infeas = std::max(Infnorm(zinfeas), Infnorm(yinfeas));
    info->objval = pobj;
}

void Model::PostsolveBasis(const std::vector<Int>& basic_status_solver,
                           Int* cbasis_user, Int* vbasis_user) const {
    const Int m = rows();
    const Int n = cols();
    assert(basic_status_solver.size() == n+m);

    std::vector<Int> cbasis_temp(num_constr_);
    std::vector<Int> vbasis_temp(num_var_);
    DualizeBackBasis(basic_status_solver, cbasis_temp, vbasis_temp);
    ScaleBackBasis(cbasis_temp, vbasis_temp);
    if (cbasis_user)
        std::copy(std::begin(cbasis_temp), std::end(cbasis_temp), cbasis_user);
    if (vbasis_user)
        std::copy(std::begin(vbasis_temp), std::end(vbasis_temp), vbasis_user);
}

// Checks if the vectors are valid LP data vectors. Returns 0 if OK and a
// negative value if a vector is invalid.
static int CheckVectors(Int m, Int n, const double* rhs,const char* constr_type,
                        const double* obj, const double* lb, const double* ub) {
    for (Int i = 0; i < m; i++)
        if (!std::isfinite(rhs[i]))
            return -1;
    for (Int j = 0; j < n; j++)
        if (!std::isfinite(obj[j]))
            return -2;
    for (Int j = 0; j < n; j++) {
        if (!std::isfinite(lb[j]) && lb[j] != -INFINITY)
            return -3;
        if (!std::isfinite(ub[j]) && ub[j] != INFINITY)
            return -3;
        if (lb[j] > ub[j])
            return -3;
    }
    for (Int i = 0; i < m; i++)
        if (constr_type[i] != '=' && constr_type[i] != '<' &&
            constr_type[i] != '>')
            return -4;
    return 0;
}

// Checks if A is a valid m-by-n matrix in CSC format. Returns 0 if OK and a
// negative value otherwise.
Int CheckMatrix(Int m, Int n, const Int *Ap, const Int *Ai, const double *Ax) {
    if (Ap[0] != 0)
        return -5;
    for (Int j = 0; j < n; j++)
        if (Ap[j] > Ap[j+1])
            return -5;
    for (Int p = 0; p < Ap[n]; p++)
        if (!std::isfinite(Ax[p]))
            return -6;
    // Test for out of bound indices and duplicates.
    std::vector<Int> marked(m, -1);
    for (Int j = 0; j < n; j++) {
        for (Int p = Ap[j]; p < Ap[j+1]; p++) {
            Int i = Ai[p];
            if (i < 0 || i >= m)
                return -7;
            if (marked[i] == j)
                return -8;
            marked[i] = j;
        }
    }
    return 0;
}

Int Model::CopyInput(Int num_constr, Int num_var, const Int* Ap, const Int* Ai,
                     const double* Ax, const double* rhs,
                     const char* constr_type, const double* obj,
                     const double* lbuser, const double* ubuser) {
    if (!(Ap && Ai && Ax && rhs && constr_type && obj && lbuser && ubuser)) {
        return IPX_ERROR_argument_null;
    }
    if (num_constr < 0 || num_var <= 0) {
        return IPX_ERROR_invalid_dimension;
    }
    if (CheckVectors(num_constr, num_var, rhs, constr_type, obj, lbuser, ubuser)
        != 0) {
        return IPX_ERROR_invalid_vector;
    }
    if (CheckMatrix(num_constr, num_var, Ap, Ai, Ax) != 0) {
        return IPX_ERROR_invalid_matrix;
    }
    num_constr_ = num_constr;
    num_eqconstr_ = std::count(constr_type, constr_type+num_constr, '=');
    num_var_ = num_var;
    num_entries_ = Ap[num_var];
    num_free_var_ = 0;
    boxed_vars_.clear();
    for (Int j = 0; j < num_var; j++) {
        if (std::isinf(lbuser[j]) && std::isinf(ubuser[j]))
            num_free_var_++;
        if (std::isfinite(lbuser[j]) && std::isfinite(ubuser[j]))
            boxed_vars_.push_back(j);
    }
    constr_type_ = std::vector<char>(constr_type, constr_type+num_constr);
    scaled_obj_ = Vector(obj, num_var);
    scaled_rhs_ = Vector(rhs, num_constr);
    scaled_lbuser_ = Vector(lbuser, num_var);
    scaled_ubuser_ = Vector(ubuser, num_var);
    A_.LoadFromArrays(num_constr, num_var, Ap, Ap+1, Ai, Ax);
    norm_obj_ = Infnorm(scaled_obj_);
    norm_rhs_ = Infnorm(scaled_rhs_);
    for (double x : scaled_lbuser_)
        if (std::isfinite(x))
            norm_rhs_ = std::max(norm_rhs_, std::abs(x));
    for (double x : scaled_ubuser_)
        if (std::isfinite(x))
            norm_rhs_ = std::max(norm_rhs_, std::abs(x));
    return 0;
}

void Model::ScaleModel(const Control& control) {
    flipped_vars_.clear();
    for (Int j = 0; j < num_var_; j++) {
        if (std::isfinite(scaled_ubuser_[j]) && std::isinf(scaled_lbuser_[j])) {
            scaled_lbuser_[j] = -scaled_ubuser_[j];
            scaled_ubuser_[j] = INFINITY;
            ScaleColumn(A_, j, -1.0);
            scaled_obj_[j] *= -1.0;
            flipped_vars_.push_back(j);
        }
    }
    colscale_.resize(0);
    rowscale_.resize(0);

    // Choose scaling method.
    if (control.scale() >= 1)
        EquilibrateMatrix();

    // Apply scaling to vectors.
    if (colscale_.size() > 0) {
        assert(colscale_.size() == num_var_);
        scaled_obj_ *= colscale_;
        scaled_lbuser_ /= colscale_;
        scaled_ubuser_ /= colscale_;
    }
    if (rowscale_.size() > 0) {
        assert(rowscale_.size() == num_constr_);
        scaled_rhs_ *= rowscale_;
    }
}

void Model::LoadPrimal() {
    num_rows_ = num_constr_;
    num_cols_ = num_var_;
    dualized_ = false;

    // Copy A and append identity matrix.
    AI_ = A_;
    for (Int i = 0; i < num_constr_; i++) {
        AI_.push_back(i, 1.0);
        AI_.add_column();
    }
    assert(AI_.cols() == num_var_+num_constr_);

    // Copy vectors and set bounds on slack variables.
    b_ = scaled_rhs_;
    c_.resize(num_var_+num_constr_);
    c_ = 0.0;
    std::copy_n(std::begin(scaled_obj_), num_var_, std::begin(c_));
    lb_.resize(num_rows_+num_cols_);
    std::copy_n(std::begin(scaled_lbuser_), num_var_, std::begin(lb_));
    ub_.resize(num_rows_+num_cols_);
    std::copy_n(std::begin(scaled_ubuser_), num_var_, std::begin(ub_));
    for (Int i = 0; i < num_constr_; i++) {
        switch(constr_type_[i]) {
        case '=':
            lb_[num_var_+i] = 0.0;
            ub_[num_var_+i] = 0.0;
            break;
        case '<':
            lb_[num_var_+i] = 0.0;
            ub_[num_var_+i] = INFINITY;
            break;
        case '>':
            lb_[num_var_+i] = -INFINITY;
            ub_[num_var_+i] = 0.0;
            break;
        }
    }
}

void Model::LoadDual() {
    num_rows_ = num_var_;
    num_cols_ = num_constr_ + boxed_vars_.size();
    dualized_ = true;

    // Check that every variable with finite scaled_ubuser_ has finite
    // scaled_lbuser_ (must be the case after scaling).
    for (Int j = 0; j < num_var_; j++) {
        if (std::isfinite(scaled_ubuser_[j]))
            assert(std::isfinite(scaled_lbuser_[j]));
    }

    // Build AI.
    AI_ = Transpose(A_);
    for (Int j = 0; j < num_var_; j++) {
        if (std::isfinite(scaled_ubuser_[j])) {
            AI_.push_back(j, -1.0);
            AI_.add_column();
        }
    }
    assert(AI_.cols() == num_cols_);
    for (Int i = 0; i < num_rows_; i++) {
        AI_.push_back(i, 1.0);
        AI_.add_column();
    }

    // Build vectors.
    b_ = scaled_obj_;
    c_.resize(num_cols_+num_rows_);
    Int put = 0;
    for (double x : scaled_rhs_)
        c_[put++] = -x;
    for (double x : scaled_ubuser_)
        if (std::isfinite(x))
            c_[put++] = x;
    assert(put == num_cols_);
    for (double x : scaled_lbuser_)
        // If x is negative infinity, then the variable will be fixed and we can
        // give it any (finite) cost.
        c_[put++] = std::isfinite(x) ? -x : 0.0;
    lb_.resize(num_cols_+num_rows_);
    ub_.resize(num_cols_+num_rows_);
    for (Int i = 0; i < num_constr_; i++)
        switch(constr_type_[i]) {
        case '=':
            lb_[i] = -INFINITY;
            ub_[i] = INFINITY;
            break;
        case '<':
            lb_[i] = -INFINITY;
            ub_[i] = 0.0;
            break;
        case '>':
            lb_[i] = 0.0;
            ub_[i] = INFINITY;
            break;
        }
    for (Int j = num_constr_; j < num_cols_; j++) {
        lb_[j] = 0.0;
        ub_[j] = INFINITY;
    }
    for (Int j = 0; j < num_var_; j++) {
        lb_[num_cols_+j] = 0.0;
        ub_[num_cols_+j] = std::isfinite(scaled_lbuser_[j]) ? INFINITY : 0.0;
    }
}

// Returns a power-of-2 factor s such that s*2^exp becomes closer to the
// interval [2^expmin, 2^expmax].
static double EquilibrationFactor(int expmin, int expmax, int exp) {
    if (exp < expmin)
        return std::ldexp(1.0, (expmin-exp+1)/2);
    if (exp > expmax)
        return std::ldexp(1.0, -((exp-expmax+1)/2));
    return 1.0;
}

void Model::EquilibrateMatrix() {
    const Int m = A_.rows();
    const Int n = A_.cols();
    const Int* Ap = A_.colptr();
    const Int* Ai = A_.rowidx();
    double* Ax = A_.values();

    colscale_.resize(0);
    rowscale_.resize(0);

    // The absolute value of each nonzero entry of AI can be written as x*2^exp
    // with x in the range [0.5,1). We consider AI well scaled if each entry is
    // such that expmin <= exp <= expmax for parameters expmin and expmax.
    // For example,
    //
    //  expmin  expmax  allowed range
    //  -----------------------------
    //       0       2    [0.5,    4)
    //      -1       3    [0.25,   8)
    //       1       7    [1.0,  128)
    //
    // If AI is not well scaled, a recursive row and column equilibration is
    // applied. In each iteration, the scaling factors are roughly 1/sqrt(max),
    // where max is the maximum row or column entry. However, the factors are
    // truncated to powers of 2, so that no round-off errors occur.

    constexpr int expmin = 0;
    constexpr int expmax = 3;
    constexpr Int maxround = 10;

    // Quick return if entries are within the target range.
    bool out_of_range = false;
    for (Int p = 0; p < Ap[n]; p++) {
        int exp;
        std::frexp(std::abs(Ax[p]), &exp);
        if (exp < expmin || exp > expmax) {
            out_of_range = true;
            break;
        }
    }
    if (!out_of_range)
        return;

    colscale_.resize(n);
    rowscale_.resize(m);
    colscale_ = 1.0;
    rowscale_ = 1.0;
    Vector colmax(n), rowmax(m);

    for (Int round = 0; round < maxround; round++) {
        // Compute infinity norm of each row and column.
        rowmax = 0.0;
        for (Int j = 0; j < n; j++) {
            colmax[j] = 0.0;
            for (Int p = Ap[j]; p < Ap[j+1]; p++) {
                Int i = Ai[p];
                double xa = std::abs(Ax[p]);
                colmax[j] = std::max(colmax[j], xa);
                rowmax[i] = std::max(rowmax[i], xa);
            }
        }
        // Replace rowmax and colmax entries by scaling factors from this round.
        bool out_of_range = false;
        for (Int i = 0; i < m; i++) {
            int exp;
            std::frexp(rowmax[i], &exp);
            rowmax[i] = EquilibrationFactor(expmin, expmax, exp);
            if (rowmax[i] != 1.0) {
                out_of_range = true;
                rowscale_[i] *= rowmax[i];
            }
        }
        for (Int j = 0; j < n; j++) {
            int exp;
            std::frexp(colmax[j], &exp);
            colmax[j] = EquilibrationFactor(expmin, expmax, exp);
            if (colmax[j] != 1.0) {
                out_of_range = true;
                colscale_[j] *= colmax[j];
            }
        }
        if (!out_of_range)
            break;
        // Rescale A.
        for (Int j = 0; j < n; j++) {
            for (Int p = Ap[j]; p < Ap[j+1]; p++) {
                Ax[p] *= colmax[j];     // column scaling
                Ax[p] *= rowmax[Ai[p]]; // row scaling
            }
        }
    }
}

void Model::FindDenseColumns() {
    num_dense_cols_ = 0;
    nz_dense_ = rows() + 1;

    std::vector<Int> colcount(num_cols_);
    for (Int j = 0; j < num_cols_; j++)
        colcount[j] = AI_.end(j)-AI_.begin(j);
    std::sort(colcount.begin(), colcount.end());

    for (Int j = 1; j < num_cols_; j++) {
        if (colcount[j] >
              std::max((ipx::Int)40l, (ipx::Int)10l*colcount[j-1])) {
            // j is the first dense column
            num_dense_cols_ = num_cols_ - j;
            nz_dense_ = colcount[j];
            break;
        }
    }

    if (num_dense_cols_ > 1000) {
        num_dense_cols_ = 0;
        nz_dense_ = rows() + 1;
    }
}

void Model::PrintCoefficientRange(const Control& control) const {
    double amin = INFINITY;
    double amax = 0.0;
    for (Int j = 0; j < A_.cols(); j++) {
        for (Int p = A_.begin(j); p < A_.end(j); p++) {
            double x = A_.value(p);
            if (x != 0.0) {
                amin = std::min(amin, std::abs(x));
                amax = std::max(amax, std::abs(x));
            }
        }
    }
    if (amin == INFINITY)       // no nonzero entries in A_
        amin = 0.0;
    control.Log()
        << Textline("Matrix range:")
        << "[" << Scientific(amin, 5, 0) << ", "
        << Scientific(amax, 5, 0) << "]\n";

    double rhsmin = INFINITY;
    double rhsmax = 0.0;
    for (double x : scaled_rhs_) {
        if (x != 0.0) {
            rhsmin = std::min(rhsmin, std::abs(x));
            rhsmax = std::max(rhsmax, std::abs(x));
        }
    }
    if (rhsmin == INFINITY)     // no nonzero entries in rhs
        rhsmin = 0.0;
    control.Log()
        << Textline("RHS range:")
        << "[" << Scientific(rhsmin, 5, 0) << ", "
        << Scientific(rhsmax, 5, 0) << "]\n";

    double objmin = INFINITY;
    double objmax = 0.0;
    for (double x : scaled_obj_) {
        if (x != 0.0) {
            objmin = std::min(objmin, std::abs(x));
            objmax = std::max(objmax, std::abs(x));
        }
    }
    if (objmin == INFINITY)     // no nonzero entries in obj
        objmin = 0.0;
    control.Log()
        << Textline("Objective range:")
        << "[" << Scientific(objmin, 5, 0) << ", "
        << Scientific(objmax, 5, 0) << "]\n";

    double boundmin = INFINITY;
    double boundmax = 0.0;
    for (double x : scaled_lbuser_) {
        if (x != 0.0 && std::isfinite(x)) {
            boundmin = std::min(boundmin, std::abs(x));
            boundmax = std::max(boundmax, std::abs(x));
        }
    }
    for (double x : scaled_ubuser_) {
        if (x != 0.0 && std::isfinite(x)) {
            boundmin = std::min(boundmin, std::abs(x));
            boundmax = std::max(boundmax, std::abs(x));
        }
    }
    if (boundmin == INFINITY)   // no finite nonzeros entries in bounds
        boundmin = 0.0;
    control.Log()
        << Textline("Bounds range:")
        << "[" << Scientific(boundmin, 5, 0) << ", "
        << Scientific(boundmax, 5, 0) << "]\n";
}

void Model::PrintPreprocessingLog(const Control& control) const {
    // Find the minimum and maximum scaling factor.
    double minscale = INFINITY;
    double maxscale = 0.0;
    if (colscale_.size() > 0) {
        auto minmax = std::minmax_element(std::begin(colscale_),
                                          std::end(colscale_));
        minscale = std::min(minscale, *minmax.first);
        maxscale = std::max(maxscale, *minmax.second);
    }
    if (rowscale_.size() > 0) {
        auto minmax = std::minmax_element(std::begin(rowscale_),
                                          std::end(rowscale_));
        minscale = std::min(minscale, *minmax.first);
        maxscale = std::max(maxscale, *minmax.second);
    }
    if (minscale == INFINITY)
        minscale = 1.0;
    if (maxscale == 0.0)
        maxscale = 1.0;

    control.Log()
        << "Preprocessing\n"
        << Textline("Dualized model:") << (dualized() ? "yes" : "no") << '\n'
        << Textline("Number of dense columns:") << num_dense_cols() << '\n';
    if (control.scale() > 0) {
        control.Log()
            << Textline("Range of scaling factors:") << "["
            << Scientific(minscale, 8, 2) << ", "
            << Scientific(maxscale, 8, 2) << "]\n";
    }
}

void Model::ScalePoint(Vector& x, Vector& slack, Vector& y, Vector& z) const {
    if (colscale_.size() > 0) {
        x /= colscale_;
        z *= colscale_;
    }
    if (rowscale_.size() > 0) {
        y /= rowscale_;
        slack *= rowscale_;
    }
    for (Int j : flipped_vars_) {
        x[j] *= -1.0;
        z[j] *= -1.0;
    }
}

void Model::ScalePoint(Vector& x, Vector& xl, Vector& xu, Vector& slack,
                       Vector& y, Vector& zl, Vector& zu) const {
    if (colscale_.size() > 0) {
        x /= colscale_;
        xl /= colscale_;
        xu /= colscale_;
        zl *= colscale_;
        zu *= colscale_;
    }
    if (rowscale_.size() > 0) {
        y /= rowscale_;
        slack *= rowscale_;
    }
    for (Int j : flipped_vars_) {
        x[j] *= -1.0;
        xl[j] = xu[j];
        xu[j] = INFINITY;
        zl[j] = zu[j];
        zu[j] = 0.0;
    }
}

void Model::ScaleBackInteriorSolution(Vector& x, Vector& xl, Vector& xu,
                                      Vector& slack, Vector& y, Vector& zl,
                                      Vector& zu) const {
    if (colscale_.size() > 0) {
        x *= colscale_;
        xl *= colscale_;
        xu *= colscale_;
        zl /= colscale_;
        zu /= colscale_;
    }
    if (rowscale_.size() > 0) {
        y *= rowscale_;
        slack /= rowscale_;
    }
    for (Int j : flipped_vars_) {
        assert(std::isfinite(xl[j]));
        assert(std::isinf(xu[j]));
        assert(zu[j] == 0.0);
        x[j] *= -1.0;
        xu[j] = xl[j];
        xl[j] = INFINITY;
        zu[j] = zl[j];
        zl[j] = 0.0;
    }
}

void Model::ScaleBackResiduals(Vector& rb, Vector& rc, Vector& rl,
                               Vector& ru) const {
    if (colscale_.size() > 0) {
        rc /= colscale_;
        rl *= colscale_;
        ru *= colscale_;
    }
    if (rowscale_.size() > 0)
        rb /= rowscale_;
    for (Int j : flipped_vars_) {
        rc[j] *= -1.0;
        assert(ru[j] == 0.0);
        ru[j] = -rl[j];
        rl[j] = 0.0;
    }
}

void Model::ScaleBackBasicSolution(Vector& x, Vector& slack, Vector& y,
                                   Vector& z) const {
    if (colscale_.size() > 0) {
        x *= colscale_;
        z /= colscale_;
    }
    if (rowscale_.size() > 0) {
        y *= rowscale_;
        slack /= rowscale_;
    }
    for (Int j : flipped_vars_) {
        x[j] *= -1.0;
        z[j] *= -1.0;
    }
}

void Model::ScaleBackBasis(std::vector<Int>& cbasis,
                           std::vector<Int>& vbasis) const {
    for (Int j : flipped_vars_) {
        assert(vbasis[j] != IPX_nonbasic_ub);
        if (vbasis[j] == IPX_nonbasic_lb)
            vbasis[j] = IPX_nonbasic_ub;
    }
}

void Model::DualizeBasicSolution(const Vector& x_user,
                                 const Vector& slack_user,
                                 const Vector& y_user,
                                 const Vector& z_user,
                                 Vector& x_solver,
                                 Vector& y_solver,
                                 Vector& z_solver) const {
    const Int m = rows();
    const Int n = cols();

    if (dualized_) {
        assert(num_var_ == m);
        assert(num_constr_ + boxed_vars_.size() == n);

        // Build dual solver variables from primal user variables.
        y_solver = -x_user;
        for (Int i = 0; i < num_constr_; i++)
            z_solver[i] = -slack_user[i];
        for (Int k = 0; k < boxed_vars_.size(); k++) {
            Int j = boxed_vars_[k];
            z_solver[num_constr_+k] = c(num_constr_+k) + y_solver[j];
        }
        for (Int i = 0; i < m; i++)
            z_solver[n+i] = c(n+i)-y_solver[i];

        // Build primal solver variables from dual user variables.
        std::copy_n(std::begin(y_user), num_constr_, std::begin(x_solver));
        std::copy_n(std::begin(z_user), num_var_, std::begin(x_solver) + n);
        for (Int k = 0; k < boxed_vars_.size(); k++) {
            Int j = boxed_vars_[k];
            if (x_solver[n+j] < 0.0) {
                // j is a boxed variable and z_user[j] < 0
                x_solver[num_constr_+k] = -x_solver[n+j];
                x_solver[n+j] = 0.0;
            } else {
                x_solver[num_constr_+k] = 0.0;
            }
        }
    }
    else {
        assert(num_constr_ == m);
        assert(num_var_ == n);
        std::copy_n(std::begin(x_user), n, std::begin(x_solver));
        std::copy_n(std::begin(slack_user), m, std::begin(x_solver) + n);
        std::copy_n(std::begin(y_user), m, std::begin(y_solver));
        std::copy_n(std::begin(z_user), n, std::begin(z_solver));
        for (Int i = 0; i < m; i++)
            z_solver[n+i] = c(n+i)-y_solver[i];
    }
}

void Model::DualizeIPMStartingPoint(const Vector& x_user,
                                    const Vector& xl_user,
                                    const Vector& xu_user,
                                    const Vector& slack_user,
                                    const Vector& y_user,
                                    const Vector& zl_user,
                                    const Vector& zu_user,
                                    Vector& x_solver,
                                    Vector& xl_solver,
                                    Vector& xu_solver,
                                    Vector& y_solver,
                                    Vector& zl_solver,
                                    Vector& zu_solver) const {
    const Int m = rows();
    const Int n = cols();

    if (dualized_) {
        // Not implemented.
        assert(false);
    }
    else {
        assert(num_constr_ == m);
        assert(num_var_ == n);

        std::copy_n(std::begin(x_user), num_var_, std::begin(x_solver));
        std::copy_n(std::begin(slack_user), num_constr_,
                    std::begin(x_solver) + n);
        std::copy_n(std::begin(xl_user), num_var_, std::begin(xl_solver));
        std::copy_n(std::begin(xu_user), num_var_, std::begin(xu_solver));
        std::copy_n(std::begin(y_user), num_constr_, std::begin(y_solver));
        std::copy_n(std::begin(zl_user), num_var_, std::begin(zl_solver));
        std::copy_n(std::begin(zu_user), num_var_, std::begin(zu_solver));

        for (Int i = 0; i < m; i++) {
            switch (constr_type_[i]) {
            case '=':
                assert(lb_[n+i] == 0.0 && ub_[n+i] == 0.0);
                // For a fixed slack variable xl, xu, zl and zu won't be used
                // by the IPM. Just put them to zero.
                xl_solver[n+i] = 0.0;
                xu_solver[n+i] = 0.0;
                zl_solver[n+i] = 0.0;
                zu_solver[n+i] = 0.0;
                break;
            case '<':
                assert(lb_[n+i] == 0.0 && ub_[n+i] == INFINITY);
                xl_solver[n+i] = slack_user[i];
                xu_solver[n+i] = INFINITY;
                zl_solver[n+i] = -y_user[i];
                zu_solver[n+i] = 0.0;
                break;
            case '>':
                assert(lb_[n+i] == -INFINITY && ub_[n+i] == 0.0);
                xl_solver[n+i] = INFINITY;
                xu_solver[n+i] = -slack_user[i];
                zl_solver[n+i] = 0.0;
                zu_solver[n+i] = y_user[i];
                break;
            }
        }
    }
}

void Model::DualizeBackInteriorSolution(const Vector& x_solver,
                                        const Vector& xl_solver,
                                        const Vector& xu_solver,
                                        const Vector& y_solver,
                                        const Vector& zl_solver,
                                        const Vector& zu_solver,
                                        Vector& x_user,
                                        Vector& xl_user,
                                        Vector& xu_user,
                                        Vector& slack_user,
                                        Vector& y_user,
                                        Vector& zl_user,
                                        Vector& zu_user) const {
    const Int m = rows();
    const Int n = cols();

    if (dualized_) {
        assert(num_var_ == m);
        assert(num_constr_ + boxed_vars_.size() == n);
        x_user = -y_solver;

        // If the solution from the solver would be exact, we could copy the
        // first num_constr_ entries from x_solver into y_user. However, to
        // satisfy the sign condition on y_user even if the solution is not
        // exact, we have to use the xl_solver and xu_solver entries for
        // inequality constraints.
        for (Int i = 0; i < num_constr_; i++) {
            switch (constr_type_[i]) {
            case '=':
                y_user[i] = x_solver[i];
                break;
            case '<':
                y_user[i] = -xu_solver[i];
                break;
            case '>':
                y_user[i] = xl_solver[i];
                break;
            }
            assert(std::isfinite(y_user[i]));
        }

        // Dual variables associated with lbuser <= x in the scaled user model
        // are the slack variables from the solver. For an exact solution we
        // would have x_solver[n+1:n+m] == xl_solver[n+1:n+m]. Using xl_solver
        // guarantees that zl_user >= 0 in any case. If variable j has no lower
        // bound in the scaled user model (i.e. is free), then the j-th slack
        // variable was fixed at zero in the solver model, but the IPM solution
        // may not satisfy this. Hence we must set zl_user[j] = 0 explicitly.
        std::copy_n(std::begin(xl_solver) + n, num_var_, std::begin(zl_user));
        for (Int j = 0; j < num_var_; j++)
            if (!std::isfinite(scaled_lbuser_[j]))
                zl_user[j] = 0.0;

        // Dual variables associated with x <= ubuser in the scaled user model
        // are the primal variables that were added for boxed variables in the
        // solver model.
        zu_user = 0.0;
        Int k = num_constr_;
        for (Int j : boxed_vars_)
            zu_user[j] = xl_solver[k++];
        assert(k == n);

        // xl in the scaled user model is zl[n+1:n+m] in the solver model or
        // infinity.
        for (Int i = 0; i < m; i++) {
            if (std::isfinite(scaled_lbuser_[i]))
                xl_user[i] = zl_solver[n+i];
            else
                xl_user[i] = INFINITY;
        }

        // xu in the scaled user model are the entries in zl for columns of the
        // negative identity matrix (that were added for boxed variables).
        xu_user = INFINITY;
        k = num_constr_;
        for (Int j : boxed_vars_)
            xu_user[j] = zl_solver[k++];
        assert(k == n);

        for (Int i = 0; i < num_constr_; i++) {
            switch (constr_type_[i]) {
            case '=':
                slack_user[i] = 0.0;
                break;
            case '<':
                slack_user[i] = zu_solver[i];
                break;
            case '>':
                slack_user[i] = -zl_solver[i];
                break;
            }
        }
    }
    else {
        assert(num_constr_ == m);
        assert(num_var_ == n);
        std::copy_n(std::begin(x_solver), num_var_, std::begin(x_user));

        // Instead of copying y_solver into y_user, we use the entries from
        // zl_solver and zu_solver for inequality constraints, so that the sign
        // condition on y_user is satisfied.
        for (Int i = 0; i < m; i++) {
            assert(lb_[n+i] == 0.0 || lb_[n+i] == -INFINITY);
            assert(ub_[n+i] == 0.0 || ub_[n+i] ==  INFINITY);
            assert(lb_[n+i] == 0.0 || ub_[n+i] == 0.0);
            switch (constr_type_[i]) {
            case '=':
                y_user[i] = y_solver[i];
                break;
            case '<':
                y_user[i] = -zl_solver[n+i];
                break;
            case '>':
                y_user[i] = zu_solver[n+i];
                break;
            }
            assert(std::isfinite(y_user[i]));
        }
        std::copy_n(std::begin(zl_solver), num_var_, std::begin(zl_user));
        std::copy_n(std::begin(zu_solver), num_var_, std::begin(zu_user));
        std::copy_n(std::begin(xl_solver), num_var_, std::begin(xl_user));
        std::copy_n(std::begin(xu_solver), num_var_, std::begin(xu_user));

        // If the solution would be exact, slack_user were given by the entries
        // of x_solver corresponding to slack columns. To satisfy the sign
        // condition in any case, we build the slack for inequality constraints
        // from xl_solver and xu_solver and set the slack for equality
        // constraints to zero.
        for (Int i = 0; i < m; i++) {
            switch (constr_type_[i]) {
            case '=':
                slack_user[i] = 0.0;
                break;
            case '<':
                slack_user[i] = xl_solver[n+i];
                break;
            case '>':
                slack_user[i] = -xu_solver[n+i];
                break;
            }
            assert(std::isfinite(slack_user[i]));
        }
    }
}

void Model::DualizeBackBasicSolution(const Vector& x_solver,
                                     const Vector& y_solver,
                                     const Vector& z_solver,
                                     Vector& x_user,
                                     Vector& slack_user,
                                     Vector& y_user,
                                     Vector& z_user) const {
    const Int m = rows();
    const Int n = cols();

    if (dualized_) {
        assert(num_var_ == m);
        assert(num_constr_ + boxed_vars_.size() == n);
        x_user = -y_solver;
        for (Int i = 0; i < num_constr_; i++)
            slack_user[i] = -z_solver[i];
        std::copy_n(std::begin(x_solver), num_constr_, std::begin(y_user));
        std::copy_n(std::begin(x_solver) + n, num_var_, std::begin(z_user));
        Int k = num_constr_;
        for (Int j : boxed_vars_)
            z_user[j] -= x_solver[k++];
        assert(k == n);
    }
    else {
        assert(num_constr_ == m);
        assert(num_var_ == n);
        std::copy_n(std::begin(x_solver), num_var_, std::begin(x_user));
        std::copy_n(std::begin(x_solver) + n, num_constr_,
                    std::begin(slack_user));
        std::copy_n(std::begin(y_solver), num_constr_, std::begin(y_user));
        std::copy_n(std::begin(z_solver), num_var_, std::begin(z_user));
    }
}

void Model::DualizeBackBasis(const std::vector<Int>& basic_status_solver,
                             std::vector<Int>& cbasis_user,
                             std::vector<Int>& vbasis_user) const {
    const Int m = rows();
    const Int n = cols();

    if (dualized_) {
        assert(num_var_ == m);
        assert(num_constr_ + boxed_vars_.size() == n);
        for (Int i = 0; i < num_constr_; i++) {
            if (basic_status_solver[i] == IPX_basic)
                cbasis_user[i] = IPX_nonbasic;
            else
                cbasis_user[i] = IPX_basic;
        }
        for (Int j = 0; j < num_var_; j++) {
            // slack cannot be superbasic
            assert(basic_status_solver[n+j] != IPX_superbasic);
            if (basic_status_solver[n+j] == 0)
                vbasis_user[j] = std::isfinite(scaled_lbuser_[j]) ?
                    IPX_nonbasic_lb : IPX_superbasic;
            else
                vbasis_user[j] = IPX_basic;
        }
        Int k = num_constr_;
        for (Int j : boxed_vars_)
            if (basic_status_solver[k++] == IPX_basic) {
                assert(vbasis_user[j] == IPX_basic);
                vbasis_user[j] = IPX_nonbasic_ub;
            }
    }
    else {
        assert(num_constr_ == m);
        assert(num_var_ == n);
        for (Int i = 0; i < num_constr_; i++) {
            // slack cannot be superbasic
            assert(basic_status_solver[n+i] != IPX_superbasic);
            if (basic_status_solver[n+i] == IPX_basic)
                cbasis_user[i] = IPX_basic;
            else
                cbasis_user[i] = IPX_nonbasic;
        }
        for (Int j = 0; j < num_var_; j++)
            vbasis_user[j] = basic_status_solver[j];
    }
}

void Model::CorrectScaledBasicSolution(Vector& x, Vector& slack, Vector& y,
                                       Vector& z,
                                       const std::vector<Int> cbasis,
                                       const std::vector<Int> vbasis) const {
    for (Int j = 0; j < num_var_; j++) {
        if (vbasis[j] == IPX_nonbasic_lb)
            x[j] = scaled_lbuser_[j];
        if (vbasis[j] == IPX_nonbasic_ub)
            x[j] = scaled_ubuser_[j];
        if (vbasis[j] == IPX_basic)
            z[j] = 0.0;
    }
    for (Int i = 0; i < num_constr_; i++) {
        if (cbasis[i] == IPX_nonbasic)
            slack[i] = 0.0;
        if (cbasis[i] == IPX_basic)
            y[i] = 0.0;
    }
}

void Model::MultiplyWithScaledMatrix(const Vector& rhs, double alpha,
                                     Vector& lhs, char trans) const {
    if (trans == 't' || trans == 'T') {
        assert(rhs.size() == num_constr_);
        assert(lhs.size() == num_var_);
        if (dualized())
            for (Int i = 0; i < num_constr_; i++)
                ScatterColumn(AI_, i, alpha*rhs[i], lhs);
        else
            for (Int j = 0; j < num_var_; j++)
                lhs[j] += alpha * DotColumn(AI_, j, rhs);
    }
    else {
        assert(rhs.size() == num_var_);
        assert(lhs.size() == num_constr_);
        if (dualized())
            for (Int i = 0; i < num_constr_; i++)
                lhs[i] += alpha * DotColumn(AI_, i, rhs);
        else
            for (Int j = 0; j < num_var_; j++)
                ScatterColumn(AI_, j, alpha*rhs[j], lhs);
    }
}

double PrimalInfeasibility(const Model& model, const Vector& x) {
    const Vector& lb = model.lb();
    const Vector& ub = model.ub();
    assert(x.size() == lb.size());

    double infeas = 0.0;
    for (Int j = 0; j < x.size(); j++) {
        infeas = std::max(infeas, lb[j]-x[j]);
        infeas = std::max(infeas, x[j]-ub[j]);
    }
    return infeas;
}

double DualInfeasibility(const Model& model, const Vector& x,
                                const Vector& z) {
    const Vector& lb = model.lb();
    const Vector& ub = model.ub();
    assert(x.size() == lb.size());
    assert(z.size() == lb.size());

    double infeas = 0.0;
    for (Int j = 0; j < x.size(); j++) {
        if (x[j] > lb[j])
            infeas = std::max(infeas, z[j]);
        if (x[j] < ub[j])
            infeas = std::max(infeas, -z[j]);
    }
    return infeas;
}

double PrimalResidual(const Model& model, const Vector& x) {
    const SparseMatrix& AIt = model.AIt();
    const Vector& b = model.b();
    assert(x.size() == AIt.rows());

    double res = 0.0;
    for (Int i = 0; i < b.size(); i++) {
        double r = b[i] - DotColumn(AIt, i, x);
        res = std::max(res, std::abs(r));
    }
    return res;
}

double DualResidual(const Model& model, const Vector& y, const Vector& z) {
    const SparseMatrix& AI = model.AI();
    const Vector& c = model.c();
    assert(y.size() == AI.rows());
    assert(z.size() == AI.cols());

    double res = 0.0;
    for (Int j = 0; j < c.size(); j++) {
        double r = c[j] - z[j] - DotColumn(AI, j, y);
        res = std::max(res, std::abs(r));
    }
    return res;
}

}  // namespace ipx