### petsc_sywang_2011

```Use PETSC and SLEPC package to solve
large scale linear system
---quick start
Sheng Yi Wang, Department of Mathematics, National Taiwan University
2011/07/29
Outline
PETSC
---Vector
---Matrix
---linear system solver
SLEPC
---eigenvalue problem solver


PBS

2
PETSC & SLEPC
2011/07/29
What is PETSC ??

PETSC = Portable Extensible Toolkit for Scientific
Computation

PETSC is intended for use in large-scale application
projects.

PETSC support MPI, that can parallel execute.

PETSC can interfaces many external software
ex: Matlab, Mathematica, MUMPS…etc
3
PETSC & SLEPC
2011/07/29
What can PETSC do??

Vector operation

Matrix operation

Linear system ( sparse or dense) Ax=b

Nonlinear solver

ODE & PDE solve
( steady state or time dependent)
4
PETSC & SLEPC
2011/07/29
User’s background

Some knowledge of parallel (MPI command)

You are familiar with C/Fortran language.

You are familiar with Linux environment.
PS: PETSC can be installed in Windows
(but you have to install many other packages)
5
PETSC & SLEPC
2011/07/29
How install??

Step1 : Set environment variables PETSC_DIR
ex: PETSC_DIR=/home/sywang/petsc-3.1-p8; export
PETSC_DIR

Step 2. Configure (to check your environment)
---BLAS, C compiler, (MPI) must be installed
6
PETSC & SLEPC
2011/07/29

7
PETSC & SLEPC
2011/07/29

Step 2: make all

Step 3: make test

Set PETSC_DIR and PETSC_ARCH to ~/.bashrc

If you install PETSC successfully, install SLEPC is very easy
8
PETSC & SLEPC
2011/07/29
Don’t be afraid, install the PETSC is the most difficult
process, after that , all you need to do is the following:
1. Call function
2. Set parameter
3. Show the output and explain the results

9
PETSC & SLEPC
2011/07/29
Start using PETSC
10
PETSC & SLEPC
2011/07/29
11
PETSC & SLEPC
2011/07/29
Declare Variables
Vec
Mat
KSP
PetscInt
PetscScalar
PetscScalar
PetscErrorCode
12
sol, rhs;
Mtx_A;
ksp;
ii= 10;
value[3];
val_rhs;
ierr;
PETSC & SLEPC
2011/07/29
Rough coding
PetscInitialize();
ObjCreate(MPI_comm,&obj);
ObjSetType(obj, );
ObjSetFromOptions(obj, );
ObjSolve(obj, );
ObjGetxxx(obj, );
ObjDestroy(obj);
PetscFinalize()
13
PETSC & SLEPC
2011/07/29
Vector
14
PETSC & SLEPC
2011/07/29
I. Vector

Vec a

VecCreate(MPI_Comm comm,Vec* x)
 VecCreate(PETSC_COMM_WORLD, &a)

VecSetSizes(Vec v, PetscInt n, PetscInt N)
 VecSetSizes(a,PETSC_DECIDE,20);
15
PETSC & SLEPC
2011/07/29
I. Vector

VecSet(Vec x,PetscScalar alpha)
 VecSet(a,1.0)

VecSetValues(Vec x,PetscInt ni,const PetscInt ix[],const
PetscScalar y[],InsertMode iora)
 VecSetValues(a,1,0,-3, INSERT_VALUES)
16
PETSC & SLEPC
2011/07/29
I. Vector

VecSetRandom(Vec x,PetscRandom rctx)
 VecSetRandom(b,r)

VecSetFromOptions(Vec vec)
 VecSetFromOptions(a)

VecDuplicate(Vec v,Vec *newv)
 VecDuplicate(a,&b)
17
PETSC & SLEPC
2011/07/29
I. Vector

VecView(a,PETSC_VIEWER_STDOUT_WORLD);


VecAssemblyBegin(a);
VecAssemblyEnd(a);

VecDestroy(a);
18
PETSC & SLEPC
2011/07/29
19
PETSC & SLEPC
2011/07/29
Example presentation

Example 1:
Create Vec and use some basic vector operation
20
PETSC & SLEPC
2011/07/29
Matrix
21
PETSC & SLEPC
2011/07/29
2. Matrix

Mat

MatCreate(MPI_Comm comm,Mat* A)
 MatCreate(PETSC_SOMM_WORLD,&mtx_a)

MatSetSizes(Mat A,int m,int n,int M,int N)
 MatSetSizes(mtx_a,PETSC_DECIDE,
PETSC_DECIDE,10,10)

MatSetFromOptions(Mat A)
 MatSetFromOptions(mtx_a)
22
mtx_a
PETSC & SLEPC
2011/07/29
2. Matrix

MatSetValues(Mat A,PetscInt m,const
PetscInt idxm[],PetscInt n,const PetscInt
idxn[],const PetscScalar v[],InsertMode
 MatSetValues(mtx_a,1,1,1,1,2.0,INSERT_VALUE)
23
PETSC & SLEPC
2011/07/29
Example

I have a small matrix s_a [1 2 ;3 4]
but I want to insert s_a into global matrix mtx_a and the
index is (1,2)
MatSetValues(mtx_a,2,1,2,2,s_a,INSERT_VALUE)

00000
00120
00340
00000
24
PETSC & SLEPC
2011/07/29
2. Matrix

MatAssemblyBegin(Mat, MAT_FINAL_ASSEMBLY)
MatAssemblyEnd(Mat ,MAT_FINAL_ASSEMBLY)

MatDestroy(Mat )

25
PETSC & SLEPC
2011/07/29
26
PETSC & SLEPC
2011/07/29
**2. Matrix: Matrix-Free

If you don’t want to create all matrix ( the matrix is too
large),PETSC allow users write their own matrix operator

extern int mult(Mat,Vec,Vec);
MatCreateShell(comm,m,n,M,N,ctx,&mat);
MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))
mult);
MatDestroy(mat);
27
PETSC & SLEPC
2011/07/29

In many case, someone has created the matrix, so we
don’t rewrite the matrix, we just want to read the matrix
and solve them quickly, how to do that??

If you have matrix which is ASCII format(just store the
nonzero element , i , j, aij)
we can call PETSC function to transfer ASCII to binary file
( PETSC can read them very quickly)
28
PETSC & SLEPC
2011/07/29
Binary file

PetscViewerBinaryOpen(MPI_Comm comm,const char
name[],PetscViewerFileType type,PetscViewer *binv)

29
PetscViewerBinaryOpen(PETSC_COMM_WORLD,fileout,
FILE_MODE_WRITE,&view)
PETSC & SLEPC
2011/07/29
Example presentation

Example 2:
Read matrix in ASCII format and transfer them to PETSC
binary file

Example 3:
Read matrix from PETSC binary file and create a vector
and use a basic matrix operation
30
PETSC & SLEPC
2011/07/29
Linear System (KSP)
31
PETSC & SLEPC
2011/07/29
3.Linear system solver : KSP

Linear system problem: give matrix A and vector b 
solve Ax=b

The dimension is too large to find inverse, and the matrix
is sparse, so we need to use iterative method to solve
them.

The basic method to solve linear system is Krylov
subspace methods.
32
PETSC & SLEPC
2011/07/29
Linear Solvers in PETSC
33
PETSC & SLEPC
2011/07/29
3. Linear system solver : KSP
First of all, set the matrix A and rhs vector b
Declare Variables
KSP
ksp

KSPCreate(MPI_Comm comm,KSP *ksp)
 KSPCreate(PETSC_COMM_WORLD,&ksp)

KSPSetOperators(KSP ksp,Mat Amat,Mat
Pmat,MatStructure flag)

34
KSPSetOperators(ksp,A,A, SAME_NONZERO_PATTERN)
PETSC & SLEPC
2011/07/29
3. Linear system solver : KSP

KSPSetType(KSP ksp, const KSPType type)
 KSPSetType(ksp,KSPGMRES)
-ksp_type

KSPSetFromOptions(KSP ksp)
 KSPSetFromOptions(ksp)

KSPSolve(ksp,Vec b,Vec x)
 KSPSolve(ksp,b,x)

KSPGetIterationNumber(KSP ksp,PetscInt *its)
 KSPGetIterationNumber(ksp,&it)
35
PETSC & SLEPC
2011/07/29
Preconditioner

In many case, the matrix A has large condition number, so
we need use a proconditoner to reduce condition number

PETSC provide many preconditioner, and some of them
can use MPI to parallel.
36
PETSC & SLEPC
2011/07/29
Precondioner in PETSC
37
PETSC & SLEPC
2011/07/29
3. Linear system solver : KSP

KSPSetFromOptions(KSP ksp)
 KSPSetFromOptions(ksp)

KSPGetPC(KSP ksp,PC *pc)
 KSPGetPC(ksp,&pc)

PCSetType(PC pc, PCType type)
 PCSetType(pc,PCBJACOBI) -pc_type
38
PETSC & SLEPC
2011/07/29
How to Parallel ??

MPI－ Message Passing Interface
PETSC and SLEPC support MPI and user don’t have to call
MPI function




mpiexec -np 4 yourprogram
mpirun -np 8 -machinefile mf yourprogram
PETSC_COMM_SELF PETSC_COMM_WORLD
PETSC and SLEPC initial will include mpi.h, so if you want
to use MPI function,you can use them, too
39
PETSC & SLEPC
2011/07/29
Assign parameters at run time

Serial :




./ex4 -ksp_type bcgs -pc_type lu -ksp_rtol 1e-4
./ex4 -ksp_type bcgs -pc_type lu -ksp_rtol 1e-4
./ex4 -ksp_type gmres -pc_type asm -ksp_rtol 1e-8 -ksp_max_it 20
Parallel:
In one node :

mpiexec –np 4 ./ex4 -ksp_type cg -pc_type sor - pc_sor_omega 1.8 ksp_rtol 1e-4 -ksp_view
In multi nodes:
 mpiexec –np 16 –machinefile mf ./ex4 -ksp_type cg -pc_type
asm -ksp_rtol 1e-4 -ksp_view -ksp_monitor_true_residual
40
PETSC & SLEPC
2011/07/29
Example presentation

Example 4 :
Read matrix from a binary file, and call PETSC
function to solve linear system. And show how
to assign parameters at run time
41
PETSC & SLEPC
2011/07/29
makefile
include \${PETSC_DIR}/conf/base
linearsym : linearsym.o chkopts
\${RM} linearsym.o
42
PETSC & SLEPC
2011/07/29
Bash script
#!/bin/bash
for((i=10;i<=30;i=i+5))
do
./linearsym -n \$i -ksp_type gmres -pc_type
jacobi
-ksp_max_it 100 -ksp_view > result_gmres_\$i
done
for((i=10;i<=30;i=i+5))
do
./linearsym -n \$i -ksp_type cg -pc_type jacobi
-ksp_max_it 100 -ksp_view > result_cg_\$i
done
43
PETSC & SLEPC
2011/07/29
Perl (bash script like)

#!/usr/bin/perl
for (\$i=1; \$i<=199; \$i++) {
\$sor_omega=0.01*\$i;
system("~/program/openmpi2/bin/mpiexec
-np \$i -machinefile mf10 ex4 -fin m22103 –
ksp_type cg –pc_type sor –pc_sor_omega
\$sor_omega -log_summary >
m22103_cg_sor_omega\$sor_omega");
system("echo \$i");
sleep(3);}
44
PETSC & SLEPC
2011/07/29
Eigenvalue (EPS)
45
PETSC & SLEPC
2011/07/29
4.Eigenvalue Problem Solver: SLEPC

SLEPC the Scalable Library for Eigenvalue Problem
Computations

Standard Eigenvalue problem :
give matrix A, want to find unknown vector x and value k
 Ax=kx

General Eigenvalue problem :
give matrix A and matrix B, want to find unknown vector
x and value k  Ax=kBx

Still , the matrix is very large and sparse.
46
PETSC & SLEPC
2011/07/29
How to install SLEPC

Export SLEPC_DIR =/home/sywang/petsc-3.1-p6;

./configure (they will follow the PETSC environmental
setting)

Make all

Make test
47
PETSC & SLEPC
2011/07/29
48
PETSC & SLEPC
2011/07/29
Declare Variables
EPS eps;
Mat A,B ;
PetscInt
ii,nn = 10,col[3];
PetscScalar value[3];
PetscScalar kr,ki;
PetscErrorCode ierr;
49
PETSC & SLEPC
2011/07/29
4.Eigenvalue Problem Solver: SLEPC

EPSCreate(MPI_Comm comm,EPS *eps)
 EPSCreate(PETSC_COMM_WORLD,eps)

EPSSetOperators(EPS eps,Mat A,Mat B)
 EPSSetOperators(eps,A,B)
 EPSSetOperators(eps,A,PETSC_NULL)

 Ax=kBx
 Ax=kx
EPSSetProblemType(EPS eps,EPSProblemType type)
 EPSSetProblemType(eps,EPS_HEP)
-eps_hermitian
 EPSSetProblemType(eps,EPS_NHEP)
eps_non_hermitian
50
PETSC & SLEPC
2011/07/29
4.Eigenvalue Problem Solver: SLEPC

EPSSetType(EPS eps,const EPSType type)
 EPSSetType(eps,EPSJD)
-eps_type jd

EPSSetTolerances(EPS eps,PetscReal tol,PetscInt maxits)
 EPSSetTolerances(eps,1e-8,200)
-eps_tol 1e-8 –eps_max_it 200
51
PETSC & SLEPC
2011/07/29
EPS Type
52
PETSC & SLEPC
2011/07/29
53
PETSC & SLEPC
2011/07/29
4.Eigenvalue Problem Solver: SLEPC

EPSSetDimensions(EPS eps,PetscInt nev,PetscInt ncv)
 EPSSetDimensions(eps,2,18) -eps_nev 2 -eps_ncv 18

EPSCreate(MPI_Comm comm,EPS *eps)
 EPSSetWhichEigenpairs(eps, EPS_SMALLEST_REAL)
-eps_smallest_real
54
PETSC & SLEPC
2011/07/29
EPS which
55
PETSC & SLEPC
2011/07/29
4.Eigenvalue Problem Solver: SLEPC

EPSSolve(eps)

EPSView(eps,PETSC_VIEWER_STDOUT_WORLD)
-eps_view
EPSGetIterationNumber(eps, &its)


EPSGetOperationCounters(EPS eps,PetscInt*
ops,PetscInt* dots,PetscInt* lits)
 EPSGetOperationCounters(eps,PETSC_NULL,PETSC_
NULL,&lits)
56
PETSC & SLEPC
2011/07/29
4.Eigenvalue Problem Solver: SLEPC

EPSGetType(eps,&type)

EPSGetDimensions(eps,&nev,&ncv)

EPSGetTolerances(eps,&tol,&maxit)

EPSGetConverged(eps,&nconv)
57
PETSC & SLEPC
2011/07/29
4.Eigenvalue Problem Solver: SLEPC

EPSGetEigenpair(EPS eps, PetscInt i, PetscScalar *eigr,
PetscScalar *eigi, Vec Vr,Vec Vi)
 EPSGetEigenpair(eps,i,&kr,&ki,PETSC_NULL,PETSC_NUL
L)

EPSComputeRelativeError(eps,i,&error)

EPSDestroy(eps)
58
PETSC & SLEPC
2011/07/29
Spectral Transformation
shift & invert

In many case, the eigenvalue we want is the smallet (but
nonzero),if we don’t use shift and invert, it takes many
time to find eigenvalue

But when you use invert, you will need to solve linear
system, they are done by calling PETSC KSP solver
59
PETSC & SLEPC
2011/07/29
Spectral Transformation

EPSCreate(MPI_Comm comm,EPS *eps)
 EPSCreate(PETSC_COMM_WORLD,eps)
ST
st;
PetscScalar shift = 2.0;

EPSGetST(EPS eps,ST* st);

STSetShift(ST st,PetscScalar shift);
 STSetShift(st,2.0);
-st_shift

STSetType(ST st,STType type);
 STSetType(st, STSHIFT) -st_type shift


60
PETSC & SLEPC
2011/07/29
61
PETSC & SLEPC
2011/07/29
Spectral Transformation

Shift-invert need to solve linear system. We need to call
PETSC function

In general we, set parameter in command line
ex:
-st_type sinvert -st_ksp_type cg st_pc_type asm -st_ksp_rtol 1e-10
-st_type sinvert -st_ksp_type bcgs st_pc_type sor -st_pc_sor_omega 1.8 eps_monitor_draw
62
PETSC & SLEPC
2011/07/29
makefile

all: eigensym
include
\${SLEPC_DIR}/conf/slepc_common
eigensym : eigensym.o chkopts
eigensym.o{SLEPC_LIB}
\${RM} eigensym.o
63
PETSC & SLEPC
2011/07/29
How to get info from PETSC and SLEPC

./linearsym -ksp_view

./eigensym -eps_view

./linearsym -ksp_monitor

./eigensym -eps_monitor

./eigensym -eps_monitor_draw_all

./linearsym -log_summary
64
PETSC & SLEPC
2011/07/29
Example

Example5 :
Read a matrix from a binary file, and find its
eigenvalue
65
PETSC & SLEPC
2011/07/29
PBS

PBS = Portable Batch System

PBS is first developed by NASA

You need to write a bash script and submit job to the
cluster
66
PETSC & SLEPC
2011/07/29
Script





#PBS
#PBS
#PBS
#PBS
-N
-q
-o
-e
Job_name
queue_name (according to cluster)
result_file
error_message
#PBS -l node=2:ppn=4
#PBS -l nodes=node01:ppn=4+node02:ppn=4
67
PETSC & SLEPC
2011/07/29
info

Submit job : qsub your_sheell

Check job

Check node : pbsmodes
pbsnodes –l free
68
: qstat
qstat -n
PETSC & SLEPC
2011/07/29
Reference




Hands-On Exercises for SLEPC
http://www.grycap.upv.es/slepc/handson/
http://www.mcs.anl.gov/petsc/petsc-as/
user guide:
http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsccurrent/docs/manual.pdf
http://www.grycap.upv.es/slepc/
user guide:
http://www.grycap.upv.es/slepc/documentation/slepc.pdf
Matrix market:
http://math.nist.gov/MatrixMarket/
69
PETSC & SLEPC
2011/07/29
Thank You
70
PETSC & SLEPC
2011/07/29
```