static char help[] = "Block Jacobi preconditioner for solving a linear system in parallel with KSP.\n\ The code indicates the\n\ procedures for setting the particular block sizes and for using different\n\ linear solvers on the individual blocks.\n\n"; /* Note: This example focuses on ways to customize the block Jacobi preconditioner. See ex1.c and ex2.c for more detailed comments on the basic usage of KSP (including working with matrices and vectors). Recall: The block Jacobi method is equivalent to the ASM preconditioner with zero overlap. */ /*T Concepts: KSP^customizing the block Jacobi preconditioner 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; /* KSP context */ KSP *subksp; /* array of local KSP contexts on this processor */ PC pc; /* PC context */ PC subpc; /* PC context for subdomain */ PetscReal norm; /* norm of solution error */ PetscErrorCode ierr; PetscInt i,j,Ii,J,*blks,m = 4,n; PetscMPIInt rank,size; PetscInt its,nlocal,first,Istart,Iend; PetscScalar v,one = 1.0,none = -1.0; PetscBool isbjacobi; ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr); ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); n = m+2; /* ------------------------------------------------------------------- Compute the matrix and right-hand-side vector that define the linear system, Ax = b. ------------------------------------------------------------------- */ /* Create and assemble parallel matrix */ 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 = MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);CHKERRQ(ierr); ierr = MatSeqAIJSetPreallocation(A,5,NULL);CHKERRQ(ierr); ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr); for (Ii=Istart; Ii0) {J = Ii - n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} if (i0) {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} if (j */ ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); ierr = PCSetType(pc,PCBJACOBI);CHKERRQ(ierr); /* ------------------------------------------------------------------- Define the problem decomposition ------------------------------------------------------------------- */ /* Call PCBJacobiSetTotalBlocks() to set individually the size of each block in the preconditioner. This could also be done with the runtime option -pc_bjacobi_blocks Also, see the command PCBJacobiSetLocalBlocks() to set the local blocks. Note: The default decomposition is 1 block per processor. */ ierr = PetscMalloc1(m,&blks);CHKERRQ(ierr); for (i=0; i -sub_ksp_type -sub_ksp_rtol 1.e-4 */ /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Advanced method, setting different solvers for various blocks. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Note that each block's KSP context is completely independent of the others, and the full range of uniprocessor KSP options is available for each block. The following section of code is intended to be a simple illustration of setting different linear solvers for the individual blocks. These choices are obviously not recommended for solving this particular problem. */ ierr = PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&isbjacobi);CHKERRQ(ierr); if (isbjacobi) { /* Call KSPSetUp() to set the block Jacobi data structures (including creation of an internal KSP context for each block). Note: KSPSetUp() MUST be called before PCBJacobiGetSubKSP(). */ ierr = KSPSetUp(ksp);CHKERRQ(ierr); /* Extract the array of KSP contexts for the local blocks */ ierr = PCBJacobiGetSubKSP(pc,&nlocal,&first,&subksp);CHKERRQ(ierr); /* Loop over the local blocks, setting various KSP options for each block. */ for (i=0; i