Actual source code: ex35.c
petsc-3.5.4 2015-05-23
2: /*
3: Used for testing AIJ matrix with all zeros.
4: */
6: static char help[] = "Used for Solving a linear system where the matrix has all zeros.\n\n";
7: /*
8: Example: mpiexec -n <np> ./ex35 -dof 2 -mat_view -check_final_residual
9: */
11: #include <petscdm.h>
12: #include <petscdmda.h>
13: #include <petscksp.h>
17: int main(int argc,char **argv)
18: {
20: KSP ksp;
21: PC pc;
22: Vec x,b;
23: DM da;
24: Mat A;
25: PetscInt dof=1;
26: PetscBool flg;
27: PetscScalar zero=0.0;
29: PetscInitialize(&argc,&argv,(char*)0,help);
30: PetscOptionsGetInt(NULL,"-dof",&dof,NULL);
32: DMDACreate(PETSC_COMM_WORLD,&da);
33: DMDASetDim(da,3);
34: DMDASetBoundaryType(da,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE);
35: DMDASetStencilType(da,DMDA_STENCIL_STAR);
36: DMDASetSizes(da,3,3,3);
37: DMDASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
38: DMDASetDof(da,dof);
39: DMDASetStencilWidth(da,1);
40: DMDASetOwnershipRanges(da,NULL,NULL,NULL);
41: DMSetFromOptions(da);
42: DMSetUp(da);
44: DMCreateGlobalVector(da,&x);
45: DMCreateGlobalVector(da,&b);
46: DMSetMatType(da,MATAIJ);
47: DMCreateMatrix(da,&A);
48: VecSet(b,zero);
50: /* Test sbaij matrix */
51: flg = PETSC_FALSE;
52: PetscOptionsGetBool(NULL,"-test_sbaij",&flg,NULL);
53: if (flg) {
54: Mat sA;
55: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
56: MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);
57: MatDestroy(&A);
58: A = sA;
59: }
61: KSPCreate(PETSC_COMM_WORLD,&ksp);
62: KSPSetFromOptions(ksp);
63: KSPSetOperators(ksp,A,A);
64: KSPGetPC(ksp,&pc);
65: PCSetDM(pc,(DM)da);
67: KSPSolve(ksp,b,x);
69: /* check final residual */
70: flg = PETSC_FALSE;
71: PetscOptionsGetBool(NULL, "-check_final_residual", &flg,NULL);
72: if (flg) {
73: Vec b1;
74: PetscReal norm;
75: KSPGetSolution(ksp,&x);
76: VecDuplicate(b,&b1);
77: MatMult(A,x,b1);
78: VecAXPY(b1,-1.0,b);
79: VecNorm(b1,NORM_2,&norm);
80: PetscPrintf(PETSC_COMM_WORLD,"Final residual %g\n",norm);
81: VecDestroy(&b1);
82: }
84: KSPDestroy(&ksp);
85: VecDestroy(&x);
86: VecDestroy(&b);
87: MatDestroy(&A);
88: DMDestroy(&da);
89: PetscFinalize();
90: return 0;
91: }