static char help[] = "Solves a tridiagonal linear system.\n\n"; /*T Concepts: KSP^basic parallel example; 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 Note: The corresponding uniprocessor example is ex1.c */ #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 */ PC pc; /* preconditioner context */ PetscReal norm,tol=1000.*PETSC_MACHINE_EPSILON; /* norm of solution error */ PetscErrorCode ierr; PetscInt i,n = 10,col[3],its,rstart,rend,nlocal; PetscScalar one = 1.0,value[3]; ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Compute the matrix and right-hand-side vector that define the linear system, Ax = b. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Create vectors. Note that we form 1 vector from scratch and then duplicate as needed. For this simple case let PETSc decide how many elements of the vector are stored on each processor. The second argument to VecSetSizes() below causes PETSc to decide. */ ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr); ierr = VecSetSizes(x,PETSC_DECIDE,n);CHKERRQ(ierr); ierr = VecSetFromOptions(x);CHKERRQ(ierr); ierr = VecDuplicate(x,&b);CHKERRQ(ierr); ierr = VecDuplicate(x,&u);CHKERRQ(ierr); /* Identify the starting and ending mesh points on each processor for the interior part of the mesh. We let PETSc decide above. */ ierr = VecGetOwnershipRange(x,&rstart,&rend);CHKERRQ(ierr); ierr = VecGetLocalSize(x,&nlocal);CHKERRQ(ierr); /* Create matrix. When using MatCreate(), the matrix format can be specified at runtime. Performance tuning note: For problems of substantial size, preallocation of matrix memory is crucial for attaining good performance. See the matrix chapter of the users manual for details. We pass in nlocal as the "local" size of the matrix to force it to have the same parallel layout as the vector created above. */ ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); ierr = MatSetSizes(A,nlocal,nlocal,n,n);CHKERRQ(ierr); ierr = MatSetFromOptions(A);CHKERRQ(ierr); ierr = MatSetUp(A);CHKERRQ(ierr); /* Assemble matrix. The linear system is distributed across the processors by chunks of contiguous rows, which correspond to contiguous sections of the mesh on which the problem is discretized. For matrix assembly, each processor contributes entries for the part that it owns locally. */ if (!rstart) { rstart = 1; 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);CHKERRQ(ierr); } if (rend == n) { rend = n-1; i = n-1; col[0] = n-2; col[1] = n-1; value[0] = -1.0; value[1] = 2.0; ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr); } /* Set entries corresponding to the mesh interior */ value[0] = -1.0; value[1] = 2.0; value[2] = -1.0; for (i=rstart; i -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 the linear system - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Solve linear system */ ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); /* View solver info; we could instead use the option -ksp_view to print this info to the screen at the conclusion of KSPSolve(). */ ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Check solution and clean up - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Check the error */ ierr = VecAXPY(x,-1.0,u);CHKERRQ(ierr); ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr); ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); if (norm > tol) { ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);CHKERRQ(ierr); } /* Free work space. All PETSc objects should be destroyed when they are no longer needed. */ ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&u);CHKERRQ(ierr); ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr); ierr = KSPDestroy(&ksp);CHKERRQ(ierr); /* Always call PetscFinalize() before exiting a program. This routine - finalizes the PETSc libraries as well as MPI - provides summary and diagnostic information if certain runtime options are chosen (e.g., -log_view). */ ierr = PetscFinalize(); return ierr; }