wiki:SolverDesign

solutions/solve.h must be fixed and enriched. Pascal made a proposal in this direction that I can't put my finger on right now. Among other things he wanted to provide left and right versions. An alternative to that is to ask the user to transpose her matrix when wanting a left solution. However, here I'm going with left and right formss. Here is my proposal.

  1. = Linbox 2.0 Linear System Solver Functions =

For each of these functions another exists with "Right" replaced by "Left". Also, each of these has the suite of versions as described in SolutionDesign.

  • solveRight(x, A, b): provides a solution to Ax=b. It throws an error if the system is inconsistent. When the system has more than one solution, an arbitrary one (implementor's choice) is provided.
  • solveRight(x, A, b, m): same but inconsistency is flagged in the method m, no error is thrown.
  • nullspaceVectorRight(x, A): A random right nullspace vector is returned in x. Random means that when the nullspace dimension is k, k calls are likely to produce a nullspace basis. This function is important because nullspace bases can be so large and expensive to compute. In fact, a single integer matrix null space vector or solution vector often occupies more storage than the input matrix.
  • nullspaceBasisRight(D, A): The columns of dense matrix D are the nullspace basis. Generically this can be done with k calls to nullspaceVectorRight, where k is the nullspace dimension. However, both elimination and blackbox methods can be more efficient when asked for a full basis at once.
  • solveRightRandom(x, A, b): same as solveRight, but a random solution vector is promised. This can be more efficient than adding a solveRight and a nullspaceVectorRight.
  • solveRightMoorePenrose(x, A, b): returns least squares fit of least norm when IntegerTag, When ModularTag, returns x = Mb, when the Moore-Penrose inverse M exists and throws an error when not. As with solve, when the optional method is present, the failure is marked in the method object rather than an error thrown.
  • solveRightNonsingular(): I'd prefer not to provide this and offer the user the opportunity to promise nonsingularity (and gain some speed) by calling solveRight(x, A, b, Method::Hybrid("nonsingular")). If someone can show me the nonsingularity must used at compile time...

This handling of random and arbitrary solutions conforms to Chen, et al.

  1. JGD: Generic Lifting code revamping --> generic with respect to a "SolveModp?" class 2a Blas Based would do:
    • Precompute QLUP mod p using FFPACK
    • Give dense Q, L, U, P to a BlasSolve? class which provide a solve mod p method using FFPACK triangul solve
    • pass this solve class to the generic lifting
    2b Sparse elim would do:
    • Precompute QLUP mod p using GaussDomain?
    • Give sparse Q, L, U, P to a SparseSolve? class which provide a solve mod p method using GaussDomain? triangul solve
    • pass this solve class to the generic lifting
    2c BB Lanczos would do:
    • no precomputation
    • Give A to a LanczosSolve? class which provide a solve mod p method
    • pass this solve class to the generic lifting
    2d BB Wiedemann would do:
    • Precompute minpoly
    • Give minpoly to a WiedSolve? class which provide a solve mod p method
    • pass this solve class to the generic lifting

etc.