/* Usage: mpiexec ex16 [-help] [all PETSc options] */ static char help[] = "Solves a sequence of linear systems with different right-hand-side vectors.\n\ Input parameters include:\n\ -ntimes : number of linear systems to solve\n\ -view_exact_sol : write exact solution vector to stdout\n\ -m : number of mesh points in x-direction\n\ -n : number of mesh points in y-direction\n\n"; /*T Concepts: KSP^repeatedly solving linear systems; Concepts: KSP^Laplacian, 2d Concepts: Laplacian, 2d Processors: n T*/ /* Include "petscksp.h" so that we can use KSP solvers. Note that this file automatically includes: petscsys.h - base PETSc routines petscvec.h - vectors petscmat.h - matrices petscis.h - index sets petscksp.h - Krylov subspace methods petscviewer.h - viewers petscpc.h - preconditioners */ #include int main(int argc,char **args) { Vec x,b,u; /* approx solution, RHS, exact solution */ Mat A; /* linear system matrix */ KSP ksp; /* linear solver context */ PetscReal norm; /* norm of solution error */ PetscErrorCode ierr; PetscInt ntimes,i,j,k,Ii,J,Istart,Iend; PetscInt m = 8,n = 7,its; PetscBool flg = PETSC_FALSE; PetscScalar v,one = 1.0,rhs; ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr); ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Compute the matrix for use in solving a series of linear systems of the form, A x_i = b_i, for i=1,2,... - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Create parallel matrix, specifying only its global dimensions. When using MatCreate(), the matrix format can be specified at runtime. Also, the parallel partitioning of the matrix is determined by PETSc at runtime. */ ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr); ierr = MatSetFromOptions(A);CHKERRQ(ierr); ierr = MatSetUp(A);CHKERRQ(ierr); /* Currently, all PETSc parallel matrix formats are partitioned by contiguous chunks of rows across the processors. Determine which rows of the matrix are locally owned. */ ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr); /* Set matrix elements for the 2-D, five-point stencil in parallel. - Each processor needs to insert only elements that it owns locally (but any non-local elements will be sent to the appropriate processor during matrix assembly). - Always specify global rows and columns of matrix entries. */ for (Ii=Istart; Ii0) {J = Ii - n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);} if (i0) {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);} if (j -pc_type -ksp_monitor -ksp_rtol These options will override those specified above as long as KSPSetFromOptions() is called _after_ any other customization routines. */ ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Solve several linear systems of the form A x_i = b_i I.e., we retain the same matrix (A) for all systems, but change the right-hand-side vector (b_i) at each step. In this case, we simply call KSPSolve() multiple times. The preconditioner setup operations (e.g., factorization for ILU) be done during the first call to KSPSolve() only; such operations will NOT be repeated for successive solves. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ntimes = 2; ierr = PetscOptionsGetInt(NULL,NULL,"-ntimes",&ntimes,NULL);CHKERRQ(ierr); for (k=1; k