Square Linear Solver

[bknueven]

Summary

Add a solver to Pyomo for square linear systems

Rationale

Pyomo can be used to build systems of linear equations expressively, but as far as I am aware you need to call out to an optimization solver (or one of the nonlinear scipy solvers) to actually solve one.

Description

I have a bare-bones implementation utilizing the LinearStandardFormCompiler:

from pyomo.repn.plugins.standard_form import LinearStandardFormCompiler
from scipy.sparse.linalg import spsolve

class ScipySquareLinearSolver:

    @staticmethod
    def solve(model, **kwds):
        repn = LinearStandardFormCompiler().write(model, slack_form=True, set_sense=None)

        # system must be square
        assert len(repn.rows) == len(repn.columns)

        # TODO: check variable bounds, integer constraints, etc.

        x = spsolve(repn.A, repn.rhs)
    
        for idx, var in enumerate(repn.columns):
            var.set_value(x[idx], skip_validation=True)

        # TODO: return a solution status
        return None

Additional information

In many ways this seems silly to me, but the motivating model for this has ~8 million variables and constraints, so performance is somewhat of a concern -- calling out to an optimization solver (esp. a freely available one) bears a lot of unneeded overhead and complexity. I'm happy to let this live in our codebase and not Pyomo, but maybe somebody else would find it useful?