Lisandro Dalcin avatar Lisandro Dalcin committed 1501485

Add petsc4py slides

Comments (0)

Files changed (6)

 \titlepage
 \end{frame}
 
-\section{Introduction}
+\section*{Outline}
+\begin{frame}
+  \frametitle{Outline}
+  \tableofcontents
+\end{frame}
+
+\section{Overview}
 
 \begin{frame}\uncover<+->{}
   \frametitle{Motivation}
   \end{itemize}
 \end{frame}
 
+\subsection{Python}
+
 \begin{frame}\uncover<+->{}
   \frametitle{Why Python?}
   \begin{itemize}\itemsep=1.5ex
   \end{columns}
 \end{frame}
 
-\section{MPI for Python}
+\subsection{MPI}
 
-\begin{frame}
+\begin{frame}\uncover<+->{}
   \frametitle{What is MPI?}
+  \begin{block}{}
+    \centering
+    \emph{Message Passing Interface}\\
+    \url{http://www.mpi-forum.org}
+  \end{block}
   \begin{itemize}\itemsep=1.5ex
   \item Standardized \alert{message passing} system
-    \begin{itemize}\itemsep=1.5ex
+    \begin{itemize}
     \item platform-independent (POSIX, Windows)
-    \item many implementations and vendors\\
+    \item many implementations and vendors
       \begin{itemize}
       \item MPICH2, Open MPI
       \item HP, Intel, Oracle, Microsoft
       \end{itemize}
     \end{itemize}
   \item Specifies semantics of a set of \alert{library routines}
-    \begin{itemize}\itemsep=1.5ex
-    \item No especial compiler support
+    \begin{itemize}
+    \item No special compiler support
     \item Language-neutral (C, \Cpp, Fortran 90)
     \end{itemize}
   \item Standard versions (backward-compatible)
-    \begin{itemize}\itemsep=1.5ex
+    \begin{itemize}
     \item \alert{MPI-1} (1994, 2008)
     \item \alert{MPI-2} (1996, 2009)
     \item \alert{MPI-3} (under development)
   \end{itemize}
 \end{frame}
 
+\subsection{PETSc}
+
+\begin{frame}\uncover<+->{}
+  \frametitle{What is PETSc?}
+  \begin{block}{}
+    \centering
+    \emph{Portable, Extensible Toolkit for Scientific Computation}\\
+    \url{http://www.mcs.anl.gov/petsc}
+  \end{block}
+  \bigskip
+  \begin{itemize}\itemsep=1.5ex
+  \item PETSc is a suite of \alert<.>{algorithms} and
+    \alert<.>{data structures} \\for the numerical solution of
+    \begin{itemize}
+    \item problems in \alert{science and engineering}
+    \item based on \alert{partial differential equations} models
+    \item discretized with \alert{finite differences/volumes/elements}
+    \item leading to \alert{large scale applications}
+    \end{itemize}
+  \item PETSc employs the MPI standard for parallelism
+  \item PETSc has an \alert{OO design}, it is implemented in \alert<.>{C},\\
+    can be used from \alert<.>{\Cpp}, provides a \alert<.>{Fortran 90}
+    interface
+  \end{itemize}
+\end{frame}
+
+
+\section{MPI for Python}
+
 \begin{frame}
   \frametitle{MPI for Python (\textbf{mpi4py})}
   \begin{itemize}\itemsep=1.5ex
-  \item Full-featured Python bindings for
+  \item Python bindings for
     \href{http://www.mpi-forum.org}{\textbf{MPI}}
   \item API based on the standard MPI-2 \Cpp bindings
   \item Supports all MPI features
   \inputminted[firstline=3]{cython}{mpi4py-cython.pyx}
 \end{frame}
 
+\subsection{Features and examples}
+
 \begin{frame}
   \frametitle{\textbf{[mpi4py]} Features -- MPI-1}
   \begin{itemize}\itemsep=1.5ex
   \inputminted[firstline=1]{python}{pingpong_v2.py}
 \end{frame}
 
+\subsection{Performance}
+
 \begin{frame}
   \begin{center}
     Point to Point Throughput -- Gigabit Ethernet
   \includegraphics[scale=0.5]{PingPong_SM.pdf}
 \end{frame}
 
+\section{PETSc for Python}
 
+\begin{frame}
+  \frametitle{PETSc for Python (\textbf{petsc4py})}
+  \begin{itemize}\itemsep=1.5ex
+  \item Python bindings for \textbf{PETSc}.
+  \end{itemize}
+\end{frame}
+
+\subsection{Features and examples}
+
+\begin{frame}
+  \frametitle{\textbf{[petsc4py]} Features -- PETSc components}
+  \begin{itemize}\itemsep=1.5ex
+  \item \textbf{Index Sets}: permutations, indexing, renumbering.
+  \item \textbf{Vectors}: sequential and distributed.
+  \item \textbf{Matrices}: sequential and distributed, sparse and dense.
+  \item \textbf{Distributed Arrays}: regular grid-based problems.
+  \item \textbf{Linear Solvers}: Krylov subspace methods.
+  \item \textbf{Preconditioners}: sparse direct solvers, multigrid
+  \item \textbf{Nonlinear Solvers}: line search, trust region, matrix-free.
+  \item \textbf{Timesteppers}: time-dependent, linear and nonlinear PDE's.
+  \end{itemize}
+\end{frame}
+
+\begin{frame}
+  \frametitle{\textbf{[petsc4py]} Vectors (\texttt{Vec}) -- CG Method}
+  \begin{columns}[t]
+    \begin{column}{.45\textwidth}
+      \tiny
+      \begin{equation*}
+        \begin{split}
+          & cg(A,x,b,i_{max},\epsilon): \\
+          & \quad i \Leftarrow 0 \\
+          & \quad r \Leftarrow b - A x \\
+          & \quad d \Leftarrow r \\
+          & \quad \delta_{0} \Leftarrow r^T r \\
+          & \quad \delta_{   } \Leftarrow \delta_{0} \\
+          & \quad \text{while}\;\; i < i_{max} \text{ and } \\
+          & \quad\quad\qquad  \delta_{   } > \delta_{0} \epsilon^2 \text{ do} :\\
+          & \quad\quad\quad  q \Leftarrow Ad \\
+          & \quad\quad\quad  \alpha \Leftarrow \frac{\delta_{   }}{d^T q} \\
+          & \quad\quad\quad  x \Leftarrow x + \alpha d\\
+          & \quad\quad\quad  r \Leftarrow r - \alpha q\\
+          & \quad\quad\quad  \delta_{old} \Leftarrow \delta_{   } \\
+          & \quad\quad\quad  \delta_{   } \Leftarrow r^T r \\
+          & \quad\quad\quad  \beta \Leftarrow \frac{\delta_{   }}{\delta_{old}} \\
+          & \quad\quad\quad  d \Leftarrow r + \beta d\\
+          & \quad\quad\quad  i \Leftarrow i + 1
+        \end{split}
+      \end{equation*}
+    \end{column}
+    \begin{column}{.45\textwidth}
+      \tiny\inputminted{python}{petsc4py_vec.py}
+    \end{column}
+  \end{columns}
+\end{frame}
+
+\begin{frame}
+  \frametitle{\textbf{[petsc4py]} Matrices (\texttt{Mat}) [1]}
+  \small\inputminted[lastline=19]{python}{petsc4py_mat.py}
+\end{frame}
+
+\begin{frame}
+  \frametitle{\textbf{[petsc4py]} Matrices (\texttt{Mat}) [2]}
+  \small\inputminted[firstline=20]{python}{petsc4py_mat.py}
+\end{frame}
+
+\begin{frame}
+  \frametitle{\textbf{[petsc4py]} Linear Solvers (\texttt{KSP+PC})}
+  \small\inputminted{python}{petsc4py_ksp.py}
+\end{frame}
+
+\subsection{Interoperability}
+
+\begin{frame}
+  \frametitle{\textbf{[petsc4py]} Interoperability}
+  Support for wrapping other PETSc-based C/\Cpp/F90 codes
+  \bigskip
+  \begin{itemize}\itemsep=1.5ex
+  \item using \href{http://cython.org}{\textbf{Cython}}
+    (\texttt{cimport} statement)
+  \item using \href{http://swig.org}{\textbf{SWIG}}
+    (\textsl{typemaps} provided)
+  \item using \textbf{F2Py} (\texttt{fortran} attribute)
+  \end{itemize}
+\end{frame}
+
+\begin{frame}
+  \frametitle{\textbf{[petsc4py]} Interoperability -- \textbf{SWIG} }
+  \small
+  \inputminted{c}{petsc4py-swig.i}
+\end{frame}
+
+\begin{frame}
+  \frametitle{\textbf{[petsc4py]} Interoperability -- \textbf{F2Py} }
+  \scriptsize
+  \inputminted{fortran}{petsc4py-f2py.pyf}
+\end{frame}
+
+\section{Applications}
+
+\begin{frame}
+\end{frame}
+
+%\section{Closing}
 
 \end{document}
 

petsc4py-f2py.pyf

+python module MyPDE
+interface
+
+   subroutine FormInitGuess(dm, x, params, ierr)
+     integer,      intent(in)  :: dm, x
+     real(kind=8), intent(in)  :: params(3)
+     integer,      intent(out) :: ierr
+   end subroutine FormInitGuess
+
+   subroutine FormFunction(dm, x, F, params, ierr)
+     integer,      intent(in)  :: dm, x, F
+     real(kind=8), intent(in)  :: params(3)
+     integer,      intent(out) :: ierr
+   end subroutine FormFunction
+
+   subroutine FormJacobian(dm, x, J, params, ierr)
+     integer,      intent(in)  :: dm, x, J
+     real(kind=8), intent(in)  :: params(3)
+     integer,      intent(out) :: ierr
+   end subroutine FormJacobian
+   
+end interface
+end python module MyPDE
+%module MyPDE
+%include petsc4py/petsc4py.i
+%{
+#include "MyPDE.h"
+%}
+
+typedef struct Params {
+  double alpha;
+  double beta;
+  double gamma;
+} Params;
+
+PetscErrorCode FormInitGuess(DM dm, Vec x, Params *p);
+PetscErrorCode FormFunction(DM dm, Vec x, Vec F, Params *p);
+PetscErrorCode FormJacobian(DM dm, Vec x, Mat J, Params *p);
+# create linear solver,
+ksp = PETSc.KSP()
+ksp.create(PETSc.COMM_WORLD)
+
+# use conjugate gradients method
+ksp.setType('cg')
+# and incomplete Cholesky
+ksp.getPC().setType('icc')
+
+# obtain sol & rhs vectors
+x, b = A.getVecs()
+x.set(0)
+b.set(1)
+
+# and next solve
+ksp.setOperators(A)
+ksp.setFromOptions()
+ksp.solve(b, x)
+from petsc4py import PETSc
+
+# grid size and spacing
+m, n  = 32, 32
+hx = 1.0/(m-1)
+hy = 1.0/(n-1)
+
+# create sparse matrix
+A = PETSc.Mat()
+A.create(PETSc.COMM_WORLD)
+A.setSizes([m*n, m*n])
+A.setType('aij') # sparse
+
+# precompute values for setting
+# diagonal and non-diagonal entries
+diagv = 2.0/hx**2 + 2.0/hy**2
+offdx = -1.0/hx**2
+offdy = -1.0/hy**2
+
+# loop over owned block of rows on this
+# processor and insert entry values
+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
+
+# communicate off-processor values
+# and setup internal data structures
+# for performing parallel operations
+A.assemblyBegin()
+A.assemblyEnd()
+def cg(A, b, x, imax=50, eps=1e-6):
+    """
+    A, b, x  : matrix, rhs, solution
+    imax     : maximum iterations
+    eps      : relative tolerance
+    """
+    # allocate work vectors
+    r = b.duplicate()
+    d = b.duplicate()
+    q = b.duplicate()
+    # initialization
+    i = 0
+    A.mult(x, r)
+    r.aypx(-1, b)
+    r.copy(d)
+    delta_0 = r.dot(r)
+    delta = delta_0
+    # enter iteration loop
+    while (i < imax and 
+           delta > delta_0 * eps**2):
+        A.mult(d, q)
+        alpha = delta / d.dot(q)
+        x.axpy(+alpha, d)
+        r.axpy(-alpha, q)
+        delta_old = delta
+        delta = r.dot(r)
+        beta = delta / delta_old
+        d.aypx(beta, r)
+        i = i + 1
+    return i, delta**0.5
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.