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