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