Return wrong results

[jingxu94]

Hi @haasad, thank you for this project.
I've been using PyPardiso recently and found that sometimes the result $x$ is wrong. I checked carefully and found that there may be a bug in the memory usage of the right-hand side of the equation. My code loops through each column of a 2D array as the rhs, at which point pypardiso gets the wrong solution. Here is a simple code snippet that reproduces the problem, where $x_1$ is the wrong result ($A \times x_1 \neq b_1, A \times x_2 = b_2 = b_1 $).

ndim = 10
A = sp.rand(ndim, ndim, density=0.5, format='csr')
b = np.random.rand(ndim, ndim)
b1 = b[:, 0]
b2 = b[:, 0].copy()
x1 = pypardiso.spsolve(A, b1)
x2 = pypardiso.spsolve(A, b2)

An instance of Class <ndarray> is composed of a contiguous one-dimensional memory segment and an indexing scheme. I guess the current code needs the rhs to point a contiguous memory slice, because as code snippet above, using the copy method or take each row of the two-dimensional array $b$ as the rhs can get the correct result.

[haasad]

Hi @jingxu94,

thanks for the report, I'll look into it.

[haasad]

Hi @jingxu94,

thanks again for the detailed bug report. You were absolutely right; there was a bug related to the memory layout for one-dimensional slices:
7998452

The pardiso solver expects fortran contiguous memory layout for the rhs and the memory layout wasn't changed in the case of one-dimensional slices.

I fixed it and released pypardiso version 0.4.1 on pypi and conda-forge.

[jingxu94]

Hi @haasad,

Nice work! Thanks again for developing and maintaining this project.

PyPardiso allows me to focus on other aspects of my python code without worrying about the efficiency of solving equations with large sparse matrices.

[haasad]

Happy to hear to pypardiso is useful to you!

As a side note:

My code loops through each column of a 2D array as the rhs, at which point pypardiso gets the wrong solution.

Are you aware that you can pass to whole 2D array to the solver and the solutions will be the columns of the returned array? This would be faster than looping.

[jingxu94]

Thank you for tips. In fact, I'm building a matrix where each column depends on the solution of the equation with the previous column as the rhs. This is an iterative process with no parallelism.

[haasad]

Ah yes, that makes sense 👍 just wanted to be sure