Data Objects

1 / 30

# Data Objects - PowerPoint PPT Presentation

##### Data Objects

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - - -
##### Presentation Transcript

1. Data Objects • Vectors (Vec) • focus: field data arising in nonlinear PDEs • Matrices (Mat) • focus: linear operators arising in nonlinear PDEs (i.e., Jacobians) • Object creation • Object assembly • Setting options • Viewing • User-defined customizations beginner beginner intermediate intermediate advanced tutorial outline: data objects

2. Vectors proc 0 • Fundamental objects for storing field solutions, right-hand sides, etc. • VecCreateMPI(...,Vec *) • MPI_Comm - processors that share the vector • number of elements local to this processor • total number of elements • Each process locally owns a subvector of contiguously numbered global indices proc 1 proc 2 proc 3 proc 4 data objects: vectors beginner

3. Vectors: Example #include "vec.h" int main(int argc,char **argv) { Vec x,y,w; /* vectors */ Vec *z; /* array of vectors */ double norm,v,v1,v2;int n = 20,ierr; Scalar one = 1.0,two = 2.0,three = 3.0,dots[3],dot; PetscInitialize(&argc,&argv,PETSC_NULL,PETSC_NULL); ierr = OptionsGetInt(PETSC_NULL,"n",&n,PETSC_NULL); CHKERRA(ierr); ierr =VecCreate(PETSC_COMM_WORLD,PETSC_DECIDE,n,&x); CHKERRA(ierr); ierr = VecSetFromOptions(x);CHKERRA(ierr); ierr = VecDuplicate(x,&y);CHKERRA(ierr); ierr = VecDuplicate(x,&w);CHKERRA(ierr); ierr = VecSet(&one,x);CHKERRA(ierr); ierr = VecSet(&two,y);CHKERRA(ierr); ierr = VecSet(&one,z[0]);CHKERRA(ierr); ierr = VecSet(&two,z[1]);CHKERRA(ierr); ierr = VecSet(&three,z[2]);CHKERRA(ierr); ierr = VecDot(x,x,&dot);CHKERRA(ierr); ierr = VecMDot(3,x,z,dots);CHKERRA(ierr);

4. Vectors: Example ierr = PetscPrintf(PETSC_COMM_WORLD,"Vector length %d\n",(int)dot); CHKERRA(ierr); ierr = VecScale(&two,x);CHKERRA(ierr); ierr = VecNorm(x,NORM_2,&norm);CHKERRA(ierr); v = norm-2.0*sqrt((double)n); if (v > -1.e-10 && v < 1.e-10) v = 0.0; ierr = PetscPrintf(PETSC_COMM_WORLD,"VecScale %g\n",v); CHKERRA(ierr); ierr = VecDestroy(x);CHKERRA(ierr); ierr = VecDestroy(y);CHKERRA(ierr); ierr = VecDestroy(w);CHKERRA(ierr); ierr = VecDestroyVecs(z,3);CHKERRA(ierr); PetscFinalize(); return 0; }

5. Vector Assembly • VecSetValues(Vec,…) • number of entries to insert/add • indices of entries • values to add • mode: [INSERT_VALUES,ADD_VALUES] • VecAssemblyBegin(Vec) • VecAssemblyEnd(Vec) data objects: vectors beginner

6. Parallel Matrix and Vector Assembly • Processors may generate any entries in vectors and matrices • Entries need not be generated on the processor on which they ultimately will be stored • PETSc automatically moves data during the assembly process if necessary data objects: vectors and matrices beginner

7. Selected Vector Operations

8. Sparse Matrices • Fundamental objects for storing linear operators (e.g., Jacobians) • MatCreateMPIAIJ(…,Mat *) • MPI_Comm - processors that share the matrix • number of local rows and columns • number of global rows and columns • optional storage pre-allocation information data objects: matrices beginner

9. proc 1 proc 2 } proc 3 proc 3: locally owned rows proc 4 Parallel Matrix Distribution • Each process locally owns a submatrix of contiguously numbered global rows. MatGetOwnershipRange(Mat A, int *rstart, int *rend) • rstart: first locally owned row of global matrix • rend -1: last locally owned row of global matrix proc 0 data objects: matrices beginner

10. Matrix Assembly • MatSetValues(Mat,…) • number of rows to insert/add • indices of rows and columns • number of columns to insert/add • values to add • mode: [INSERT_VALUES,ADD_VALUES] • MatAssemblyBegin(Mat) • MatAssemblyEnd(Mat) data objects: matrices beginner

11. Viewer Concepts • Information about PETSc objects • runtime choices for solvers, nonzero info for matrices, etc. • Data for later use in restarts or external tools • vector fields, matrix contents • various formats (ASCII, binary) • Visualization • simple x-window graphics • vector fields • matrix sparsity structure beginner viewers

