Use PETSC and SLEPC package to solve large scale linear system
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
ex1: ./configure PETSC_ARCH=intel-64-complex -with-cc=icc --with-fc=ifort --with-debugging=0 -download-c-blas-lapack=1 --download-mpich=1 --withscalar-type=complex
6
PETSC & SLEPC
2011/07/29
Ex2: ./configure PETSC_ARCH=intel-64-O3-real --withcc=gcc --with-fc=gfortran --with-debugging=1 --with-blaslapack-dir=/opt/intel/mkl/10.2.5.035 --with-mpidir=/opt/mpich2 --with-scalar-type=real --downloadmumps=1
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)
InsertMode : INSERT_VALUES, ADD_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
addv)
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
2. Read matrix from file
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
-${CLINKER} –o linearsym linearsym.o{PETSC_LIB}
${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
-${CLINKER} -o eigensym
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/
PETSC home page :
http://www.mcs.anl.gov/petsc/petsc-as/
user guide:
http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsccurrent/docs/manual.pdf
SLEPC home page :
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
© Copyright 2024