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
|
# ex5.py: Simple quadratic eigenvalue problem
# ===========================================
#
# This example solves a polynomial eigenvalue problem of degree 2,
# which is the most commonly found in applications. Larger degree
# polynomials are handled similarly. The coefficient matrices in
# this example do not come from an application, they are just simple
# matrices that are easy to build, just to illustrate how it works.
# The eigenvalues in this example are purely imaginary and come in
# conjugate pairs.
#
# The full source code for this demo can be `downloaded here
# <../_static/ex5.py>`__.
# Initialization is similar to previous examples.
import sys, slepc4py
slepc4py.init(sys.argv)
from petsc4py import PETSc
from slepc4py import SLEPc
Print = PETSc.Sys.Print
# A function to build the matrices. The lowest degree coefficient is
# the 2-D Laplacian, the highest degree one is the identity matrix,
# and the other matrix is set to zero (which means that this problem
# could have been solved as a linear eigenproblem).
def construct_operators(m,n):
Print("Quadratic Eigenproblem, N=%d (%dx%d grid)"% (m*n, m, n))
# K is the 2-D Laplacian
K = PETSc.Mat().create()
K.setSizes([n*m, n*m])
K.setFromOptions()
Istart, Iend = K.getOwnershipRange()
for I in range(Istart,Iend):
v = -1.0; i = I//n; j = I-i*n;
if i>0:
J=I-n; K[I,J] = v
if i<m-1:
J=I+n; K[I,J] = v
if j>0:
J=I-1; K[I,J] = v
if j<n-1:
J=I+1; K[I,J] = v
v=4.0; K[I,I] = v
K.assemble()
# C is the zero matrix
C = PETSc.Mat().create()
C.setSizes([n*m, n*m])
C.setFromOptions()
C.assemble()
# M is the identity matrix
M = PETSc.Mat().createConstantDiagonal([n*m, n*m], 1.0)
return M, C, K
# The polynomial eigenvalue solver is similar to the linear eigensolver
# used in previous examples. The main difference is that we must provide
# a list of matrices, from lowest to highest degree.
def solve_eigensystem(M, C, K):
# Setup the eigensolver
Q = SLEPc.PEP().create()
Q.setOperators([K, C, M])
Q.setDimensions(6)
Q.setProblemType(SLEPc.PEP.ProblemType.GENERAL)
Q.setFromOptions()
# Solve the eigensystem
Q.solve()
# Create the result vectors
xr, xi = K.createVecs()
its = Q.getIterationNumber()
Print("Number of iterations of the method: %i" % its)
sol_type = Q.getType()
Print("Solution method: %s" % sol_type)
nev, ncv, mpd = Q.getDimensions()
Print("")
Print("Number of requested eigenvalues: %i" % nev)
tol, maxit = Q.getTolerances()
Print("Stopping condition: tol=%.4g, maxit=%d" % (tol, maxit))
nconv = Q.getConverged()
Print("Number of converged approximate eigenpairs: %d" % nconv)
if nconv > 0:
Print("")
Print(" k ||(k^2M+Ck+K)x||/||kx|| ")
Print("-------------------- -------------------------")
for i in range(nconv):
k = Q.getEigenpair(i, xr, xi)
error = Q.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 two user-defined command-line options
# (the dimensions of the mesh) and calls the other two functions.
if __name__ == '__main__':
opts = PETSc.Options()
m = opts.getInt('m', 32)
n = opts.getInt('n', m)
M, C, K = construct_operators(m,n)
solve_eigensystem(M, C, K)
M = C = K = None
|
