Factorize changes the input matrix

[germa89]

Trying to fix a unit test in the math module, I realized that the assert done on the m2 fail.

As @FredAns told me it seems that mm.factorize does change m2, which in my opinion doesn't make sense if we are returning an object (solver in this case).

Explaining it with code

dim = 1000

m2 = mm.rand(dim, dim)
# factorize do changes inplace in m2, so we
# need a copy to later compare.
m3 = m2.copy()

solver = mm.factorize(m2)

v = mm.ones(dim)
C = solver.solve(v)

assert np.allclose(C, np.linalg.solve(m2.asarray(), v.asarray())) # Assertion Error
assert np.allclose(C, np.linalg.solve(m3.asarray(), v.asarray())) # True

Current implementation

As we saw in the example above, the method factorize return an object Solver which we can use it to solve the linear system (A*x=B).
However, this method also change the matrix m2, meaning that if you try to use m2 to solve a linear system using Numpy for example, you won't get the right answer (That is why the first assertion fail).

I believe this is not reasonable, if we use m2 as an argument of a function (mm.factorize) and it does return an object, to me the changes in m2 are not obvious.

Solutions

1. If we want to keep the memory consumption to minimum, I would probably implement the factorize feature as a method of m2 like:

m2.factorize()

I think this way, it is more obvious that the changes are affecting m2.

It might be an issue changing object nature on m2 because solving in implemented on AnsSolver class, and m2 is an AnsMat, but both are derived from ApdlMathObj.

2. Using copy. The factorize method will internally copy the m2 matrix, and perform factorization on this copy. Then you could do something like this:

solver = mm.factorize(m2)
assert np.allclose(C, np.linalg.solve(m2.asarray(), v.asarray())) # True

Assuming the object Solver is intrinsically different (not only different in the class definition) than the factorized version of m2, it will generate a copy of m2 which will be deleted once the function has finished. Hence this is not very performant compared with the above one. However, this is the type of notation or way of working someone would expect in Numpy.

Pinging @akaszynski @koubaa @mikerife @FredAns

[FredAns]

Just copy/paste here some comments I've exchanged with German by chat:

If we do a copy we need to warn the user. It's really critical in a general case.
Copy a matrix in numpy or matlab is something we can accept, because usually we handle matrices of size 100 x 100 or 1000 x 1000 max, for R&D purpose, and that's a few Mb max.
But APDLMath is an HPC Math Toolbox, and we work with matrices of size 100.000.000 x 100.000.000. User needs to know if we do a copy, we cannot hide that, especially on the cloud, where it's easy to not know what mem is used and how much free memory remains. As I said copy such a matrix can be a problem and limit what you can do after.
For a static linear analysis, we don't care about the K matrix after factorization, we just build it to be factorized. So we need to be able to not copy it.

and the m2.factorize() and m2.solve() design is not really what we tend to do, at least on a daily C++ context, where we prefer to specialize small objects instead of having big objects with everything inside. That's a bit orthogonal with the componentization goal.

a 10.000 x 10.000 dense matrix is already 800Mb

Memory allocations is one of the most critical point we see with big customers (who are running very big models). We try to save every Gb of memory we can. These solves can require several Tb of RAM these days.

[FredAns]

I would vote, as German suggested to me, to add an option to the factorize method, to do a copy of the matrix is required, so if copy=yes this copy of the orig matrix is modified by the factorization call, but keeping the default behavior to 'not copy'
And I agree we need to enhance the doc to warn about this overwrite of data in the input matrix ( or to add a warning when we call factorize)

[koubaa]

Copying is something the user should opt into, as it is one of the most expensive operations you can do. I do not think an option to copy is necessary, the user should be able to copy it themselves before calling factorize if that is what they want. The function should document clearly what it does, so there's no surprise there.

@FredAns it is not clear to me what the solver return value intended for.

[FredAns]

@koubaa we can have several factorized operators in the same time, so we need to have one object to handle it.

[FredAns]

dim = 1000
m2 = mm.rand(dim, dim)
engine = mm.DSP(mtx=m2) # __init__.py documents that it will mutate m2 and own it.
engine.solve()

This design hard codes the list of available engine types in the python code.
By using a string we keep the possibility to access new types available in Mapdl without the need to change the python Client API

[koubaa]

@FredAns so how about:

dim = 1000
m2 = mm.rand(dim, dim)
engine = mm.Engine(name="DSP", mtx=m2)
engine.solve()

[FredAns]

Yep I agree with this design, assuming you just forgot the args of the solve call:

dim = 1000
m2 = mm.rand(dim, dim)
engine = mm.Engine(name="DSP", mtx=m2)
engine.solve( B, X)

Pehaps it could be good to put the factorization step out of the Engine creation, so the user master this call ( which the most time consuming in this process):

dim = 1000
m2 = mm.rand(dim, dim)
engine = mm.Engine(name="DSP", mtx=m2)
engine.factorize()
engine.solve( B, X)

So the code of the Engine function can be inspired from the existing solver function ( for instance the way we select 'LAPACK' or 'DSP' Engine type wrt the Matrix type ( Dense or Sparse) when no engine name is specified.

[germa89]

The current implementation is:

dim = 1000
m2 = mm.rand(dim, dim)
engine = mm.AnsSolver(name="DSP")
engine.factorize(m2,)
engine.solve( B, X)

I don't see we are proposing any change (just Engine for AnsSolver), so I will keep the implementation as it is.

[germa89]

CHANGED THE DISCUSSION BODY

[germa89]

I would personally go for both proposed solutions; implement factorize in the array object and use a copy when using mm.factorize (we can implement an option for do the changes in place so we keep the current behavior, but the default should be a copy).

I believe the current implementation is prone to confusion.

[mikerife]

Hi @germa89 it makes more sense to me to keep the behavior; unless we have really good documentation that highlights the difference in this and the straight APDL. In either case I'm no fan of a default copy being made. I'm not sure of the 1 TB sized models that @FredAns is talking about (that's a big model!). But ~100 GB is pretty standard these days for just the run-of-the-mill medium sized models. See those all the time.
Mike

[germa89]

Ok, Default will be no copy then.

[germa89]

It seems there is consensus on keep the default behavior, add a copy option and improve docs.

[koubaa]

I personally don't feel a copy option is necessary

[germa89]

I kind of agree... most of the people won't use the copy (if you know what are you doing).

But having that option/argument also reinforce the idea that the user should be careful about the matrix being changed in factorize, because this behavior is somehow more present, in the function signature and the docs.

As I said, it is not needed, but I don't think it is totally useless.

[koubaa]

@germa89 I agree it could be useful, but necessary is a higher bar

[germa89]

I'm OK with that. :)