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.

Initialization is similar to previous examples.

import sys
import 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
        gi = i // n  # map row number to
        gj = i - gi * n  # grid coordinates
        if gi > 0:
            j = i - n
            A[i, j] = offdx
        if gi < m - 1:
            j = i + n
            A[i, j] = offdx
        if gj > 0:
            j = i - 1
            A[i, j] = offdy
        if gj < 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(f'Number of iterations of the method: {its}')
    sol_type = E.getType()
    Print(f'Solution method: {sol_type}')
    nev, _ncv, _mpd = E.getDimensions()
    Print(f'Number of requested eigenvalues: {nev}')
    tol, maxit = E.getTolerances()
    Print(f'Stopping condition: tol={tol:.4g}, maxit={maxit}')
    nconv = E.getConverged()
    Print(f'Number of converged eigenpairs: {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(f' {k.real:9f}{k.imag:+9f} j  {error:12g}')
            else:
                Print(f' {k.real:12f}       {error:12g}')
        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(f'Symmetric-definite Eigenproblem, N={m * n} ({m}x{n} grid)')
    A = Laplacian2D(m, n)
    B = QuasiDiagonal(m * n)
    solve_eigensystem(A, B)


if __name__ == '__main__':
    main()