Use of base Solver class

Can you add a base Solver class for later?

[dpo]

That should probably be in a separate repository.

Yes, in POCS or similar.

I'm trying my own version of this generic base solver class:

https://github.com/PythonOptimizers/SuiteSparse.py/blob/develop/suitesparse/solver.cpx

It's quite basic but it already gives you an idea:

• already index/element typed;
• constructor takes only "one" fixed argument: the matrix to solve itself. All other arguments are named arguments. This is on purpose and is a feature!
• the reference to the matrix is self.__A.
• solve, analyze and factorize are completely generic.

What do you think?

@dpo @syarra

[dpo]

Sounds good. I'm thinking the class should keep track of whether analyze and factorize have been performed. If analyze hasn't been called and the user calls factorize, then it should call analyze first automatically. Similarly, solve should call factorize if it hasn't been called already.

The sophisticated user will want to call analyze/factorize separately to look at the ordering or input their own, or modify the factorization if it's permitted.

It seems those features could be implemented generically.

Another simple addition could be to implement __call__ and __mul__, as aliases to solve.

Yes for all except __mul__?

Something like this?

A = ...
solver = Solver(A)
# what would this code do?
solver * b

[dpo]

Yes, except I would write

LU = LUFactorization(A)  # for example
x = LU * b  # x solves the system Ax=b

OK, makes sense.

Callbacks to the rescue!

Here is an example:

from suitesparse.umfpack.umfpack_solver import UmfpackSolver
from cysparse.sparse.ll_mat import *

A = LLSparseMatrix(ncol=3, nrow=4)

s = UmfpackSolver(A)

print s.solver_name

b = 3

x = s * b

print "call factorize for the first time"
s.factorize()

print "=" * 80

print "recall factorize without effect..."
s.factorize()

print "=" * 80
print "recall factorize with force_factorize..."
s.factorize(force_factorize=True)

print "=" * 80
print "recall factorize with force_analyze..."
s.factorize(force_analyze=True)

print "=" * 80
print "recall factorize with force_analyze and force_factorize..."
s.factorize(force_analyze=True, force_factorize=True)

Output:

UMFPACK
Calling real solve with  3
call factorize for the first time
Called real analyze
Called real factorize
================================================================================
recall factorize without effect...
================================================================================
recall factorize with force_factorize...
Called real factorize
================================================================================
recall factorize with force_analyze...
Called real analyze
================================================================================
recall factorize with force_analyze and force_factorize...
Called real analyze
Called real factorize

See magic here

By the way, just to let you know: CySparse can now be completely imported with one "single" import:

from cysparse.sparse.ll_mat import *

No need to import anything else. And yes it is in the doc, and no you are not penalized.

[syarra]

I like it!
The only thing missing is that A might not be a CySparse matrix.
It could be a sparse matrix in the triplet format (arow, acol, aval).
For instance, in qr_mumps.py the triplet format is the base class and a derived class is made for CySparse matrices.
What do you think about it?

No problem. Will fix this. ;-)

By the way, A is a Python object... could be a CySparse matrix or a NumPy matrix or whatever...

[dpo]

It looks good, but force_analyze=True together with force_factorize=True is redundant. You can't factorize without analyzing first. So the syntax should just be

s.factorize(force=True)

and that would perform the analysis followed by the numerical factorization.

I was also thinking of something like

class DirectSolver(object):
  def __init__(self, ...):
    self.__analyzed = False
    self.__factorized = False
    ...

  def analyze(self, force=False):
    if self.__analyzed and not force:
      return
    self.__analyze() # perform analysis
    self.__analyzed = True

  def factorize(self, force=False):
    if self.__factorized and not force:
      return
    self.analyze(force=force)
    self.__factorize() # perform factorization
    self.__factorized = True