1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
|
ex3.py: Matrix-free eigenproblem for the 2-D Laplacian
======================================================
This example solves the eigenproblem for the 2-D discrete Laplacian
without building the matrix explicitly.
The full source code for this demo can be `downloaded here
<../_static/ex3.py>`__.
Initialization is similar to previous examples.
::
import sys, slepc4py
slepc4py.init(sys.argv)
from petsc4py import PETSc
from slepc4py import SLEPc
import numpy as np
In this case, the program cannot be run in parallel, so we check that
the number of MPI processes is 1. In order to enable parallelism, we
should implement a parallel matrix-vector operation ourselves, which
is not done in this example.
::
assert PETSc.COMM_WORLD.getSize() == 1
Print = PETSc.Sys.Print
This function computes the matrix-vector product f = L*x where the
Laplacian L is not built explicitly, and the vectors x,f are viewed
as two-dimensional arrays associated to grid points.
::
def laplace2d(U, x, f):
U[:,:] = 0
U[1:-1, 1:-1] = x
# Grid spacing
m, n = x.shape
hx = 1.0/(m-1) # x grid spacing
hy = 1.0/(n-1) # y grid spacing
# Setup 5-points stencil
u = U[1:-1, 1:-1] # center
uN = U[1:-1, :-2] # north
uS = U[1:-1, 2: ] # south
uW = U[ :-2, 1:-1] # west
uE = U[2:, 1:-1] # east
# Apply Laplacian
f[:,:] = \
(2*u - uE - uW) * (hy/hx) \
+ (2*u - uN - uS) * (hx/hy) \
For a matrix-free solution in slepc4py we have to create a class that
wraps the matrix-vector operation and optionally other operations of
the matrix. In this case, we provide the constructor and the ``mult``
operation, that simply calls the ``laplace2d`` function above.
::
class Laplacian2D(object):
def __init__(self, m, n):
self.m, self.n = m, n
scalar = PETSc.ScalarType
self.U = np.zeros([m+2, n+2], dtype=scalar)
def mult(self, A, x, y):
m, n = self.m, self.n
xx = x.getArray(readonly=1).reshape(m,n)
yy = y.getArray(readonly=0).reshape(m,n)
laplace2d(self.U, xx, yy)
In this example, building the matrix amounts to creating an object of
the class defined above, and passing it to a special petsc4py matrix
with `createPython() <petsc4py.PETSc.Mat.createPython>`.
::
def construct_operator(m, n):
# Create shell matrix
context = Laplacian2D(m,n)
A = PETSc.Mat().createPython([m*n,m*n], context)
return A
This function receives the matrix and the problem type, then solves the
eigenvalue problem and prints information about the computed solution.
Although we know that eigenvalues and eigenvectors are real in this
example, the function is prepared to solve it as a non-symmetric problem,
by passing `SLEPc.EPS.ProblemType.NHEP`, that is why the code handles
possibly complex eigenvalues and eigenvectors.
::
def solve_eigensystem(A, problem_type=SLEPc.EPS.ProblemType.HEP):
# Create the result vectors
xr, xi = A.createVecs()
# Setup the eigensolver
E = SLEPc.EPS().create()
E.setOperators(A,None)
E.setDimensions(3,PETSc.DECIDE)
E.setProblemType(problem_type)
E.setFromOptions()
# Solve the eigensystem
E.solve()
Print("")
its = E.getIterationNumber()
Print("Number of iterations of the method: %i" % its)
sol_type = E.getType()
Print("Solution method: %s" % sol_type)
nev, ncv, mpd = E.getDimensions()
Print("Number of requested eigenvalues: %i" % nev)
tol, maxit = E.getTolerances()
Print("Stopping condition: tol=%.4g, maxit=%d" % (tol, maxit))
nconv = E.getConverged()
Print("Number of converged eigenpairs: %d" % nconv)
if nconv > 0:
Print("")
Print(" k ||Ax-kx||/||kx|| ")
Print("----------------- ------------------")
for i in range(nconv):
k = E.getEigenpair(i, xr, xi)
error = E.computeError(i)
if k.imag != 0.0:
Print(" %9f%+9f j %12g" % (k.real, k.imag, error))
else:
Print(" %12f %12g" % (k.real, error))
Print("")
The main program simply processes three user-defined command-line options
and calls the other two functions.
::
def main():
opts = PETSc.Options()
N = opts.getInt('N', 32)
m = opts.getInt('m', N)
n = opts.getInt('n', m)
Print("Symmetric Eigenproblem (matrix-free), "
"N=%d (%dx%d grid)" % (m*n, m, n))
A = construct_operator(m,n)
solve_eigensystem(A)
if __name__ == '__main__':
main()
|