12. Viewing Vector Fields Solution components, using runtime option -snes_vecmonitor • VecView(Vec x,Viewer v); • Default viewers • ASCII (sequential): VIEWER_STDOUT_SELF • ASCII (parallel): VIEWER_STDOUT_WORLD • X-windows: VIEWER_DRAW_WORLD • Default ASCII formats • VIEWER_FORMAT_ASCII_DEFAULT • VIEWER_FORMAT_ASCII_MATLAB • VIEWER_FORMAT_ASCII_COMMON • VIEWER_FORMAT_ASCII_INFO • etc. velocity: v velocity: u temperature: T vorticity: z beginner viewers

13. Viewing Matrix Data • MatView(Mat A, Viewer v); • Runtime options available after matrix assembly • -mat_view_info • info about matrix assembly • -mat_view_draw • sparsity structure • -mat_view • data in ASCII • etc. beginner viewers

14. Linear Solvers Goal: Support the solution of linear systems, Ax=b, particularly for sparse, parallel problems arising within PDE-based models User provides: • Code to evaluate A, b solvers:linear beginner

15. Linear Solver Object SLES Object KSP Object Coeff. Matrix - Method: CG, GMRES etc. - Convergence tolerance - ... Precond. Matrix PC object - Pctype: ILU, LU etc.

16. Basic Linear Solver Code (C/C++) SLES sles; /* linear solver context */ Mat A; /* matrix */ Vec x, b; /* solution, RHS vectors */ int n, its; /* problem dimension, number of iterations */ MatCreate(MPI_COMM_WORLD,n,n,&A); /* assemble matrix */ VecCreate(MPI_COMM_WORLD,n,&x); VecDuplicate(x,&b); /* assemble RHS vector */ SLESCreate(MPI_COMM_WORLD,&sles); SLESSetOperators(sles,A,A,DIFFERENT_NONZERO_PATTERN); SLESSetFromOptions(sles); SLESSolve(sles,b,x,&its); solvers:linear beginner

17. Basic Linear Solver Code (Fortran) SLES sles Mat A Vec x, b integer n, its, ierr call MatCreate(MPI_COMM_WORLD,n,n,A,ierr) call VecCreate(MPI_COMM_WORLD,n,x,ierr) call VecDuplicate(x,b,ierr) call SLESCreate(MPI_COMM_WORLD,sles,ierr) call SLESSetOperators(sles,A,A,DIFFERENT_NONZERO_PATTERN,ierr) call SLESSetFromOptions(sles,ierr) call SLESSolve(sles,b,x,its,ierr) C then assemble matrix and right-hand-side vector solvers:linear beginner

18. Customization Options • Procedural Interface • Provides a great deal of control on a usage-by-usage basis inside a single code • Gives full flexibility inside an application • Command Line Interface • Applies same rule to all queries via a database • Enables the user to have complete control at runtime, with no extra coding solvers:linear intermediate

19. Setting Solver Options within Code • SLESGetKSP(SLES sles,KSP *ksp) • KSPSetType(KSP ksp,KSPType type) • KSPSetTolerances(KSP ksp,double rtol,double atol,double dtol, int maxits) • ... • SLESGetPC(SLES sles,PC *pc) • PCSetType(PC pc,PCType) • PCASMSetOverlap(PC pc,int overlap) • ... solvers:linear intermediate

20. 2 1 beginner intermediate Setting Solver Options at Runtime • -ksp_type [cg,gmres,bcgs,tfqmr,…] • -pc_type [lu,ilu,jacobi,sor,asm,…] • -ksp_max_it <max_iters> • -ksp_gmres_restart <restart> • -pc_asm_overlap <overlap> • -pc_asm_type [basic,restrict,interpolate,none] • etc ... 1 2 solvers:linear

21. SLES: Review of Basic Usage • SLESCreate( ) - Create SLES context • SLESSetOperators( ) - Set linear operators • SLESSetFromOptions( ) - Set runtime solver options for [SLES,KSP,PC] • SLESSolve( ) - Run linear solver • SLESView( ) - View solver options actually used at runtime (alternative: -sles_view) • SLESDestroy( ) - Destroy solver solvers:linear beginner

22. 2 1 beginner intermediate SLES: Review of Selected Preconditioner Options 1 2 And many more options... solvers: linear: preconditioners

23. 2 1 beginner intermediate SLES: Review of Selected Krylov Method Options 1 2 And many more options... solvers: linear: Krylov methods

24. SLES: Runtime Script Example solvers:linear intermediate

25. Viewing SLES Runtime Options solvers:linear intermediate

26. 2 1 3 beginner intermediate advanced SLES: Example Programs Location:www.mcs.anl.gov/petsc/src/sles/examples/tutorials/ • ex1.c, ex1f.F - basic uniprocessor codes • ex2.c, ex2f.F - basic parallel codes • ex11.c - using complex numbers • ex4.c - using different linear system and preconditioner matrices • ex9.c - repeatedly solving different linear systems • ex15.c - setting a user-defined preconditioner 1 2 3 • And many more examples ... solvers:linear

