Skip to content

Allow re-use of analysis phase / expose _call_pardiso as non-private call_pardiso #88

@Huite

Description

@Huite

Thanks for wrapping PARDISO and making it available!

Issue

I'm solving a weakly nonlinear problem with Picard iteration where the sparsity pattern stays constant but matrix coefficients change. Currently, pypardiso uses is_already_factorized to choose between phase 33 (solve only) or phase 13 (analysis + factorization + solve), but I need phase 23 (factorization + solve, skipping analysis).

Performance impact: For my problem (phase 13 = 4.4s total):

  • Analysis: ~3.0s (70%)
  • Factorization: ~1.2s
  • Solve: ~0.2s

Skipping the repeated analysis step would significantly speed up iteration.

Current Workaround

I can work around this issue pretty easily by calling the private _call_pardiso method directly:

def picard(A, b, x_initial, maxiter):
    solver = pypardiso.PyPardisoSolver()
    solver.set_phase(11)
    solver._call_pardiso(A, b)
     
    for i in range(maxiter):
        update(A, b, x)
        solver.set_phase(23)
        x = solver._call_pardiso(A, b)
        if converged(x, b):
            return i, x
    
    return maxiter, x

This makes me slightly uncomfortable since I'm guessing _call_pardiso is private, although the remove_stored_factorization method mentions it expliclitly.

Suggested Solutions

  1. Add an analyze() method (similar to existing factorize()) that stores the sparsity pattern separately
  2. Make _call_pardiso public as it's not too error-prone with proper documentation

The second option seems most straightfoward, also because remove_stored_factorization mentions the method (suggesting it should be public).
Adding analyze has the difficulty of analyzing the sparsity pattern; checking indices and indptr is insufficient since nonstructural zeros may be inserted without showing up in the indices.
(Side note: the current hashing can be optimized by using the short-circuiting nature of Python's and -- I'll suggest this is in another issue since #1 mentioned performance.)

Would any of these approaches work for official support? (Or both.) Happy to submit a PR.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions