Actual source code: bench_kspsolve.c

  1: /*
  2: Poisson in 3D. Modeled by the PDE:

  4:   - delta u(x,y,z) = f(x,y,z)

  6: With exact solution:

  8:    u(x,y,z) = 1.0

 10: Example usage:

 12:   Run on GPU (requires respective backends installed):
 13:     ./bench_kspsolve -mat_type aijcusparse
 14:     ./bench_kspsolve -mat_type aijhipsparse
 15:     ./bench_kspsolve -mat_type aijkokkos

 17:   Test only MatMult:
 18:     ./bench_kspsolve -matmult

 20:   Test MatMult over 1000 iterations:
 21:     ./bench_kspsolve -matmult -its 1000

 23:   Change size of problem (e.g., use a 128x128x128 grid):
 24:     ./bench_kspsolve -n 128
 25: */
 26: static char help[] = "Solves 3D Laplacian with 27-point finite difference stencil.\n";

 28: #include <petscksp.h>

 30: typedef struct {
 31:   PetscMPIInt rank, size;
 32:   PetscInt    n;        /* global size in each dimension */
 33:   PetscInt    dim;      /* global size */
 34:   PetscInt    nnz;      /* local nonzeros */
 35:   PetscBool   matmult;  /* Do MatMult() only */
 36:   PetscBool   debug;    /* Debug flag for PreallocateCOO() */
 37:   PetscBool   splitksp; /* Split KSPSolve and PCSetUp */
 38:   PetscInt    its;      /* No of matmult_iterations */
 39:   PetscInt    Istart, Iend;
 40: } AppCtx;

 42: static PetscErrorCode PreallocateCOO(Mat A, void *ctx)
 43: {
 44:   AppCtx  *user = (AppCtx *)ctx;
 45:   PetscInt n = user->n, n2 = n * n, n1 = n - 1;
 46:   PetscInt xstart, ystart, zstart, xend, yend, zend, nm2 = n - 2, idx;
 47:   PetscInt nnz[] = {8, 12, 18, 27}; /* nnz for corner, edge, face, and center. */

 49:   PetscFunctionBeginUser;
 50:   xstart = user->Istart % n;
 51:   ystart = (user->Istart / n) % n;
 52:   zstart = user->Istart / n2;
 53:   xend   = (user->Iend - 1) % n;
 54:   yend   = ((user->Iend - 1) / n) % n;
 55:   zend   = (user->Iend - 1) / n2;

 57:   /* bottom xy-plane */
 58:   user->nnz = 0;
 59:   idx       = !zstart ? 0 : 1;
 60:   if (zstart == zend && (!zstart || zstart == n1)) idx = 0;
 61:   if (zstart == zend && (zstart && zstart < n1)) idx = 1;
 62:   if (!xstart && !ystart) // bottom left
 63:     user->nnz += 4 * nnz[idx] + 4 * nm2 * nnz[idx + 1] + nm2 * nm2 * nnz[idx + 2];
 64:   else if (xstart && xstart < n1 && !ystart) // bottom
 65:     user->nnz += 3 * nnz[idx] + (3 * nm2 + n1 - xstart) * nnz[idx + 1] + nm2 * nm2 * nnz[idx + 2];
 66:   else if (xstart == n1 && !ystart) // bottom right
 67:     user->nnz += 3 * nnz[idx] + (3 * nm2) * nnz[idx + 1] + nm2 * nm2 * nnz[idx + 2];
 68:   else if (!xstart && ystart && ystart < n1) // left
 69:     user->nnz += 2 * nnz[idx] + (nm2 + 2 * (n1 - ystart)) * nnz[idx + 1] + nm2 * (n1 - ystart) * nnz[idx + 2];
 70:   else if (xstart && xstart < n1 && ystart && ystart < n1) // center
 71:     user->nnz += 2 * nnz[idx] + (nm2 + (n1 - ystart) + (nm2 - ystart)) * nnz[idx + 1] + (nm2 * (nm2 - ystart) + (n1 - xstart)) * nnz[idx + 2];
 72:   else if (xstart == n1 && ystart && ystart < n1) // right
 73:     user->nnz += 2 * nnz[idx] + (nm2 + n1 - ystart + nm2 - ystart) * nnz[idx + 1] + nm2 * (nm2 - ystart) * nnz[idx + 2];
 74:   else if (ystart == n1 && !xstart) // top left
 75:     user->nnz += 2 * nnz[idx] + nm2 * nnz[idx + 1];
 76:   else if (ystart == n1 && xstart && xstart < n1) // top
 77:     user->nnz += nnz[idx] + (n1 - xstart) * nnz[idx + 1];
 78:   else if (xstart == n1 && ystart == n1) // top right
 79:     user->nnz += nnz[idx];

 81:   /* top and middle plane the same */
 82:   if (zend == zstart) user->nnz = user->nnz - 4 * nnz[idx] - 4 * nm2 * nnz[idx + 1] - nm2 * nm2 * nnz[idx + 2];

 84:   /* middle xy-planes */
 85:   if (zend - zstart > 1) user->nnz = user->nnz + (zend - zstart - 1) * (4 * nnz[1] + 4 * nnz[2] * nm2 + nnz[3] * nm2 * nm2);

 87:   /* top xy-plane */
 88:   idx = zend == n1 ? 0 : 1;
 89:   if (zstart == zend && (!zend || zend == n1)) idx = 0;
 90:   if (zstart == zend && (zend && zend < n1)) idx = 1;
 91:   if (!xend && !yend) // bottom left
 92:     user->nnz += nnz[idx];
 93:   else if (xend && xend < n1 && !yend) // bottom
 94:     user->nnz += nnz[idx] + xend * nnz[idx + 1];
 95:   else if (xend == n1 && !yend) // bottom right
 96:     user->nnz += 2 * nnz[idx] + nm2 * nnz[idx + 1];
 97:   else if (!xend && yend && yend < n1) // left
 98:     user->nnz += 2 * nnz[idx] + (nm2 + yend + yend - 1) * nnz[idx + 1] + nm2 * (yend - 1) * nnz[idx + 2];
 99:   else if (xend && xend < n1 && yend && yend < n1) // center
100:     user->nnz += 2 * nnz[idx] + (nm2 + yend + yend - 1) * nnz[idx + 1] + (nm2 * (yend - 1) + xend) * nnz[idx + 2];
101:   else if (xend == n1 && yend && yend < n1) // right
102:     user->nnz += 2 * nnz[idx] + (nm2 + 2 * yend) * nnz[idx + 1] + (nm2 * yend) * nnz[idx + 2];
103:   else if (!xend && yend == n1) // top left
104:     user->nnz += 3 * nnz[idx] + (3 * nm2) * nnz[idx + 1] + (nm2 * nm2) * nnz[idx + 2];
105:   else if (xend && xend < n1 && yend == n1) // top
106:     user->nnz += 3 * nnz[idx] + (3 * nm2 + xend) * nnz[idx + 1] + (nm2 * nm2) * nnz[idx + 2];
107:   else if (xend == n1 && yend == n1) // top right
108:     user->nnz += 4 * nnz[idx] + (4 * nm2) * nnz[idx + 1] + (nm2 * nm2) * nnz[idx + 2];
109:   if (user->debug)
110:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "rank %d: xyzstart = %" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ", xvzend = %" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ", nnz = %" PetscInt_FMT "\n", user->rank, xstart, ystart, zstart, xend, yend, zend,
111:                           user->nnz));
112:   PetscFunctionReturn(PETSC_SUCCESS);
113: }