27. Example: ex1.c (Sequential Tridiagonal Solver) int main(int argc,char **args) { Vec x, b, u; /* approx solution, RHS, exact solution */ Mat A; /* linear system matrix */ SLES sles; /* linear solver context */ PC pc; /* preconditioner context */ KSP ksp; /* Krylov subspace method context */ double norm; /* norm of solution error */ int ierr, i, n = 10, col[3], its, flg, size; Scalar neg_one = -1.0, one = 1.0, value[3]; PetscInitialize(&argc,&args,(char *)0,help); MPI_Comm_size(PETSC_COMM_WORLD,&size); if (size != 1) SETERRA(1,0,"This is a uniprocessor example only!"); ierr = OptionsGetInt(PETSC_NULL,"-n",&n,&flg); CHKERRA(ierr); ierr = MatCreate(PETSC_COMM_WORLD,n,n,&A); CHKERRA(ierr); value[0] = -1.0; value[1] = 2.0; value[2] = -1.0; for (i=1; i<n-1; i++ ) { col[0] = i-1; col[1] = i; col[2] = i+1; ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES); CHKERRA(ierr); } i = n - 1; col[0] = n - 2; col[1] = n - 1; ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES); CHKERRA(ierr); i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0; ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES); CHKERRA(ierr); ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRA(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRA(ierr);

28. Example: ex1.c (cont’d) ierr = VecCreate(PETSC_COMM_WORLD,PETSC_DECIDE,n,&x); ierr = VecSetFromOptions(x);CHKERRA(ierr); ierr = VecDuplicate(x,&b); CHKERRA(ierr); ierr = VecDuplicate(x,&u); CHKERRA(ierr); ierr = VecSet(&one,u); CHKERRA(ierr); ierr = MatMult(A,u,b); CHKERRA(ierr); ierr = SLESCreate(PETSC_COMM_WORLD,&sles); CHKERRA(ierr); ierr =SLESSetOperators(sles,A,A,DIFFERENT_NONZERO_PATTERN); ierr = SLESGetKSP(sles,&ksp); CHKERRA(ierr); ierr = SLESGetPC(sles,&pc); CHKERRA(ierr); ierr = PCSetType(pc,PCJACOBI); CHKERRA(ierr); ierr = KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT, PETSC_DEFAULT); CHKERRA(ierr); /* Set runtime options, e.g., -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol> */ ierr = SLESSetFromOptions(sles); CHKERRA(ierr); ierr = SLESSolve(sles,b,x,&its); CHKERRA(ierr);

29. Example: ex1.c (cont’d) ierr = SLESSolve(sles,b,x,&its); CHKERRA(ierr); ierr = SLESView(sles,VIEWER_STDOUT_WORLD); CHKERRA(ierr); ierr = VecAXPY(&neg_one,u,x); CHKERRA(ierr); ierr = VecNorm(x,NORM_2,&norm); CHKERRA(ierr); if (norm > 1.e-12) PetscPrintf(PETSC_COMM_WORLD,"Norm= %g, Iters= %d\n",norm,its); else PetscPrintf(PETSC_COMM_WORLD,"Norm < 1.e-12, Iters= %d\n",its); ierr = VecDestroy(x); CHKERRA(ierr); ierr = VecDestroy(u); CHKERRA(ierr); ierr = VecDestroy(b); CHKERRA(ierr); ierr = MatDestroy(A);CHKERRA(ierr); ierr = SLESDestroy(sles); CHKERRA(ierr); PetscFinalize(); return 0; }

30. Example: ex2.c (Parallel Laplace Solver) int main(int argc,char **args) { Vec x, b, u; /* approx solution, RHS, exact solution */ Mat A; /* linear system matrix */ SLES sles; /* linear solver context */ PetscRandom rctx; /* random number generator context */ double norm; /* norm of solution error */ int i, j, I, J, Istart, Iend, ierr, m = 8, n = 7, its, flg; Scalar v, one = 1.0, neg_one = -1.0; KSP ksp; PetscInitialize(&argc,&args,(char *)0,help); ierr = OptionsGetInt(PETSC_NULL,"-m",&m,&flg); CHKERRA(ierr); ierr = OptionsGetInt(PETSC_NULL,"-n",&n,&flg); CHKERRA(ierr); ierr = MatGetOwnershipRange(A,&Istart,&Iend); CHKERRA(ierr); for ( I=Istart; I<Iend; I++ ) { v = -1.0; i = I/n; j = I - i*n; if ( i>0 ) {J = I - n; MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES); } if ( i<m-1 ) {J = I + n; MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES); } if ( j>0 ) {J = I - 1; MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES); } if ( j<n-1 ) {J = I + 1; MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES); } v = 4.0; MatSetValues(A,1,&I,1,&I,&v,INSERT_VALUES); } ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRA(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRA(ierr);