Actual source code: ex27.c

  1: static char help[] = "Reads a PETSc matrix and vector from a file and solves the normal equations.\n\n";

  3: /*
  4:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
  5:   automatically includes:
  6:      petscsys.h       - base PETSc routines   petscvec.h - vectors
  7:      petscmat.h - matrices
  8:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
  9:      petscviewer.h - viewers               petscpc.h  - preconditioners
 10: */
 11: #include <petscksp.h>
 12: #include <petscviewerhdf5.h>

 14: static PetscErrorCode VecLoadIfExists_Private(Vec b, PetscViewer fd, PetscBool *has)
 15: {
 16:   PetscBool hdf5 = PETSC_FALSE;

 18:   PetscFunctionBeginUser;
 19:   PetscCall(PetscObjectTypeCompare((PetscObject)fd, PETSCVIEWERHDF5, &hdf5));
 20:   if (hdf5) {
 21: #if defined(PETSC_HAVE_HDF5)
 22:     PetscCall(PetscViewerHDF5HasObject(fd, (PetscObject)b, has));
 23:     if (*has) PetscCall(VecLoad(b, fd));
 24: #else
 25:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "PETSc must be configured with HDF5 to use this feature");
 26: #endif
 27:   } else {
 28:     PetscErrorCode ierrp;
 29:     PetscCall(PetscPushErrorHandler(PetscReturnErrorHandler, NULL));
 30:     ierrp = VecLoad(b, fd);
 31:     PetscCall(PetscPopErrorHandler());
 32:     *has = ierrp ? PETSC_FALSE : PETSC_TRUE;
 33:   }
 34:   PetscFunctionReturn(PETSC_SUCCESS);
 35: }

 37: int main(int argc, char **args)
 38: {
 39:   KSP         ksp;                                                       /* linear solver context */
 40:   Mat         A, N;                                                      /* matrix */
 41:   Vec         x, b, r, Ab, v[2];                                         /* approx solution, RHS, residual */
 42:   PetscViewer fd;                                                        /* viewer */
 43:   char        file[PETSC_MAX_PATH_LEN]    = "";                          /* input file name */
 44:   char        file_x0[PETSC_MAX_PATH_LEN] = "";                          /* name of input file with initial guess */
 45:   char        A_name[128] = "A", b_name[128] = "b", x0_name[128] = "x0"; /* name of the matrix, RHS and initial guess */
 46:   KSPType     ksptype;
 47:   PetscBool   has;
 48:   PetscInt    its, n, m;
 49:   PetscReal   norm;
 50:   PetscBool   nonzero_guess      = PETSC_TRUE;
 51:   PetscBool   solve_normal       = PETSC_FALSE;
 52:   PetscBool   solve_augmented    = PETSC_FALSE;
 53:   PetscBool   truncate           = PETSC_FALSE;
 54:   PetscBool   explicit_transpose = PETSC_FALSE;
 55:   PetscBool   hdf5               = PETSC_FALSE;
 56:   PetscBool   test_custom_layout = PETSC_FALSE;
 57:   PetscBool   sbaij              = PETSC_FALSE;
 58:   PetscMPIInt rank, size;

 60:   PetscFunctionBeginUser;
 61:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 62:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 63:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 64:   /*
 65:      Determine files from which we read the linear system
 66:      (matrix, right-hand side and initial guess vector).
 67:   */
 68:   PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), NULL));
 69:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-truncate", &truncate, NULL));
 70:   if (!truncate) PetscCall(PetscOptionsGetString(NULL, NULL, "-f_x0", file_x0, sizeof(file_x0), NULL));
 71:   PetscCall(PetscOptionsGetString(NULL, NULL, "-A_name", A_name, sizeof(A_name), NULL));
 72:   PetscCall(PetscOptionsGetString(NULL, NULL, "-b_name", b_name, sizeof(b_name), NULL));
 73:   PetscCall(PetscOptionsGetString(NULL, NULL, "-x0_name", x0_name, sizeof(x0_name), NULL));
 74:   /*
 75:      Decide whether to solve the original system (-solve_normal 0)
 76:      or the normal equation (-solve_normal 1).
 77:   */
 78:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-solve_normal", &solve_normal, NULL));
 79:   if (!solve_normal) PetscCall(PetscOptionsGetBool(NULL, NULL, "-solve_augmented", &solve_augmented, NULL));
 80:   if (solve_augmented) {
 81:     PetscCall(PetscOptionsGetBool(NULL, NULL, "-explicit_transpose", &explicit_transpose, NULL));
 82:     PetscCall(PetscOptionsGetBool(NULL, NULL, "-sbaij", &sbaij, NULL));
 83:   }
 84:   /*
 85:      Decide whether to use the HDF5 reader.
 86:   */
 87:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-hdf5", &hdf5, NULL));
 88:   /*
 89:      Decide whether custom matrix layout will be tested.
 90:   */
 91:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_custom_layout", &test_custom_layout, NULL));

 93:   /* -----------------------------------------------------------
 94:                   Beginning of linear solver loop
 95:      ----------------------------------------------------------- */
 96:   /*
 97:      Loop through the linear solve 2 times.
 98:       - The intention here is to preload and solve a small system;
 99:         then load another (larger) system and solve it as well.
100:         This process preloads the instructions with the smaller
101:         system so that more accurate performance monitoring (via
102:         -log_view) can be done with the larger one (that actually
103:         is the system of interest).
104:   */
105:   PetscPreLoadBegin(PETSC_FALSE, "Load system");

107:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
108:                          Load system
109:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

111:   /*
112:      Open binary file.  Note that we use FILE_MODE_READ to indicate
113:      reading from this file.
114:   */
115:   if (hdf5) {
116: #if defined(PETSC_HAVE_HDF5)
117:     PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
118:     PetscCall(PetscViewerPushFormat(fd, PETSC_VIEWER_HDF5_MAT));
119: #else
120:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "PETSc must be configured with HDF5 to use this feature");
121: #endif
122:   } else {
123:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
124:   }

126:   /*
127:      Load the matrix.
128:      Matrix type is set automatically but you can override it by MatSetType() prior to MatLoad().
129:      Do that only if you really insist on the given type.
130:   */
131:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
132:   PetscCall(PetscObjectSetName((PetscObject)A, A_name));
133:   PetscCall(MatSetFromOptions(A));
134:   PetscCall(MatLoad(A, fd));
135:   if (truncate) {
136:     Mat      P, B;
137:     PetscInt M, N;
138:     PetscCall(MatGetLocalSize(A, &m, &n));
139:     PetscCall(MatGetSize(A, &M, &N));
140:     PetscCall(MatCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, m, PETSC_DECIDE, M, N / 1.5, &P));
141:     PetscCall(MatGetOwnershipRangeColumn(P, &m, &n));
142:     for (; m < n; ++m) PetscCall(MatSetValue(P, m, m, 1.0, INSERT_VALUES));
143:     PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
144:     PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
145:     PetscCall(MatShift(P, 1.0));
146:     PetscCall(MatMatMult(A, P, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &B));
147:     PetscCall(MatDestroy(&P));
148:     PetscCall(MatDestroy(&A));
149:     A = B;
150:   }
151:   if (test_custom_layout && size > 1) {
152:     /* Perturb the local sizes and create the matrix anew */
153:     PetscInt m1, n1;
154:     PetscCall(MatGetLocalSize(A, &m, &n));
155:     m = rank ? m - 1 : m + size - 1;
156:     n = (rank == size - 1) ? n + size - 1 : n - 1;
157:     PetscCall(MatDestroy(&A));
158:     PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
159:     PetscCall(PetscObjectSetName((PetscObject)A, A_name));
160:     PetscCall(MatSetSizes(A, m, n, PETSC_DECIDE, PETSC_DECIDE));
161:     PetscCall(MatSetFromOptions(A));
162:     PetscCall(MatLoad(A, fd));
163:     PetscCall(MatGetLocalSize(A, &m1, &n1));
164:     PetscCheck(m1 == m && n1 == n, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "resulting sizes differ from requested ones: %" PetscInt_FMT " %" PetscInt_FMT " != %" PetscInt_FMT " %" PetscInt_FMT, m1, n1, m, n);
165:   }
166:   PetscCall(MatGetLocalSize(A, &m, &n));

