1 #include <petsc.h>
2 #include <petscvec.h>
3 #include <petscmat.h>
4 #include <petscksp.h>
5 #include "del2mat.h"
6
7 #define DEL2MAT_MULT ((void(*)(void))Del2Mat_mult)
8 #define DEL2MAT_DIAG ((void(*)(void))Del2Mat_diag)
9
main(int argc,char ** argv)10 int main(int argc,char **argv)
11 {
12 PetscInt n;
13 PetscScalar h;
14 Del2Mat shell;
15 Mat A;
16 Vec x,b;
17 KSP ksp;
18 PC pc;
19 PetscMPIInt size;
20 /* PETSc initialization */
21 PetscInitialize(&argc, &argv, NULL, NULL);
22 MPI_Comm_size(PETSC_COMM_WORLD,&size);
23 if (size != 1) {
24 PetscPrintf(PETSC_COMM_WORLD, "This a sequential example\n");
25 PetscFinalize();
26 return 1;
27 }
28 /* number of nodes in each direction
29 * excluding those at the boundary */
30 n = 32;
31 h = 1.0/(n+1); /* grid spacing */
32 /* setup linear system (shell) matrix */
33 MatCreate(PETSC_COMM_SELF, &A);
34 MatSetSizes(A, n*n*n, n*n*n, n*n*n, n*n*n);
35 MatSetType(A, MATSHELL);
36 shell.N = n;
37 PetscMalloc((n+2)*(n+2)*(n+2)*sizeof(PetscScalar),&shell.F);
38 PetscMemzero(shell.F, (n+2)*(n+2)*(n+2)*sizeof(PetscScalar));
39 MatShellSetContext(A, &shell);
40 MatShellSetOperation(A, MATOP_MULT, DEL2MAT_MULT);
41 MatShellSetOperation(A, MATOP_MULT_TRANSPOSE, DEL2MAT_MULT);
42 MatShellSetOperation(A, MATOP_GET_DIAGONAL, DEL2MAT_DIAG);
43 MatSetUp(A);
44 /* setup linear system vectors */
45 MatCreateVecs(A, &x, &b);
46 VecSet(x, 0);
47 VecSet(b, 1);
48 /* setup Krylov linear solver */
49 KSPCreate(PETSC_COMM_SELF, &ksp);
50 KSPGetPC(ksp, &pc);
51 KSPSetType(ksp, KSPCG); /* use conjugate gradients */
52 PCSetType(pc, PCNONE); /* with no preconditioning */
53 KSPSetFromOptions(ksp);
54 /* iteratively solve linear system of equations A*x=b */
55 KSPSetOperators(ksp,A,A);
56 KSPSolve(ksp, b, x);
57 /* scale solution vector to account for grid spacing */
58 VecScale(x, h*h);
59 /* free memory and destroy objects */
60 PetscFree(shell.F);
61 VecDestroy(&x);
62 VecDestroy(&b);
63 MatDestroy(&A);
64 KSPDestroy(&ksp);
65 /* finalize PETSc */
66 PetscFinalize();
67 return 0;
68 }
69