ex6.py: Compute exp(t*A)*v for a matrix from a Markov model#

This example illustrates the functionality in slepc4py for computing matrix functions, or more precisely, the application of a matrix function on a given vector. The example works with the exponential function, which is most commonly found in applications.

The main focus of slepc4py is eigenvalue and singular value problems, but it has some codes to deal with matrix functions, which sometimes are needed in the context of eigenproblems, but have interest on their own.

The full source code for this demo can be downloaded here.

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

This function builds a matrix that implements a Markov model of a random walk on a triangular grid. The entries of the matrix represent probabilities of moving to neighboring cells in the grid.

def build_matrix(m):
    N = m*(m+1)/2
    Print("Markov y=exp(t*A)*e_1, N=%d (m=%d)"% (N, m))
    A = PETSc.Mat().create()
    A.setSizes([N, N])
    A.setFromOptions()
    Istart, Iend = A.getOwnershipRange()
    ix = 0
    cst = 0.5/(m-1)
    for i in range(1,m+1):
       jmax = m-i+1
       for j in range(1,jmax+1):
           ix = ix + 1
           if ix-1<Istart or ix>Iend:
               continue  # compute only owned rows
           if j!=jmax:
               pd = cst*(i+j-1)
               # north
               if i==1:
                   A[ix-1,ix] = 2*pd
               else:
                   A[ix-1,ix] = pd
               # east
               if j==1:
                   A[ix-1,ix+jmax-1] = 2*pd
               else:
                   A[ix-1,ix+jmax-1] = pd
           # south
           pu = 0.5 - cst*(i+j-3)
           if j>1:
               A[ix-1,ix-2] = pu
           # west
           if i>1:
               A[ix-1,ix-jmax-2] = pu

    A.assemble()
    return A

The following function solves the problem. This case is quite different from eigenproblems, and is more similar to solving a linear system of equations with KSP. To configure the problem we must provide the matrix and the function (the exponential in this case). Note how the internal FN object is extracted from the MFN solver. Also, it is often necessary to specify a scale factor, which in this case represents the time for which we want to obtain the evolved state. Once the solver is set up, we call solve() passing the right-hand side vector b and the solution vector x.

def solve_exp(t, A, b, x):
    # Setup the solver
    M = SLEPc.MFN().create()
    M.setOperator(A)
    f = M.getFN()
    f.setType(SLEPc.FN.Type.EXP)
    f.setScale(t)
    M.setTolerances(1e-7)
    M.setFromOptions()
    # Solve the problem
    M.solve(b,x)

    its = M.getIterationNumber()
    Print("Number of iterations of the method: %i" % its)
    sol_type = M.getType()
    Print("Solution method: %s" % sol_type)
    ncv = M.getDimensions()
    Print("")
    Print("Subspace dimension: %i" % ncv)
    tol, maxit = M.getTolerances()
    Print("Stopping condition: tol=%.4g, maxit=%d" % (tol, maxit))
    Print("Computed vector at time t=%.4g has norm %g" % (t.real, x.norm()))
    Print("")

The main program processes the command-line option m (size of the grid), builds the matrix and calls the solver. Note how the vectors are created from the matrix. In this case, the right-hand side vector is the first element of the canonical basis.

if __name__ == '__main__':
    opts = PETSc.Options()
    m = opts.getInt('m', 15)
    A = build_matrix(m)   # transition probability matrix
    x, b = A.createVecs()
    x.set(0)
    b.set(0)
    b[0] = 1
    b.assemble()
    t = 2
    solve_exp(t, A, b, x)    # compute x=exp(t*A)*b
    A = None
    b = x = None

This Page