1 / 30

Data Objects

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.

amanda-cherry
Download Presentation

Data Objects

An Image/Link below is provided (as is) to download presentation 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. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

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);

More Related