115: PetscErrorCode FillCOO(Mat A, void *ctx)
116: {
117:   AppCtx      *user = (AppCtx *)ctx;
118:   PetscInt     Ii, x, y, z, n = user->n, n2 = n * n, n1 = n - 1;
119:   PetscInt     count = 0;
120:   PetscInt    *coo_i, *coo_j;
121:   PetscScalar *coo_v;
122:   PetscScalar  h     = 1.0 / (n - 1);
123:   PetscScalar  vcorn = -1.0 / 13 * h; //-1.0/12.0*h;
124:   PetscScalar  vedge = -3.0 / 26 * h; //-1.0/6.0*h;
125:   PetscScalar  vface = -3.0 / 13 * h; //-1.0e-9*h;
126:   PetscScalar  vcent = 44.0 / 13 * h; //8.0/3.0*h;

128:   PetscFunctionBeginUser;
129:   PetscCall(PetscCalloc3(user->nnz, &coo_i, user->nnz, &coo_j, user->nnz, &coo_v));
130:   /* Fill COO arrays */
131:   for (Ii = user->Istart; Ii < user->Iend; Ii++) {
132:     x = Ii % n;
133:     y = (Ii / n) % n;
134:     z = Ii / n2;
135:     if (x > 0 && y > 0 && z > 0) {
136:       coo_i[count] = Ii;
137:       coo_j[count] = Ii - 1 - n - n2;
138:       coo_v[count] = vcorn;
139:       count++;
140:     }
141:     if (x > 0 && y < n1 && z > 0) {
142:       coo_i[count] = Ii;
143:       coo_j[count] = Ii - 1 + n - n2;
144:       coo_v[count] = vcorn;
145:       count++;
146:     }
147:     if (x < n1 && y > 0 && z > 0) {
148:       coo_i[count] = Ii;
149:       coo_j[count] = Ii + 1 - n - n2;
150:       coo_v[count] = vcorn;
151:       count++;
152:     }
153:     if (x < n1 && y < n1 && z > 0) {
154:       coo_i[count] = Ii;
155:       coo_j[count] = Ii + 1 + n - n2;
156:       coo_v[count] = vcorn;
157:       count++;
158:     }
159:     if (x > 0 && y > 0 && z < n1) {
160:       coo_i[count] = Ii;
161:       coo_j[count] = Ii - 1 - n + n2;
162:       coo_v[count] = vcorn;
163:       count++;
164:     }
165:     if (x > 0 && y < n1 && z < n1) {
166:       coo_i[count] = Ii;
167:       coo_j[count] = Ii - 1 + n + n2;
168:       coo_v[count] = vcorn;
169:       count++;
170:     }
171:     if (x < n1 && y > 0 && z < n1) {
172:       coo_i[count] = Ii;
173:       coo_j[count] = Ii + 1 - n + n2;
174:       coo_v[count] = vcorn;
175:       count++;
176:     }
177:     if (x < n1 && y < n1 && z < n1) {
178:       coo_i[count] = Ii;
179:       coo_j[count] = Ii + 1 + n + n2;
180:       coo_v[count] = vcorn;
181:       count++;
182:     }
183:     /* 12 edges */
184:     if (x > 0 && y > 0) {
185:       coo_i[count] = Ii;
186:       coo_j[count] = Ii - 1 - n;
187:       coo_v[count] = vedge;
188:       count++;
189:     }
190:     if (x > 0 && y < n1) {
191:       coo_i[count] = Ii;
192:       coo_j[count] = Ii - 1 + n;
193:       coo_v[count] = vedge;
194:       count++;
195:     }
196:     if (x < n1 && y > 0) {
197:       coo_i[count] = Ii;
198:       coo_j[count] = Ii + 1 - n;
199:       coo_v[count] = vedge;
200:       count++;
201:     }
202:     if (x < n1 && y < n1) {
203:       coo_i[count] = Ii;
204:       coo_j[count] = Ii + 1 + n;
205:       coo_v[count] = vedge;
206:       count++;
207:     }
208:     if (x > 0 && z > 0) {
209:       coo_i[count] = Ii;
210:       coo_j[count] = Ii - 1 - n2;
211:       coo_v[count] = vedge;
212:       count++;
213:     }
214:     if (x > 0 && z < n1) {
215:       coo_i[count] = Ii;
216:       coo_j[count] = Ii - 1 + n2;
217:       coo_v[count] = vedge;
218:       count++;
219:     }
220:     if (x < n1 && z > 0) {
221:       coo_i[count] = Ii;
222:       coo_j[count] = Ii + 1 - n2;
223:       coo_v[count] = vedge;
224:       count++;
225:     }
226:     if (x < n1 && z < n1) {
227:       coo_i[count] = Ii;
228:       coo_j[count] = Ii + 1 + n2;
229:       coo_v[count] = vedge;
230:       count++;
231:     }
232:     if (y > 0 && z > 0) {
233:       coo_i[count] = Ii;
234:       coo_j[count] = Ii - n - n2;
235:       coo_v[count] = vedge;
236:       count++;
237:     }
238:     if (y > 0 && z < n1) {
239:       coo_i[count] = Ii;
240:       coo_j[count] = Ii - n + n2;
241:       coo_v[count] = vedge;
242:       count++;
243:     }
244:     if (y < n1 && z > 0) {
245:       coo_i[count] = Ii;
246:       coo_j[count] = Ii + n - n2;
247:       coo_v[count] = vedge;
248:       count++;
249:     }
250:     if (y < n1 && z < n1) {
251:       coo_i[count] = Ii;
252:       coo_j[count] = Ii + n + n2;
253:       coo_v[count] = vedge;
254:       count++;
255:     }
256:     /* 6 faces */
257:     if (x > 0) {
258:       coo_i[count] = Ii;
259:       coo_j[count] = Ii - 1;
260:       coo_v[count] = vface;
261:       count++;
262:     }
263:     if (x < n1) {
264:       coo_i[count] = Ii;
265:       coo_j[count] = Ii + 1;
266:       coo_v[count] = vface;
267:       count++;
268:     }
269:     if (y > 0) {
270:       coo_i[count] = Ii;
271:       coo_j[count] = Ii - n;
272:       coo_v[count] = vface;
273:       count++;
274:     }
275:     if (y < n1) {
276:       coo_i[count] = Ii;
277:       coo_j[count] = Ii + n;
278:       coo_v[count] = vface;
279:       count++;
280:     }
281:     if (z > 0) {
282:       coo_i[count] = Ii;
283:       coo_j[count] = Ii - n2;
284:       coo_v[count] = vface;
285:       count++;
286:     }
287:     if (z < n1) {
288:       coo_i[count] = Ii;
289:       coo_j[count] = Ii + n2;
290:       coo_v[count] = vface;
291:       count++;
292:     }
293:     /* Center */
294:     coo_i[count] = Ii;
295:     coo_j[count] = Ii;
296:     coo_v[count] = vcent;
297:     count++;
298:   }
299:   PetscCheck(count == user->nnz, PETSC_COMM_SELF, PETSC_ERR_USER_INPUT, "Expected %" PetscInt_FMT " nonzeros but got %" PetscInt_FMT " nonzeros in COO format", user->nnz, count);
300:   PetscCall(MatSetPreallocationCOO(A, user->nnz, coo_i, coo_j));
301:   PetscCall(MatSetValuesCOO(A, coo_v, INSERT_VALUES));
302:   PetscCall(PetscFree3(coo_i, coo_j, coo_v));
303:   PetscFunctionReturn(PETSC_SUCCESS);
304: }

