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 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