Actual source code: ex35.c
petsc-3.4.2 2013-07-02
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 <petscdmda.h>
12: #include <petscksp.h>
14: extern PetscErrorCode ComputeMatrix(DM,Mat);
15: extern PetscErrorCode ComputeRHS(DM,Vec);
19: int main(int argc,char **argv)
20: {
22: KSP ksp;
23: PC pc;
24: Vec x,b;
25: DM da;
26: Mat A;
27: PetscInt dof=1;
28: PetscBool flg;
29: PetscScalar zero=0.0;
31: PetscInitialize(&argc,&argv,(char*)0,help);
32: PetscOptionsGetInt(NULL,"-dof",&dof,NULL);
34: DMDACreate(PETSC_COMM_WORLD,&da);
35: DMDASetDim(da,3);
36: DMDASetBoundaryType(da,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE);
37: DMDASetStencilType(da,DMDA_STENCIL_STAR);
38: DMDASetSizes(da,3,3,3);
39: DMDASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
40: DMDASetDof(da,dof);
41: DMDASetStencilWidth(da,1);
42: DMDASetOwnershipRanges(da,NULL,NULL,NULL);
43: DMSetFromOptions(da);
44: DMSetUp(da);
46: DMCreateGlobalVector(da,&x);
47: DMCreateGlobalVector(da,&b);
48: DMCreateMatrix(da,MATAIJ,&A);
49: VecSet(b,zero);
51: /* Test sbaij matrix */
52: flg = PETSC_FALSE;
53: PetscOptionsGetBool(NULL,"-test_sbaij",&flg,NULL);
54: if (flg) {
55: Mat sA;
56: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
57: MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);
58: MatDestroy(&A);
59: A = sA;
60: }
62: KSPCreate(PETSC_COMM_WORLD,&ksp);
63: KSPSetFromOptions(ksp);
64: KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);
65: KSPGetPC(ksp,&pc);
66: PCSetDM(pc,(DM)da);
68: KSPSolve(ksp,b,x);
70: /* check final residual */
71: flg = PETSC_FALSE;
72: PetscOptionsGetBool(NULL, "-check_final_residual", &flg,NULL);
73: if (flg) {
74: Vec b1;
75: PetscReal norm;
76: KSPGetSolution(ksp,&x);
77: VecDuplicate(b,&b1);
78: MatMult(A,x,b1);
79: VecAXPY(b1,-1.0,b);
80: VecNorm(b1,NORM_2,&norm);
81: PetscPrintf(PETSC_COMM_WORLD,"Final residual %g\n",norm);
82: VecDestroy(&b1);
83: }
85: KSPDestroy(&ksp);
86: VecDestroy(&x);
87: VecDestroy(&b);
88: MatDestroy(&A);
89: DMDestroy(&da);
90: PetscFinalize();
91: return 0;
92: }