306: int main(int argc, char **argv)
307: {
308:   AppCtx                     user;                               /* Application context */
309:   Vec                        x, b, u;                            /* approx solution, RHS, and exact solution */
310:   Mat                        A;                                  /* linear system matrix */
311:   KSP                        ksp;                                /* linear solver context */
312:   PC                         pc;                                 /* Preconditioner */
313:   PetscReal                  norm;                               /* Error norm */
314:   PetscInt                   nlocal     = PETSC_DECIDE, ksp_its; /* Number of KSP iterations */
315:   PetscInt                   global_nnz = 0;                     /* Total number of nonzeros */
316:   PetscLogDouble             time_start, time_mid1 = 0.0, time_mid2 = 0.0, time_end, time_avg, floprate;
317:   PetscBool                  printTiming = PETSC_TRUE; /* If run in CI, do not print timing result */
318:   PETSC_UNUSED PetscLogStage stage;

320:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
321:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &user.size));
322:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &user.rank));

324:   user.n        = 100;         /* Default grid points per dimension */
325:   user.matmult  = PETSC_FALSE; /* Test MatMult only */
326:   user.its      = 100;         /* Default no. of iterations for MatMult test */
327:   user.debug    = PETSC_FALSE; /* Debug PreallocateCOO() */
328:   user.splitksp = PETSC_FALSE; /* Split KSPSolve and PCSetUp */
329:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &user.n, NULL));
330:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-its", &user.its, NULL));
331:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-matmult", &user.matmult, NULL));
332:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-debug", &user.debug, NULL));
333:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-print_timing", &printTiming, NULL));
334:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-split_ksp", &user.splitksp, NULL));