168:   /*
169:      Load the RHS vector if it is present in the file, otherwise use a vector of all ones.
170:   */
171:   PetscCall(MatCreateVecs(A, &x, &b));
172:   PetscCall(PetscObjectSetName((PetscObject)b, b_name));
173:   PetscCall(VecSetFromOptions(b));
174:   PetscCall(VecLoadIfExists_Private(b, fd, &has));
175:   if (!has) {
176:     PetscScalar one = 1.0;
177:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Failed to load RHS, so use a vector of all ones.\n"));
178:     PetscCall(VecSetFromOptions(b));
179:     PetscCall(VecSet(b, one));
180:   }

182:   /*
183:      Load the initial guess vector if it is present in the file, otherwise use a vector of all zeros.
184:   */
185:   PetscCall(PetscObjectSetName((PetscObject)x, x0_name));
186:   PetscCall(VecSetFromOptions(x));
187:   if (!truncate) {
188:     /* load file_x0 if it is specified, otherwise try to reuse file */
189:     if (file_x0[0]) {
190:       PetscCall(PetscViewerDestroy(&fd));
191:       if (hdf5) {
192: #if defined(PETSC_HAVE_HDF5)
193:         PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, file_x0, FILE_MODE_READ, &fd));
194: #endif
195:       } else {
196:         PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file_x0, FILE_MODE_READ, &fd));
197:       }
198:     }
199:     PetscCall(VecLoadIfExists_Private(x, fd, &has));
200:   } else has = PETSC_FALSE;
201:   if (truncate || !has) {
202:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Failed to load initial guess, so use a vector of all zeros.\n"));
203:     PetscCall(VecSet(x, 0.0));
204:     nonzero_guess = PETSC_FALSE;
205:   }
206:   PetscCall(PetscViewerDestroy(&fd));

208:   PetscCall(VecDuplicate(x, &Ab));

210:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
211:                     Setup solve for system
212:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

214:   /*
215:      Conclude profiling last stage; begin profiling next stage.
216:   */
217:   PetscPreLoadStage("KSPSetUp");

219:   PetscCall(MatCreateNormalHermitian(A, &N));
220:   PetscCall(MatMultHermitianTranspose(A, b, Ab));

222:   /*
223:      Create linear solver; set operators; set runtime options.
224:   */
225:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));

227:   if (solve_normal) {
228:     PetscCall(KSPSetOperators(ksp, N, N));
229:   } else if (solve_augmented) {
230:     Mat       array[4], C, S;
231:     Vec       view;
232:     PetscInt  M, n;
233:     PetscReal diag;

235:     PetscCall(MatDestroy(&N));
236:     PetscCall(MatGetSize(A, &M, NULL));
237:     PetscCall(MatGetLocalSize(A, NULL, &n));
238:     PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD, m, m, M, M, -1.0, array));
239:     array[1] = A;
240:     if (!explicit_transpose) PetscCall(MatCreateHermitianTranspose(A, array + 2));
241:     else PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, array + 2));
242:     PetscCall(PetscOptionsGetReal(NULL, NULL, "-nonzero_A11", &diag, &has));
243:     if (has) PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD, n, n, PETSC_DECIDE, PETSC_DECIDE, diag, array + 3));
244:     else array[3] = NULL;
245:     PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, array, &C));
246:     if (!sbaij) PetscCall(MatNestSetVecType(C, VECNEST));
247:     PetscCall(MatCreateVecs(C, v + 1, v));
248:     if (!sbaij) {
249:       PetscCall(VecNestGetSubVec(v[0], 0, &view));
250:       PetscCall(VecCopy(b, view));
251:       PetscCall(VecNestGetSubVec(v[1], 1, &view));
252:       PetscCall(VecCopy(x, view));
253:       PetscCall(KSPSetOperators(ksp, C, C));
254:     } else {
255:       const PetscScalar *read;
256:       PetscScalar       *write;
257:       PetscCall(VecGetArrayRead(b, &read));
258:       PetscCall(VecGetArrayWrite(v[0], &write));
259:       for (PetscInt i = 0; i < m; ++i) write[i] = read[i];
260:       PetscCall(VecRestoreArrayWrite(v[0], &write));
261:       PetscCall(VecRestoreArrayRead(b, &read));
262:       PetscCall(VecGetArrayRead(x, &read));
263:       PetscCall(VecGetArrayWrite(v[1], &write));
264:       for (PetscInt i = 0; i < n; ++i) write[m + i] = read[i];
265:       PetscCall(VecRestoreArrayWrite(v[1], &write));
266:       PetscCall(VecRestoreArrayRead(x, &read));
267:       PetscCall(MatConvert(C, MATSBAIJ, MAT_INITIAL_MATRIX, &S));
268:       PetscCall(KSPSetOperators(ksp, S, S));
269:       PetscCall(MatDestroy(&S));
270:     }
271:     PetscCall(MatDestroy(&C));
272:     PetscCall(MatDestroy(array));
273:     PetscCall(MatDestroy(array + 2));
274:     PetscCall(MatDestroy(array + 3));
275:   } else {
276:     PC pc;
277:     PetscCall(KSPSetType(ksp, KSPLSQR));
278:     PetscCall(KSPGetPC(ksp, &pc));
279:     PetscCall(PCSetType(pc, PCNONE));
280:     PetscCall(KSPSetOperators(ksp, A, N));
281:   }
282:   PetscCall(KSPSetInitialGuessNonzero(ksp, nonzero_guess));
283:   PetscCall(KSPSetFromOptions(ksp));

