haasad / PyPardiso

Python interface to the Intel MKL Pardiso library to solve large sparse linear systems of equations
BSD 3-Clause "New" or "Revised" License
135 stars 20 forks source link

save factorized A to npz file #34

Closed xper0418 closed 2 years ago

xper0418 commented 2 years ago

First of all, thank you for this project. Very useful for me. Is there any way to save factorized A to scipy sparse npz or any other files, so I can load files any time and solve matrices?

haasad commented 2 years ago

Hi @xper0418, the factorization is not available as objects accessible by python, but is stored internally by the Pardiso solver. So storing the factorization from python won't be possible.

However there appears to be a function pardiso-handle-store, which can store and restore the state of the solver to a file, but I have never tried it. I think a wrapper for these functions could be added to pypardiso without too much effort.

What is your use case? I'd be interested to see the Pardiso statistics for the problem(s) you're solving. You can turn the statistical info of the solver on by running pypardiso.ps.set_statistical_info_on() before making a call to pypardiso.spsolve.

xper0418 commented 2 years ago

Thank you for your quick response. I've tried pypardiso.set_statistical_info_on(), but it looks like it doesn't work on windows, and I also ran the code on pycharm not notebook. A matrix is stiffness matrix from finite element analysis, and B matrix is load matrix in my case.

haasad commented 2 years ago

It is pypardiso.ps.set_statistical_info_on(), not pypardiso.set_statistical_info_on(). Could you maybe try that again?

In the statistical info it is shown how long the solver spends in each phase. This would help to estimate the time you could save by storing the factorization to disk.

haasad commented 2 years ago

Example ipython session:

In [1]: import pypardiso

In [2]: import numpy as np

In [3]: import scipy.sparse as sp

In [4]: pypardiso.ps.set_statistical_info_on()

In [5]: A = sp.rand(100, 100, density=0.5, format='csr')

In [6]: b = np.random.rand(100)

In [7]: x = pypardiso.spsolve(A, b)
=== PARDISO is running in In-Core mode, because iparam(60)=0 ===

Percentage of computed non-zeros for LL^T factorization
 1 %  2 %  3 %  4 %  100 % 

=== PARDISO: solving a real nonsymmetric system ===
1-based array indexing is turned ON
PARDISO double precision computation is turned ON
METIS algorithm at reorder step is turned ON
Scaling is turned ON
Matching is turned ON
Single-level factorization algorithm is turned ON

Summary: ( starting phase is reordering, ending phase is factorization )
================

Times:
======
Time spent in calculations of symmetric matrix portrait (fulladj): 0.000161 s
Time spent in reordering of the initial matrix (reorder)         : 0.002396 s
Time spent in symbolic factorization (symbfct)                   : 0.001044 s
Time spent in data preparations for factorization (parlist)      : 0.000006 s
Time spent in copying matrix to internal data structure (A to LU): 0.000000 s
Time spent in factorization step (numfct)                        : 0.001085 s
Time spent in allocation of internal data structures (malloc)    : 0.000176 s
Time spent in matching/scaling                                   : 0.000396 s
Time spent in additional calculations                            : 0.000367 s
Total time spent                                                 : 0.005631 s

Statistics:
===========
Parallel Direct Factorization is running on 4 OpenMP

< Linear system Ax = b >
             number of equations:           100
             number of non-zeros in A:      5000
             number of non-zeros in A (%): 50.000000

             number of right-hand sides:    1

< Factors L and U >
             number of columns for each panel: 128
             number of independent subgraphs:  0
< Preprocessing with state of the art partitioning metis>
             number of supernodes:                    6
             size of largest supernode:               95
             number of non-zeros in L:                9415
             number of non-zeros in U:                385
             number of non-zeros in L+U:              9800
             gflop   for the numerical factorization: 0.000628

             gflop/s for the numerical factorization: 0.578536

=== PARDISO: solving a real nonsymmetric system ===

Summary: ( solution phase )
================

Times:
======
Time spent in direct solver at solve step (solve)                : 0.000183 s
Time spent in additional calculations                            : 0.000033 s
Total time spent                                                 : 0.000216 s

Statistics:
===========
Parallel Direct Factorization is running on 4 OpenMP

< Linear system Ax = b >
             number of equations:           100
             number of non-zeros in A:      5000
             number of non-zeros in A (%): 50.000000

             number of right-hand sides:    1

< Factors L and U >
             number of columns for each panel: 128
             number of independent subgraphs:  0
< Preprocessing with state of the art partitioning metis>
             number of supernodes:                    6
             size of largest supernode:               95
             number of non-zeros in L:                9415
             number of non-zeros in U:                385
             number of non-zeros in L+U:              9800
             gflop   for the numerical factorization: 0.000628

             gflop/s for the numerical factorization: 0.578536

In [8]: 
xper0418 commented 2 years ago

I'm sorry, ctrl+c, v is not allowed at this place, so I can show you only part of statistics.

Time spent in reordering ~ : 11.3s Time spent in factorization step : 22.3s Total time spent : 39.8s

Statistics: Parallel Direct Factorization is running on 16 OpenMP

number of equations : 1234698 number of non-zeros in A : 49453830 number of non-zeros in A(%) : 0.003244

number of non-zeros in L : 711564777 number of non-zeros in U : 682468307 number of non-zeros in L+U : 1394033084

Summary: Time spent in direct solver at solve step : 2.7s