336:   user.dim   = user.n * user.n * user.n;
337:   global_nnz = 64 + 27 * (user.n - 2) * (user.n - 2) * (user.n - 2) + 108 * (user.n - 2) * (user.n - 2) + 144 * (user.n - 2);
338:   PetscCheck(user.n >= 2, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Requires at least 2 grid points (-n 2), you specified -n %" PetscInt_FMT, user.n);
339:   PetscCheck(user.dim >= user.size, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "MPI size (%d) exceeds the grid size %" PetscInt_FMT " (-n %" PetscInt_FMT ")", user.size, user.dim, user.n);
340:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "===========================================\n"));
341:   if (user.matmult) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test: MatMult performance - Poisson\n"));
342:   else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test: KSP performance - Poisson\n"));
343:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\tInput matrix: 27-pt finite difference stencil\n"));
344:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\t-n %" PetscInt_FMT "\n", user.n));
345:   if (user.matmult) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\t-its %" PetscInt_FMT "\n", user.its));
346:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\tDoFs = %" PetscInt_FMT "\n", user.dim));
347:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\tNumber of nonzeros = %" PetscInt_FMT "\n", global_nnz));

349:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
350:   /*  Create the Vecs and Mat                                            */
351:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
352:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nStep1  - creating Vecs and Mat...\n"));
353:   PetscCall(PetscLogStageRegister("Step1  - Vecs and Mat", &stage));
354:   PetscCall(PetscLogStagePush(stage));
355:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
356:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, user.dim, user.dim));
357:   PetscCall(MatSetFromOptions(A));