285:   /*
286:      Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
287:      enable more precise profiling of setting up the preconditioner.
288:      These calls are optional, since both will be called within
289:      KSPSolve() if they haven't been called already.
290:   */
291:   PetscCall(KSPSetUp(ksp));
292:   PetscCall(KSPSetUpOnBlocks(ksp));

294:   /*
295:                          Solve system
296:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

298:   /*
299:      Begin profiling next stage
300:   */
301:   PetscPreLoadStage("KSPSolve");

303:   /*
304:      Solve linear system
305:   */
306:   if (solve_normal) {
307:     PetscCall(KSPSolve(ksp, Ab, x));
308:   } else if (solve_augmented) {
309:     KSP      *subksp;
310:     PC        pc;
311:     Mat       C;
312:     Vec       view;
313:     PetscBool flg;

315:     PetscCall(KSPGetPC(ksp, &pc));
316:     PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &flg));
317:     if (flg) {
318:       PetscCall(KSPGetOperators(ksp, &C, NULL));
319:       PetscCall(PCFieldSplitGetSubKSP(pc, NULL, &subksp));
320:       PetscCall(KSPGetPC(subksp[1], &pc));
321:       PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCHPDDM, &flg));
322:       if (flg) {
323: #if defined(PETSC_HAVE_HPDDM) && defined(PETSC_HAVE_DYNAMIC_LIBRARIES) && defined(PETSC_USE_SHARED_LIBRARIES)
324:         Mat aux, S, **array;
325:         IS  is;

327:         PetscCall(MatCreate(PETSC_COMM_SELF, &aux));
328:         PetscCall(ISCreate(PETSC_COMM_SELF, &is));
329:         PetscCall(PCHPDDMSetAuxiliaryMat(pc, is, aux, NULL, NULL)); /* dummy objects just to cover corner cases in PCSetUp() */
330:         PetscCall(ISDestroy(&is));
331:         PetscCall(MatDestroy(&aux));
332:         PetscCall(MatNestGetSubMats(C, NULL, NULL, &array));
333:         PetscCall(MatCreateSchurComplement(array[0][0], array[0][0], array[0][1], array[1][0], array[1][1], &S));
334:         PetscCall(MatSetOptionsPrefix(S, "fieldsplit_1_"));
335:         PetscCall(KSPSetOperators(subksp[1], S, S));
336:         PetscCall(MatDestroy(&S));
337:         PetscCall(PCSetFromOptions(pc));
338: #endif
339:       }
340:       PetscCall(PetscFree(subksp));
341:     }
342:     PetscCall(KSPSolve(ksp, v[0], v[1]));
343:     if (!sbaij) {
344:       PetscCall(VecNestGetSubVec(v[1], 1, &view));
345:       PetscCall(VecCopy(view, x));
346:     } else {
347:       const PetscScalar *read;
348:       PetscScalar       *write;
349:       PetscCall(MatGetLocalSize(A, &m, &n));
350:       PetscCall(VecGetArrayRead(v[1], &read));
351:       PetscCall(VecGetArrayWrite(x, &write));
352:       for (PetscInt i = 0; i < n; ++i) write[i] = read[m + i];
353:       PetscCall(VecRestoreArrayWrite(x, &write));
354:       PetscCall(VecRestoreArrayRead(v[1], &read));
355:     }
356:   } else {
357:     PetscCall(KSPSolve(ksp, b, x));
358:   }
359:   PetscCall(PetscObjectSetName((PetscObject)x, "x"));

