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