359:   /* cannot use MatGetOwnershipRange() because the matrix has yet to be preallocated; that happens in MatSetPreallocateCOO() */
360:   PetscCall(PetscSplitOwnership(PetscObjectComm((PetscObject)A), &nlocal, &user.dim));
361:   PetscCallMPI(MPI_Scan(&nlocal, &user.Istart, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A)));
362:   user.Istart -= nlocal;
363:   user.Iend = user.Istart + nlocal;

365:   PetscCall(PreallocateCOO(A, &user)); /* Determine local number of nonzeros */
366:   PetscCall(FillCOO(A, &user));        /* Fill COO Matrix */
367: #if !defined(PETSC_HAVE_HIP)           /* Due to errors in MatSolve_SeqAIJHIPSPARSE_ICC0() */
368:   PetscCall(MatSetOption(A, MAT_SPD, PETSC_TRUE));
369:   PetscCall(MatSetOption(A, MAT_SPD_ETERNAL, PETSC_TRUE));
370: #endif
371:   PetscCall(MatCreateVecs(A, &u, &b));
372:   if (!user.matmult) PetscCall(VecDuplicate(b, &x));
373:   PetscCall(VecSet(u, 1.0));   /* Exact solution */
374:   PetscCall(MatMult(A, u, b)); /* Compute RHS based on exact solution */
375:   PetscCall(PetscLogStagePop());