361:   /*
362:       Conclude profiling this stage
363:    */
364:   PetscPreLoadStage("Cleanup");

366:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
367:           Check error, print output, free data structures.
368:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

370:   /*
371:      Check error
372:   */
373:   PetscCall(VecDuplicate(b, &r));
374:   PetscCall(MatMult(A, x, r));
375:   PetscCall(VecAXPY(r, -1.0, b));
376:   PetscCall(VecNorm(r, NORM_2, &norm));
377:   PetscCall(KSPGetIterationNumber(ksp, &its));
378:   PetscCall(KSPGetType(ksp, &ksptype));
379:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSP type: %s\n", ksptype));
380:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %3" PetscInt_FMT "\n", its));
381:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g\n", (double)norm));

383:   /*
384:      Free work space.  All PETSc objects should be destroyed when they
385:      are no longer needed.
386:   */
387:   PetscCall(MatDestroy(&A));
388:   PetscCall(VecDestroy(&b));
389:   PetscCall(MatDestroy(&N));
390:   PetscCall(VecDestroy(&Ab));
391:   PetscCall(VecDestroy(&r));
392:   PetscCall(VecDestroy(&x));
393:   if (solve_augmented) {
394:     PetscCall(VecDestroy(v));
395:     PetscCall(VecDestroy(v + 1));
396:   }
397:   PetscCall(KSPDestroy(&ksp));
398:   PetscPreLoadEnd();
399:   /* -----------------------------------------------------------
400:                       End of linear solver loop
401:      ----------------------------------------------------------- */

403:   PetscCall(PetscFinalize());
404:   return 0;
405: }

