|
#!/usr/bin/env python3 |
|
""" poisson_nd( n, d ) -> sparse matrix `A` 1d n x n, 2d n^2 x n^2, 3d n^3 x n^3 |
|
poisson2( 3 ): |
|
[[ 4 -1 . -1 . . . . .] |
|
[-1 4 -1 . -1 . . . .] |
|
[ . -1 4 . . -1 . . .] |
|
|
|
[-1 . . 4 -1 . -1 . .] |
|
[ . -1 . -1 4 -1 . -1 .] |
|
[ . . -1 . -1 4 . . -1] |
|
|
|
[ . . . -1 . . 4 -1 .] |
|
[ . . . . -1 . -1 4 -1] |
|
[ . . . . . -1 . -1 4]] |
|
(A_h = this / h^2) |
|
|
|
poisson2( n ) - 4 I has n null vecs, |
|
poisson2( n ) - 4.0001 I has n near-null vecs -- a nice test for arpack and *gmres |
|
inv( poisson1( n ) - 2 I ) is all 1 0 -1 ! plot under my gists |
|
""" |
|
# see also Invop, Cheblinop |
|
|
|
from functools import partial |
|
import numpy as np |
|
from scipy import sparse |
|
|
|
#............................................................................... |
|
def poisson1( n ) -> "sparse csc tridiagonal [-1 2 1] ...": |
|
# Make Gilbert Strang's favorite matrix |
|
# http://www-math.mit.edu/~gs/PIX/cupcakematrix.jpg |
|
# w Discrete_Laplace_operator |
|
P1d = sparse.diags( [[-1]*(n-1), [2]*n, [-1]*(n-1)], |
|
[-1, 0, 1]) # ~ - f'' |
|
return P1d |
|
|
|
def poisson_nd( n, d=2 ): |
|
""" -> A sparse csc 1d n x n, 2d n^2 x n^2, 3d n^3 x n^3 """ |
|
# w Kronecker_sum_of_discrete_Laplacians |
|
# https://math.stackexchange.com/questions/2629649/understanding-high-and-low-frequency-eigenmodes |
|
# poisson2 evals: 2 (pi / n)^2 .. 8, cond ~ 0.8 n^2 |
|
n = int(n) |
|
assert n > 0, n |
|
P1d = poisson1(n) |
|
if d == 1: |
|
return P1d.tocsc() |
|
return sparse.kronsum( P1d, poisson_nd( n, d - 1 )).tocsc() |
|
# kron I A + kron B I |
|
|
|
poisson2 = partial( poisson_nd, d=2 ) |
|
poisson3 = partial( poisson_nd, d=3 ) |
|
poisson2.__name__ = "poisson2" |
|
poisson3.__name__ = "poisson3" |
|
|
|
__version__ = "2022-06-24 June" # denis-bz-py t-online.de |
|
|
|
#............................................................................... |
|
if __name__ == "__main__": # run various solvers, direct / iterative |
|
import sys |
|
from time import time |
|
from numpy.linalg import norm |
|
from scipy.sparse import linalg as sparselin |
|
|
|
print( 80 * "▄" ) |
|
# print( "▄ bp" ); breakpoint() |
|
try: |
|
from scikits import umfpack |
|
print( "umfpack:", umfpack.__version__ ) |
|
except ImportError: |
|
print( "umfpack:", None ) |
|
|
|
np.set_printoptions( edgeitems=10, threshold=10, linewidth=120, formatter = dict( |
|
float = lambda x: "%.2g" % x )) # float arrays %.2g |
|
|
|
def sparsesolver( name, A ): |
|
# https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html#solving-linear-problems |
|
if name == "dense": |
|
return lambda b: np.linalg.solve( A.A, b ) |
|
elif name == "splu": # direct, with umfpack v fast |
|
return sparselin.splu(A).solve |
|
elif name == "spsolve": # direct, v slow |
|
return lambda b: sparselin.spsolve( A, b ) |
|
else: |
|
solver = getattr( sparselin, name ) # iter lgmres ... sparsesolvers.py tl;dr |
|
return lambda b: solver( A, b, atol=1e-5 ) |
|
|
|
def nbytes( A ): |
|
return A.nbytes if hasattr( A, "nbytes" ) \ |
|
else A.data.nbytes + A.indices.nbytes + A.indptr.nbytes |
|
|
|
#........................................................................... |
|
solvers = "splu lgmres gmres" # spsolve" |
|
# solvers = "dense" |
|
d = 2 |
|
n = 100 |
|
incdiag = -4.0001 # -4: n null vecs of the n^2 |
|
save = 0 |
|
|
|
# to change these params, run this.py a=1 b=None 'c = expr' ... in sh or ipython -- |
|
for arg in sys.argv[1:]: |
|
exec( arg ) |
|
|
|
#........................................................................... |
|
A = poisson_nd( n, d=d ) # n^d x n^d |
|
N = A.shape[0] |
|
if incdiag != 0: |
|
A += incdiag * sparse.eye( N ) |
|
params = "poisson%d n %g incdiag %g A %s %.0f mbytes " % ( |
|
d, n, incdiag, A.shape, nbytes(A) / 1e6 ) |
|
print( "params:", params ) |
|
# print( "A: \n", A[:10,:10].A ) |
|
|
|
b = np.zeros( d * (n,) ) |
|
b[ d * (n // 2,) ] = 1 # midpoint |
|
b = b.reshape( -1 ) |
|
x0 = np.zeros( N ) # ? |
|
|
|
for solvername in solvers.split(): |
|
print( "\n--", solvername ) |
|
t0 = time() # wall clock, e.g. 329% of 4 cores |
|
#....................................................................... |
|
solver = sparsesolver( solvername, A ) # , b, x0 |
|
x = solver( b ) |
|
|
|
if isinstance( x, tuple ): |
|
x, info = x # iterative solvers |
|
b_Ax = norm( b - A.dot(x) ) / norm( b ) |
|
print( |
|
"time: %6.1f sec %-8s poisson%d n %d incdiag %g |b-Ax|/|b| %.2g x %g .. %g " % ( |
|
time() - t0, solvername, d, n, incdiag, |
|
b_Ax, x.min(), x.max() )) |
|
if d == 1: |
|
print( x ) |
|
elif d == 2 and n <= 20: |
|
print( x.reshape( n, n )) # incdiag: ripple |
|
|
|
if save: |
|
out = "poisson%d-n%d-%s.npz" % (d, n, solvername) |
|
print( "\n>", out, "\n" ) |
|
np.savez( out, x=x, n=n, params=params ) |
|
|
|
# $sslin/_dsolve/linsolve.py with umfpack -- |
|
# time: 1.0 sec splu poisson2 n 300 incdiag -4.0001 |
|
# time: 30.8 sec spsolve poisson2 n 300 incdiag -4.0001 |
|
# time: 6.7 sec splu poisson3 n 30 incdiag -6.0001 |
|
# time: 47.2 sec spsolve poisson3 n 30 incdiag -6.0001 |
|
|