377:   if (user.matmult) {
378:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
379:     /*  MatMult                                                            */
380:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
381:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Step2  - running MatMult() %" PetscInt_FMT " times...\n", user.its));
382:     PetscCall(PetscLogStageRegister("Step2  - MatMult", &stage));
383:     PetscCall(PetscLogStagePush(stage));
384:     PetscCall(PetscTime(&time_start));
385:     for (int i = 0; i < user.its; i++) PetscCall(MatMult(A, u, b));
386:     PetscCall(PetscTime(&time_end));
387:     PetscCall(PetscLogStagePop());
388:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
389:     /*  Calculate Performance metrics                                      */
390:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
391:     time_avg = (time_end - time_start) / ((PetscLogDouble)user.its);
392:     floprate = 2 * global_nnz / time_avg * 1e-9;
393:     if (printTiming) {
394:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n%-15s%-7.5f seconds\n", "Average time:", time_avg));
395:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%-15s%-9.3e Gflops/sec\n", "FOM:", floprate)); /* figure of merit */
396:     }
397:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "===========================================\n"));
398:   } else {
399:     if (!user.splitksp) {
400:       /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
401:       /*  Solve the linear system of equations                               */
402:       /*  Measure only time of PCSetUp() and KSPSolve()                      */
403:       /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
404:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Step2  - running KSPSolve()...\n"));
405:       PetscCall(PetscLogStageRegister("Step2  - KSPSolve", &stage));
406:       PetscCall(PetscLogStagePush(stage));
407:       PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
408:       PetscCall(KSPSetOperators(ksp, A, A));
409:       PetscCall(KSPSetFromOptions(ksp));
410:       PetscCall(PetscTime(&time_start));
411:       PetscCall(KSPSolve(ksp, b, x));
412:       PetscCall(PetscTime(&time_end));
413:       PetscCall(PetscLogStagePop());
414:     } else {
415:       /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
416:       /*  Solve the linear system of equations                               */
417:       /*  Measure only time of PCSetUp() and KSPSolve()                      */
418:       /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
419:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Step2a - running PCSetUp()...\n"));
420:       PetscCall(PetscLogStageRegister("Step2a - PCSetUp", &stage));
421:       PetscCall(PetscLogStagePush(stage));
422:       PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
423:       PetscCall(KSPSetOperators(ksp, A, A));
424:       PetscCall(KSPSetFromOptions(ksp));
425:       PetscCall(KSPGetPC(ksp, &pc));
426:       PetscCall(PetscTime(&time_start));
427:       PetscCall(PCSetUp(pc));
428:       PetscCall(PetscTime(&time_mid1));
429:       PetscCall(PetscLogStagePop());
430:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Step2b - running KSPSolve()...\n"));
431:       PetscCall(PetscLogStageRegister("Step2b - KSPSolve", &stage));
432:       PetscCall(PetscLogStagePush(stage));
433:       PetscCall(PetscTime(&time_mid2));
434:       PetscCall(KSPSolve(ksp, b, x));
435:       PetscCall(PetscTime(&time_end));
436:       PetscCall(PetscLogStagePop());
437:     }

439:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
440:     /*  Calculate error norm                                               */
441:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
442:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Step3  - calculating error norm...\n"));
443:     PetscCall(PetscLogStageRegister("Step3  - Error norm", &stage));
444:     PetscCall(PetscLogStagePush(stage));
445:     PetscCall(VecAXPY(x, -1.0, u));
446:     PetscCall(VecNorm(x, NORM_2, &norm));
447:     PetscCall(PetscLogStagePop());

449:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
450:     /*  Summary                                                            */
451:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
452:     PetscCall(KSPGetIterationNumber(ksp, &ksp_its));
453:     if (printTiming) {
454:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n%-15s%-1.3e\n", "Error norm:", (double)norm));
455:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%-15s%-3" PetscInt_FMT "\n", "KSP iters:", ksp_its));
456:       if (user.splitksp) {
457:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%-15s%-7.5f seconds\n", "PCSetUp:", time_mid1 - time_start));
458:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%-15s%-7.5f seconds\n", "KSPSolve:", time_end - time_mid2));
459:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%-15s%-7.5f seconds\n", "Total Solve:", time_end - time_start));
460:       } else {
461:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%-15s%-7.5f seconds\n", "KSPSolve:", time_end - time_start));
462:       }
463:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%-15s%-1.3e DoFs/sec\n", "FOM:", user.dim / (time_end - time_start))); /* figure of merit */
464:     }
465:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "===========================================\n"));
466:   }
467:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
468:   /*  Free up memory                                                     */
469:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
470:   if (!user.matmult) {
471:     PetscCall(KSPDestroy(&ksp));
472:     PetscCall(VecDestroy(&x));
473:   }
474:   PetscCall(VecDestroy(&u));
475:   PetscCall(VecDestroy(&b));
476:   PetscCall(MatDestroy(&A));
477:   PetscCall(PetscFinalize());
478:   return 0;
479: }

481: /*TEST

483:   testset:
484:     args: -print_timing false -matmult -its 10 -n 8
485:     nsize: {{1 3}}
486:     output_file: output/bench_kspsolve_matmult.out

488:     test:
489:       suffix: matmult

491:     test:
492:       suffix: hip_matmult
493:       requires: hip
494:       args: -mat_type aijhipsparse

496:     test:
497:       suffix: cuda_matmult
498:       requires: cuda
499:       args: -mat_type aijcusparse

501:     test:
502:       suffix: kok_matmult
503:       requires: kokkos_kernels
504:       args: -mat_type aijkokkos

506:   testset:
507:     args: -print_timing false -its 10 -n 8
508:     nsize: {{3 5}}
509:     output_file: output/bench_kspsolve_ksp.out

511:     test:
512:       suffix: ksp

514:     test:
515:       suffix: nbr
516:       requires: defined(PETSC_HAVE_MPI_PERSISTENT_NEIGHBORHOOD_COLLECTIVES)
517:       args: -sf_type neighbor -sf_neighbor_persistent {{0 1}}

519:     test:
520:       suffix: hip_ksp
521:       requires: hip
522:       args: -mat_type aijhipsparse

524:     test:
525:       suffix: cuda_ksp
526:       requires: cuda
527:       args: -mat_type aijcusparse

529:     test:
530:       suffix: kok_ksp
531:       requires: kokkos_kernels
532:       args: -mat_type aijkokkos

534:     test:
535:       suffix: kok_nbr
536:       requires: kokkos_kernels defined(PETSC_HAVE_MPI_PERSISTENT_NEIGHBORHOOD_COLLECTIVES)
537:       args: -mat_type aijkokkos -sf_type neighbor -sf_neighbor_persistent {{0 1}}

539: TEST*/