407: /*TEST

409:    test:
410:       suffix: 1
411:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
412:       args: -f ${DATAFILESPATH}/matrices/medium -ksp_view -ksp_monitor_short -ksp_max_it 100 -solve_normal

414:    test:
415:       suffix: 2
416:       nsize: 2
417:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
418:       args: -f ${DATAFILESPATH}/matrices/shallow_water1 -ksp_view -ksp_monitor_short -ksp_max_it 100 -solve_normal -pc_type none

420:    # Test handling failing VecLoad without abort
421:    testset:
422:      requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
423:      args: -ksp_type cg -ksp_view -ksp_converged_reason -ksp_monitor_short -ksp_max_it 10
424:      test:
425:         suffix: 3
426:         nsize: {{1 2}separate output}
427:         args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
428:         args: -f_x0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_x0
429:      test:
430:         suffix: 3a
431:         nsize: {{1 2}separate output}
432:         args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
433:         args: -f_x0 NONEXISTING_FILE
434:      test:
435:         suffix: 3b
436:         nsize: {{1 2}separate output}
437:         args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_with_x0 # this file includes all A, b and x0
438:      test:
439:         # Load square matrix, RHS and initial guess from HDF5 (Version 7.3 MAT-File)
440:         suffix: 3b_hdf5
441:         requires: hdf5 defined(PETSC_HDF5_HAVE_ZLIB)
442:         nsize: {{1 2}separate output}
443:         args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_with_x0.mat -hdf5

445:    # Test least-square algorithms
446:    testset:
447:      requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
448:      args: -f ${DATAFILESPATH}/matrices/rectangular_ultrasound_4889x841
449:      test:
450:         suffix: 4
451:         nsize: {{1 2 4}}
452:         args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
453:         args: -solve_normal -ksp_type cg
454:      test:
455:         suffix: 4a
456:         nsize: {{1 2 4}}
457:         args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
458:         args: -ksp_type {{cgls lsqr}separate output}
459:      test:
460:         # Test KSPLSQR-specific options
461:         suffix: 4b
462:         nsize: 2
463:         args: -ksp_converged_reason -ksp_rtol 1e-3 -ksp_max_it 200 -ksp_view
464:         args: -ksp_type lsqr -ksp_convergence_test lsqr -ksp_lsqr_monitor -ksp_lsqr_compute_standard_error -ksp_lsqr_exact_mat_norm {{0 1}separate output}
465:      test:
466:         suffix: 4c
467:         nsize: 4
468:         requires: hpddm slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
469:         filter: grep -v "shared subdomain KSP between SLEPc and PETSc" | grep -v "total: nonzeros="
470:         args: -ksp_converged_reason -ksp_rtol 1e-5 -ksp_max_it 100 -ksp_view
471:         args: -ksp_type lsqr -pc_type hpddm -pc_hpddm_define_subdomains -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_levels_1_st_share_sub_ksp {{false true}shared output}
472:         args: -pc_hpddm_levels_1_pc_asm_sub_mat_type aij -pc_hpddm_levels_1_pc_asm_type basic -pc_hpddm_levels_1_sub_pc_type cholesky
473:      test:
474:         suffix: 4d
475:         nsize: 4
476:         requires: hpddm slepc suitesparse defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
477:         filter: grep -v "shared subdomain KSP between SLEPc and PETSc"
478:         args: -ksp_converged_reason -ksp_rtol 1e-5 -ksp_max_it 100 -ksp_view
479:         args: -ksp_type lsqr -pc_type hpddm -pc_hpddm_define_subdomains -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_levels_1_st_share_sub_ksp {{false true}shared output} -pc_hpddm_levels_1_st_pc_type qr
480:         args: -pc_hpddm_levels_1_pc_asm_sub_mat_type normalh -pc_hpddm_levels_1_pc_asm_type basic -pc_hpddm_levels_1_sub_pc_type qr
481:      test:
482:         suffix: 4e
483:         nsize: 4
484:         requires: hpddm slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
485:         args: -solve_augmented -ksp_type gmres
486:         args: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition self -fieldsplit_0_pc_type jacobi -fieldsplit_ksp_type preonly
487:         args: -prefix_push fieldsplit_1_ -pc_type hpddm -pc_hpddm_schur_precondition least_squares -pc_hpddm_define_subdomains -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_levels_1_st_share_sub_ksp -pc_hpddm_levels_1_sub_pc_type cholesky -prefix_pop -fieldsplit_1_mat_schur_complement_ainv_type {{lump blockdiag}shared output}
488:      test:
489:         suffix: 4f
490:         nsize: 4
491:         requires: hpddm slepc suitesparse defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
492:         filter: sed -e "s/(1,0) : type=mpiaij/(1,0) : type=transpose/g" -e "s/hermitiantranspose/transpose/g"
493:         args: -solve_augmented -ksp_type gmres -ksp_view -explicit_transpose {{false true}shared output}
494:         args: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition self -fieldsplit_0_pc_type jacobi -fieldsplit_ksp_type preonly
495:         args: -prefix_push fieldsplit_1_ -pc_type hpddm -pc_hpddm_schur_precondition least_squares -pc_hpddm_define_subdomains -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_levels_1_st_share_sub_ksp -pc_hpddm_levels_1_sub_pc_type qr -prefix_pop
496:      test:
497:         suffix: 4f_nonzero
498:         nsize: 4
499:         requires: hpddm slepc suitesparse defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
500:         args: -solve_augmented -nonzero_A11 {{0.0 1e-14}shared output} -ksp_type gmres
501:         args: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition self -fieldsplit_0_pc_type jacobi -fieldsplit_ksp_type preonly
502:         args: -prefix_push fieldsplit_1_ -pc_type hpddm -pc_hpddm_schur_precondition least_squares -pc_hpddm_define_subdomains -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_levels_1_st_share_sub_ksp -pc_hpddm_levels_1_sub_pc_type qr -prefix_pop
503:      test:
504:         suffix: 4f_nonzero_shift
505:         nsize: 4
506:         output_file: output/ex27_4f_nonzero.out
507:         requires: hpddm slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
508:         args: -solve_augmented -nonzero_A11 {{0.0 1e-6}shared output} -ksp_type gmres
509:         args: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition self -fieldsplit_0_pc_type jacobi -fieldsplit_ksp_type preonly
510:         args: -prefix_push fieldsplit_1_ -pc_type hpddm -pc_hpddm_schur_precondition least_squares -pc_hpddm_define_subdomains -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_levels_1_st_share_sub_ksp -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_levels_1_eps_gen_non_hermitian -prefix_pop
511:      test:
512:         suffix: 4g
513:         nsize: 4
514:         requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
515:         args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
516:         args: -ksp_type lsqr -pc_type hypre
517:      test:
518:         suffix: 4h
519:         nsize: {{1 4}}
520:         args: -solve_augmented -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition self -pc_fieldsplit_detect_saddle_point -sbaij true -ksp_type fgmres

522:    test:
523:       # Load rectangular matrix from HDF5 (Version 7.3 MAT-File)
524:       suffix: 4a_lsqr_hdf5
525:       nsize: {{1 2 4 8}}
526:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) hdf5 defined(PETSC_HDF5_HAVE_ZLIB)
527:       args: -f ${DATAFILESPATH}/matrices/matlab/rectangular_ultrasound_4889x841.mat -hdf5
528:       args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
529:       args: -ksp_type lsqr
530:       args: -test_custom_layout {{0 1}}

532:    # Test for correct cgls convergence reason
533:    test:
534:       suffix: 5
535:       nsize: 1
536:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
537:       args: -f ${DATAFILESPATH}/matrices/rectangular_ultrasound_4889x841
538:       args: -ksp_converged_reason -ksp_rtol 1e-2 -ksp_max_it 100
539:       args: -ksp_type cgls

541:    # Load a matrix, RHS and solution from HDF5 (Version 7.3 MAT-File). Test immediate convergence.
542:    testset:
543:      nsize: {{1 2 4 8}}
544:      requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) hdf5 defined(PETSC_HDF5_HAVE_ZLIB)
545:      args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 10
546:      args: -ksp_type lsqr
547:      args: -test_custom_layout {{0 1}}
548:      args: -hdf5 -x0_name x
549:      test:
550:         suffix: 6_hdf5
551:         args: -f ${DATAFILESPATH}/matrices/matlab/small.mat
552:      test:
553:         suffix: 6_hdf5_rect
554:         args: -f ${DATAFILESPATH}/matrices/matlab/small_rect.mat
555:      test:
556:         suffix: 6_hdf5_dense
557:         args: -f ${DATAFILESPATH}/matrices/matlab/small_dense.mat -mat_type dense
558:      test:
559:         suffix: 6_hdf5_rect_dense
560:         args: -f ${DATAFILESPATH}/matrices/matlab/small_rect_dense.mat -mat_type dense

562:    # Test correct handling of local dimensions in PCApply
563:    testset:
564:      requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
565:      requires: hdf5 defined(PETSC_HDF5_HAVE_ZLIB)
566:      nsize: 3
567:      suffix: 7
568:      args: -f ${DATAFILESPATH}/matrices/matlab/small.mat -hdf5 -test_custom_layout 1 -ksp_type lsqr -pc_type jacobi

570:    # Test complex matrices
571:    testset:
572:      requires: double complex !defined(PETSC_USE_64BIT_INDICES)
573:      args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/nh-complex-int32-float64
574:      output_file: output/ex27_8.out
575:      filter: grep -v "KSP type"
576:      test:
577:        suffix: 8
578:        args: -solve_normal 0 -ksp_type {{lsqr cgls}}
579:      test:
580:        suffix: 8_normal
581:        args: -solve_normal 1 -ksp_type {{cg bicg}}

583:    testset:
584:      requires: double suitesparse !defined(PETSC_USE_64BIT_INDICES)
585:      args: -solve_normal {{0 1}shared output} -pc_type qr
586:      output_file: output/ex27_9.out
587:      filter: grep -v "KSP type"
588:      test:
589:        suffix: 9_real
590:        requires: !complex
591:        args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/ns-real-int32-float64
592:      test:
593:        suffix: 9_complex
594:        requires: complex
595:        args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/nh-complex-int32-float64

597:    test:
598:      suffix: 10
599:      requires: !complex double suitesparse !defined(PETSC_USE_64BIT_INDICES)
600:      nsize: 2
601:      args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/ns-real-int32-float64 -pc_type bjacobi -sub_pc_type qr

603:    test:
604:      suffix: 11
605:      nsize: 4
606:      requires: datafilespath double complex !defined(PETSC_USE_64BIT_INDICES) hpddm slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
607:      args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -truncate
608:      args: -ksp_converged_reason -ksp_rtol 1e-5 -ksp_max_it 100
609:      args: -ksp_type lsqr -pc_type hpddm -pc_hpddm_define_subdomains -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_levels_1_eps_threshold_absolute 1e-6
610:      args: -pc_hpddm_levels_1_pc_asm_sub_mat_type aij -pc_hpddm_levels_1_pc_asm_type basic -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_coarse_pc_type lu
611:      filter: sed -e "s/ 10/ 9/g"

613: TEST*/