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
|
# ------------------------------------------------------------------------
# Standard symmetric eigenproblem for the Laplacian operator in 1-D
# ------------------------------------------------------------------------
import sys, slepc4py
slepc4py.init(sys.argv)
from petsc4py import PETSc
from slepc4py import SLEPc
import numpy
opts = PETSc.Options()
n = opts.getInt('n', 30)
A = PETSc.Mat(); A.create()
A.setSizes([n, n])
A.setFromOptions()
rstart, rend = A.getOwnershipRange()
# first row
if rstart == 0:
A[0, :2] = [2, -1]
rstart += 1
# last row
if rend == n:
A[n-1, -2:] = [-1, 2]
rend -= 1
# other rows
for i in range(rstart, rend):
A[i, i-1:i+2] = [-1, 2, -1]
A.assemble()
E = SLEPc.EPS(); E.create()
E.setOperators(A)
E.setProblemType(SLEPc.EPS.ProblemType.HEP)
history = []
def monitor(eps, its, nconv, eig, err):
if nconv<len(err): history.append(err[nconv])
E.setMonitor(monitor)
E.setFromOptions()
E.solve()
Print = PETSc.Sys.Print
Print()
Print("******************************")
Print("*** SLEPc Solution Results ***")
Print("******************************")
Print()
its = E.getIterationNumber()
Print( "Number of iterations of the method: %d" % its )
eps_type = E.getType()
Print( "Solution method: %s" % eps_type )
nev, ncv, mpd = E.getDimensions()
Print( "Number of requested eigenvalues: %d" % 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:
# Create the results vectors
v, _ = A.createVecs()
#
Print()
Print(" k ||Ax-kx||/||kx|| ")
Print("----------------- ------------------")
for i in range(nconv):
k = E.getEigenpair(i, v)
error = E.computeError(i)
Print( " %12f %12g" % (k, error) )
Print()
|
