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
|
ex2.py: Standard symmetric eigenproblem for the 2-D Laplacian
=============================================================
This example computes eigenvalues and eigenvectors of the discrete Laplacian
on a two-dimensional domain with finite differences.
The full source code for this demo can be `downloaded here
<../_static/ex2.py>`__.
Initialization is similar to previous examples.
::
try: range = xrange
except: pass
import sys, slepc4py
slepc4py.init(sys.argv)
from petsc4py import PETSc
from slepc4py import SLEPc
Print = PETSc.Sys.Print
In this example we have organized the code in several functions. This
one builds the finite-difference Laplacian matrix by computing the
indices of each entry. An alternative would be to use the functionality
offered by `DMDA <petsc4py.PETSc.DMDA>`.
::
def construct_operator(m, n):
# Create matrix for 2D Laplacian operator
A = PETSc.Mat().create()
A.setSizes([m*n, m*n])
A.setFromOptions()
# Fill matrix
hx = 1.0/(m-1) # x grid spacing
hy = 1.0/(n-1) # y grid spacing
diagv = 2.0*hy/hx + 2.0*hx/hy
offdx = -1.0*hy/hx
offdy = -1.0*hx/hy
Istart, Iend = A.getOwnershipRange()
for I in range(Istart, Iend) :
A[I,I] = diagv
i = I//n # map row number to
j = I - i*n # grid coordinates
if i> 0 : J = I-n; A[I,J] = offdx
if i< m-1: J = I+n; A[I,J] = offdx
if j> 0 : J = I-1; A[I,J] = offdy
if j< n-1: J = I+1; A[I,J] = offdy
A.assemble()
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 (sparse matrix), "
"N=%d (%dx%d grid)" % (m*n, m, n))
A = construct_operator(m,n)
solve_eigensystem(A)
if __name__ == '__main__':
main()
|
