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
151
152
153
154
155
156
|
ex9.py: Generalized symmetric-definite eigenproblem
===================================================
This example computes eigenvalues and eigenvectors of a generalized
symmetric-definite eigenvalue problem, where the first matrix is the
discrete Laplacian in two dimensions and the second matrix is quasi
diagonal.
The full source code for this demo can be `downloaded here
<../_static/ex9.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
This function builds the discretized Laplacian operator in 2 dimensions.
::
def Laplacian2D(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 builds a quasi-diagonal matrix. It is two times the identity
matrix except for the 2x2 leading submatrix ``[6 -1; -1 1]``.
::
def QuasiDiagonal(N):
# Create matrix
B = PETSc.Mat().create()
B.setSizes([N, N])
B.setFromOptions()
# Fill matrix
Istart, Iend = B.getOwnershipRange()
for I in range(Istart, Iend):
B[I,I] = 2.0
if Istart==0:
B[0,0] = 6.0
B[0,1] = -1.0
B[1,0] = -1.0
B[1,1] = 1.0
B.assemble()
return B
The following function receives the two matrices and solves the
eigenproblem. In this example we illustrate how to pass objects
that have been created beforehand, instead of extracting the internal
objects. We are using a spectral transformation of type `ST.Type.PRECOND`
and a Block Jacobi preconditioner. We want to compute the leftmost
eigenvalues. The selected eigensolver is LOBPCG, which is appropriate
for this use case. After the solve, we print the computed solution.
::
def solve_eigensystem(A, B, problem_type=SLEPc.EPS.ProblemType.GHEP):
# Create the results vectors
xr, xi = A.createVecs()
pc = PETSc.PC().create()
# pc.setType(pc.Type.HYPRE)
pc.setType(pc.Type.BJACOBI)
ksp = PETSc.KSP().create()
ksp.setType(ksp.Type.PREONLY)
ksp.setPC( pc )
F = SLEPc.ST().create()
F.setType(F.Type.PRECOND)
F.setKSP( ksp )
F.setShift(0)
# Setup the eigensolver
E = SLEPc.EPS().create()
E.setST(F)
E.setOperators(A,B)
E.setType(E.Type.LOBPCG)
E.setDimensions(10,PETSc.DECIDE)
E.setWhichEigenpairs(E.Which.SMALLEST_REAL)
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 functions.
::
def main():
opts = PETSc.Options()
N = opts.getInt('N', 10)
m = opts.getInt('m', N)
n = opts.getInt('n', m)
Print("Symmetric-definite Eigenproblem, N=%d (%dx%d grid)" % (m*n, m, n))
A = Laplacian2D(m,n)
B = QuasiDiagonal(m*n)
solve_eigensystem(A,B)
if __name__ == '__main__':
main()
|
