Interfacing to numpy/scipy¶
In some occasions or for some users it might be interesting to access NGSolve data from python in a fashion which is compatible with numpy and/or scipy. We give a few examples of possible use cases.
Working with small vectors and dense matrices:¶
see [Vectors and matrices](ngspy-howto-linalg)
Working with large vectors¶
You can get a "view" on an NGSolve-BaseVector using .FV() (which will give you a FlatVector) combined with .NumPy() which will give a numpy array which operates on the NGSolve-Vector-Data. For example the following works, assuming b to be an NGSolve-Vector:
b.FV().NumPy()[:] = abs(b.FV().NumPy()) - 1.0
which will give you the component-wise operation (absolute value minus one) applied on the vector b. During this operation data does not need to be copied.
Working with sparse matrices¶
You can access the sparse matrix information of a BaseMatrix using
rows,cols,vals = a.mat.COO()
Note that a bilinear form with flag symmetric==True will only give you one half of the matrix. These information can be put into a standard scipy-matrix, e.g. with
import scipy.sparse as sp A = sp.csr_matrix((vals,(rows,cols)))
You can use this, for instance, to examine the sparsity pattern of your matrix:
import matplotlib.pylab as plt plt.spy(A) plt.show()
or to compute the condition number (note that we export to a dense matrix here):
import numpy as np np.linalg.cond(A.todense())
Using iterative solvers from scipy¶
To use iterative solvers from scipy we have to wrap a LinearOperator around the NGSolve-matrix. The crucial component is the application of matrix vector product. Here is a very simple example where the scipy-cg-solver is used to solve the linear system (no preconditioner, no Dirichlet-dofs):
import scipy import scipy.sparse.linalg tmp1 = f.vec.CreateVector() tmp2 = f.vec.CreateVector() def matvec(v): tmp1.FV().NumPy()[:] = v tmp2.data = a.mat * tmp1 return tmp2.FV().NumPy() A = scipy.sparse.linalg.LinearOperator( (a.mat.height,a.mat.width), matvec) u.vec.FV().NumPy()[:], succ = scipy.sparse.linalg.cg(A, f.vec.FV().NumPy())
You can also use a sparse matrix format from python to run to previous example, see above. However, for preconditioning actions a sparse matrix is not necessarily set up such that the LinearOperator is often more useful.