Actual source code: ex5.c

  2: static char help[] = "Solves two linear systems in parallel with KSP.  The code\n\
  3: illustrates repeated solution of linear systems with the same preconditioner\n\
  4: method but different matrices (having the same nonzero structure).  The code\n\
  5: also uses multiple profiling stages.  Input arguments are\n\
  6:   -m <size> : problem size\n\
  7:   -mat_nonsym : use nonsymmetric matrix (default is symmetric)\n\n";

  9: /*T
 10:    Concepts: KSP^repeatedly solving linear systems;
 11:    Concepts: PetscLog^profiling multiple stages of code;
 12:    Processors: n
 13: T*/

 15: /* 
 16:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
 17:   automatically includes:
 18:      petsc.h       - base PETSc routines   petscvec.h - vectors
 19:      petscsys.h    - system routines       petscmat.h - matrices
 20:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 21:      petscviewer.h - viewers               petscpc.h  - preconditioners
 22: */
 23:  #include petscksp.h

 27: int main(int argc,char **args)
 28: {
 29:   KSP            ksp;             /* linear solver context */
 30:   Mat            C;                /* matrix */
 31:   Vec            x,u,b;          /* approx solution, RHS, exact solution */
 32:   PetscReal      norm;             /* norm of solution error */
 33:   PetscScalar    v,none = -1.0;
 34:   PetscInt       Ii,J,ldim,low,high,iglobal,Istart,Iend;
 36:   PetscInt       i,j,m = 3,n = 2,its;
 37:   PetscMPIInt    size,rank;
 38:   PetscTruth     mat_nonsymmetric;
 39: #if defined (PETSC_USE_LOG)
 40:   int            stages[2];
 41: #endif

 43:   PetscInitialize(&argc,&args,(char *)0,help);
 44:   PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
 45:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 46:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 47:   n = 2*size;

 49:   /*
 50:      Set flag if we are doing a nonsymmetric problem; the default is symmetric.
 51:   */
 52:   PetscOptionsHasName(PETSC_NULL,"-mat_nonsym",&mat_nonsymmetric);

 54:   /*
 55:      Register two stages for separate profiling of the two linear solves.
 56:      Use the runtime option -log_summary for a printout of performance
 57:      statistics at the program's conlusion.
 58:   */
 59:   PetscLogStageRegister(&stages[0],"Original Solve");
 60:   PetscLogStageRegister(&stages[1],"Second Solve");

 62:   /* -------------- Stage 0: Solve Original System ---------------------- */
 63:   /* 
 64:      Indicate to PETSc profiling that we're beginning the first stage
 65:   */
 66:   PetscLogStagePush(stages[0]);

 68:   /* 
 69:      Create parallel matrix, specifying only its global dimensions.
 70:      When using MatCreate(), the matrix format can be specified at
 71:      runtime. Also, the parallel partitioning of the matrix is
 72:      determined by PETSc at runtime.
 73:   */
 74:   MatCreate(PETSC_COMM_WORLD,&C);
 75:   MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
 76:   MatSetFromOptions(C);

 78:   /* 
 79:      Currently, all PETSc parallel matrix formats are partitioned by
 80:      contiguous chunks of rows across the processors.  Determine which
 81:      rows of the matrix are locally owned. 
 82:   */
 83:   MatGetOwnershipRange(C,&Istart,&Iend);

 85:   /* 
 86:      Set matrix entries matrix in parallel.
 87:       - Each processor needs to insert only elements that it owns
 88:         locally (but any non-local elements will be sent to the
 89:         appropriate processor during matrix assembly). 
 90:       - Always specify global row and columns of matrix entries.
 91:   */
 92:   for (Ii=Istart; Ii<Iend; Ii++) {
 93:     v = -1.0; i = Ii/n; j = Ii - i*n;
 94:     if (i>0)   {J = Ii - n; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
 95:     if (i<m-1) {J = Ii + n; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
 96:     if (j>0)   {J = Ii - 1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
 97:     if (j<n-1) {J = Ii + 1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
 98:     v = 4.0; MatSetValues(C,1,&Ii,1,&Ii,&v,ADD_VALUES);
 99:   }

101:   /*
102:      Make the matrix nonsymmetric if desired
103:   */
104:   if (mat_nonsymmetric) {
105:     for (Ii=Istart; Ii<Iend; Ii++) {
106:       v = -1.5; i = Ii/n;
107:       if (i>1)   {J = Ii-n-1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
108:     }
109:   } else {
110:     MatSetOption(C,MAT_SYMMETRIC);
111:     MatSetOption(C,MAT_SYMMETRY_ETERNAL);
112:   }

114:   /* 
115:      Assemble matrix, using the 2-step process:
116:        MatAssemblyBegin(), MatAssemblyEnd()
117:      Computations can be done while messages are in transition
118:      by placing code between these two statements.
119:   */
120:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
121:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

123:   /* 
124:      Create parallel vectors.
125:       - When using VecSetSizes(), we specify only the vector's global
126:         dimension; the parallel partitioning is determined at runtime. 
127:       - Note: We form 1 vector from scratch and then duplicate as needed.
128:   */
129:   VecCreate(PETSC_COMM_WORLD,&u);
130:   VecSetSizes(u,PETSC_DECIDE,m*n);
131:   VecSetFromOptions(u);
132:   VecDuplicate(u,&b);
133:   VecDuplicate(b,&x);

135:   /* 
136:      Currently, all parallel PETSc vectors are partitioned by
137:      contiguous chunks across the processors.  Determine which
138:      range of entries are locally owned.
139:   */
140:   VecGetOwnershipRange(x,&low,&high);

142:   /*
143:     Set elements within the exact solution vector in parallel.
144:      - Each processor needs to insert only elements that it owns
145:        locally (but any non-local entries will be sent to the
146:        appropriate processor during vector assembly).
147:      - Always specify global locations of vector entries.
148:   */
149:   VecGetLocalSize(x,&ldim);
150:   for (i=0; i<ldim; i++) {
151:     iglobal = i + low;
152:     v = (PetscScalar)(i + 100*rank);
153:     VecSetValues(u,1,&iglobal,&v,INSERT_VALUES);
154:   }

156:   /* 
157:      Assemble vector, using the 2-step process:
158:        VecAssemblyBegin(), VecAssemblyEnd()
159:      Computations can be done while messages are in transition,
160:      by placing code between these two statements.
161:   */
162:   VecAssemblyBegin(u);
163:   VecAssemblyEnd(u);

165:   /* 
166:      Compute right-hand-side vector
167:   */
168:   MatMult(C,u,b);
169: 
170:   /* 
171:     Create linear solver context
172:   */
173:   KSPCreate(PETSC_COMM_WORLD,&ksp);

175:   /* 
176:      Set operators. Here the matrix that defines the linear system
177:      also serves as the preconditioning matrix.
178:   */
179:   KSPSetOperators(ksp,C,C,DIFFERENT_NONZERO_PATTERN);

181:   /* 
182:      Set runtime options (e.g., -ksp_type <type> -pc_type <type>)
183:   */

185:   KSPSetFromOptions(ksp);

187:   /* 
188:      Solve linear system.  Here we explicitly call KSPSetUp() for more
189:      detailed performance monitoring of certain preconditioners, such
190:      as ICC and ILU.  This call is optional, as KSPSetUp() will
191:      automatically be called within KSPSolve() if it hasn't been
192:      called already.
193:   */
194:   KSPSetUp(ksp);
195:   KSPSolve(ksp,b,x);
196: 
197:   /* 
198:      Check the error
199:   */
200:   VecAXPY(x,none,u);
201:   VecNorm(x,NORM_2,&norm);
202:   KSPGetIterationNumber(ksp,&its);
203:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A, Iterations %D\n",norm,its);

205:   /* -------------- Stage 1: Solve Second System ---------------------- */
206:   /* 
207:      Solve another linear system with the same method.  We reuse the KSP
208:      context, matrix and vector data structures, and hence save the
209:      overhead of creating new ones.

211:      Indicate to PETSc profiling that we're concluding the first
212:      stage with PetscLogStagePop(), and beginning the second stage with
213:      PetscLogStagePush().
214:   */
215:   PetscLogStagePop();
216:   PetscLogStagePush(stages[1]);

218:   /* 
219:      Initialize all matrix entries to zero.  MatZeroEntries() retains the
220:      nonzero structure of the matrix for sparse formats.
221:   */
222:   MatZeroEntries(C);

224:   /* 
225:      Assemble matrix again.  Note that we retain the same matrix data
226:      structure and the same nonzero pattern; we just change the values
227:      of the matrix entries.
228:   */
229:   for (i=0; i<m; i++) {
230:     for (j=2*rank; j<2*rank+2; j++) {
231:       v = -1.0;  Ii = j + n*i;
232:       if (i>0)   {J = Ii - n; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
233:       if (i<m-1) {J = Ii + n; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
234:       if (j>0)   {J = Ii - 1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
235:       if (j<n-1) {J = Ii + 1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
236:       v = 6.0; MatSetValues(C,1,&Ii,1,&Ii,&v,ADD_VALUES);
237:     }
238:   }
239:   if (mat_nonsymmetric) {
240:     for (Ii=Istart; Ii<Iend; Ii++) {
241:       v = -1.5; i = Ii/n;
242:       if (i>1)   {J = Ii-n-1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
243:     }
244:   }
245:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
246:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

248:   /* 
249:      Compute another right-hand-side vector
250:   */
251:   MatMult(C,u,b);

253:   /* 
254:      Set operators. Here the matrix that defines the linear system
255:      also serves as the preconditioning matrix.
256:       - The flag SAME_NONZERO_PATTERN indicates that the
257:         preconditioning matrix has identical nonzero structure
258:         as during the last linear solve (although the values of
259:         the entries have changed). Thus, we can save some
260:         work in setting up the preconditioner (e.g., no need to
261:         redo symbolic factorization for ILU/ICC preconditioners).
262:       - If the nonzero structure of the matrix is different during
263:         the second linear solve, then the flag DIFFERENT_NONZERO_PATTERN
264:         must be used instead.  If you are unsure whether the
265:         matrix structure has changed or not, use the flag
266:         DIFFERENT_NONZERO_PATTERN.
267:       - Caution:  If you specify SAME_NONZERO_PATTERN, PETSc
268:         believes your assertion and does not check the structure
269:         of the matrix.  If you erroneously claim that the structure
270:         is the same when it actually is not, the new preconditioner
271:         will not function correctly.  Thus, use this optimization
272:         feature with caution!
273:   */
274:   KSPSetOperators(ksp,C,C,SAME_NONZERO_PATTERN);

276:   /* 
277:      Solve linear system
278:   */
279:   KSPSetUp(ksp);
280:   KSPSolve(ksp,b,x);

282:   /* 
283:      Check the error
284:   */
285:   VecAXPY(x,none,u);
286:   VecNorm(x,NORM_2,&norm);
287:   KSPGetIterationNumber(ksp,&its);
288:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A, Iterations %D\n",norm,its);

290:   /* 
291:      Free work space.  All PETSc objects should be destroyed when they
292:      are no longer needed.
293:   */
294:   KSPDestroy(ksp);
295:   VecDestroy(u);
296:   VecDestroy(x);
297:   VecDestroy(b);
298:   MatDestroy(C);

300:   /*
301:      Indicate to PETSc profiling that we're concluding the second stage 
302:   */
303:   PetscLogStagePop();

305:   PetscFinalize();
306:   return 0;
307: }