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
main(int argc,char ** argv)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