Actual source code: ex28.c
2: static char help[] = "Test procedural KSPSetFromOptions() or at runtime; Test PCREDUNDANT.\n\n";
4: #include <petscksp.h>
6: int main(int argc,char **args)
7: {
8: Vec x, b, u; /* approx solution, RHS, exact solution */
9: Mat A; /* linear system matrix */
10: KSP ksp; /* linear solver context */
11: PC pc; /* preconditioner context */
12: PetscReal norm; /* norm of solution error */
13: PetscInt i,n = 10,col[3],its,rstart,rend,nlocal;
14: PetscScalar one = 1.0,value[3];
15: PetscBool TEST_PROCEDURAL=PETSC_FALSE;
17: PetscInitialize(&argc,&args,(char*)0,help);
18: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
19: PetscOptionsGetBool(NULL,NULL,"-procedural",&TEST_PROCEDURAL,NULL);
21: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
22: Compute the matrix and right-hand-side vector that define
23: the linear system, Ax = b.
24: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
26: /*
27: Create vectors. Note that we form 1 vector from scratch and
28: then duplicate as needed. For this simple case let PETSc decide how
29: many elements of the vector are stored on each processor. The second
30: argument to VecSetSizes() below causes PETSc to decide.
31: */
32: VecCreate(PETSC_COMM_WORLD,&x);
33: VecSetSizes(x,PETSC_DECIDE,n);
34: VecSetFromOptions(x);
35: VecDuplicate(x,&b);
36: VecDuplicate(x,&u);
38: /* Identify the starting and ending mesh points on each
39: processor for the interior part of the mesh. We let PETSc decide
40: above. */
42: VecGetOwnershipRange(x,&rstart,&rend);
43: VecGetLocalSize(x,&nlocal);
45: /* Create a tridiagonal matrix. See ../tutorials/ex23.c */
46: MatCreate(PETSC_COMM_WORLD,&A);
47: MatSetSizes(A,nlocal,nlocal,n,n);
48: MatSetFromOptions(A);
49: MatSetUp(A);
50: /* Assemble matrix */
51: if (!rstart) {
52: rstart = 1;
53: i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0;
54: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
55: }
56: if (rend == n) {
57: rend = n-1;
58: i = n-1; col[0] = n-2; col[1] = n-1; value[0] = -1.0; value[1] = 2.0;
59: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
60: }
62: /* Set entries corresponding to the mesh interior */
63: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
64: for (i=rstart; i<rend; i++) {
65: col[0] = i-1; col[1] = i; col[2] = i+1;
66: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
67: }
68: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
69: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
71: /* Set exact solution; then compute right-hand-side vector. */
72: VecSet(u,one);
73: MatMult(A,u,b);
75: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76: Create the linear solver and set various options
77: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
78: KSPCreate(PETSC_COMM_WORLD,&ksp);
79: KSPSetOperators(ksp,A,A);
81: /*
82: Set linear solver defaults for this problem (optional).
83: - By extracting the KSP and PC contexts from the KSP context,
84: we can then directly call any KSP and PC routines to set
85: various options.
86: - The following statements are optional; all of these
87: parameters could alternatively be specified at runtime via
88: KSPSetFromOptions();
89: */
90: if (TEST_PROCEDURAL) {
91: /* Example of runtime options: '-pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type bjacobi' */
92: PetscMPIInt size,rank,subsize;
93: Mat A_redundant;
94: KSP innerksp;
95: PC innerpc;
96: MPI_Comm subcomm;
98: KSPGetPC(ksp,&pc);
99: PCSetType(pc,PCREDUNDANT);
100: MPI_Comm_size(PETSC_COMM_WORLD,&size);
101: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
103: PCRedundantSetNumber(pc,size-2);
104: KSPSetFromOptions(ksp);
106: /* Get subcommunicator and redundant matrix */
107: KSPSetUp(ksp);
108: PCRedundantGetKSP(pc,&innerksp);
109: KSPGetPC(innerksp,&innerpc);
110: PCGetOperators(innerpc,NULL,&A_redundant);
111: PetscObjectGetComm((PetscObject)A_redundant,&subcomm);
112: MPI_Comm_size(subcomm,&subsize);
113: if (subsize==1 && rank == 0) {
114: PetscPrintf(PETSC_COMM_SELF,"A_redundant:\n");
115: MatView(A_redundant,PETSC_VIEWER_STDOUT_SELF);
116: }
117: } else {
118: KSPSetFromOptions(ksp);
119: }
121: /* Solve linear system */
122: KSPSolve(ksp,b,x);
124: /* Check the error */
125: VecAXPY(x,-1.0,u);
126: VecNorm(x,NORM_2,&norm);
127: KSPGetIterationNumber(ksp,&its);
128: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
130: /* Free work space. */
131: VecDestroy(&x)); PetscCall(VecDestroy(&u);
132: VecDestroy(&b)); PetscCall(MatDestroy(&A);
133: KSPDestroy(&ksp);
134: PetscFinalize();
135: return 0;
136: }
138: /*TEST
140: test:
141: nsize: 3
142: output_file: output/ex28.out
144: test:
145: suffix: 2
146: args: -procedural -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type bjacobi
147: nsize: 3
149: test:
150: suffix: 3
151: args: -procedural -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type bjacobi
152: nsize: 5
154: TEST*/