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