static char help[] = "Solves a linear system in parallel with KSP.\n\n"; /*T Concepts: KSP^solving a Helmholtz equation Concepts: complex numbers; Concepts: Helmholtz equation Processors: n T*/ /* Description: Solves a complex linear system in parallel with KSP. The model problem: Solve Helmholtz equation on the unit square: (0,1) x (0,1) -delta u - sigma1*u + i*sigma2*u = f, where delta = Laplace operator Dirichlet b.c.'s on all sides Use the 2-D, five-point finite difference stencil. Compiling the code: This code uses the complex numbers version of PETSc, so configure must be run to enable this */ /* 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 */ PetscInt dim,i,j,Ii,J,Istart,Iend,n = 6,its,use_random; PetscErrorCode ierr; PetscScalar v,none = -1.0,sigma2,pfive = 0.5,*xa; PetscRandom rctx; PetscReal h2,sigma1 = 100.0; PetscBool flg = PETSC_FALSE; ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; ierr = PetscOptionsGetReal(NULL,NULL,"-sigma1",&sigma1,NULL);CHKERRQ(ierr); ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); dim = n*n; /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Compute the matrix and right-hand-side vector that define the linear system, Ax = b. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* 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,dim,dim);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 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. */ ierr = PetscOptionsGetBool(NULL,NULL,"-norandom",&flg,NULL);CHKERRQ(ierr); if (flg) use_random = 0; else use_random = 1; if (use_random) { ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr); ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); ierr = PetscRandomSetInterval(rctx,0.0,PETSC_i);CHKERRQ(ierr); } else { sigma2 = 10.0*PETSC_i; } h2 = 1.0/((n+1)*(n+1)); 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 -pc_type -ksp_monitor -ksp_rtol */ ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Solve the linear system - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Check solution and clean up - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Print the first 3 entries of x; this demonstrates extraction of the real and imaginary components of the complex vector, x. */ flg = PETSC_FALSE; ierr = PetscOptionsGetBool(NULL,NULL,"-print_x3",&flg,NULL);CHKERRQ(ierr); if (flg) { ierr = VecGetArray(x,&xa);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"The first three entries of x are:\n");CHKERRQ(ierr); for (i=0; i<3; i++) { ierr = PetscPrintf(PETSC_COMM_WORLD,"x[%D] = %g + %g i\n",i,(double)PetscRealPart(xa[i]),(double)PetscImaginaryPart(xa[i]));CHKERRQ(ierr); } ierr = VecRestoreArray(x,&xa);CHKERRQ(ierr); } /* Check the error */ ierr = VecAXPY(x,none,u);CHKERRQ(ierr); ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr); ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); if (norm < 1.e-12) { ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error < 1.e-12 iterations %D\n",its);CHKERRQ(ierr); } else { 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 = KSPDestroy(&ksp);CHKERRQ(ierr); if (use_random) {ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);} ierr = VecDestroy(&u);CHKERRQ(ierr); ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr); ierr = PetscFinalize(); return ierr; }