Actual source code: ex42.c

  1: static char help[] = "Solves a linear system in parallel with MINRES.\n\n";

  3: #include <petscksp.h>

  5: int main(int argc, char **args)
  6: {
  7:   Vec         x, b; /* approx solution, RHS */
  8:   Mat         A;    /* linear system matrix */
  9:   KSP         ksp;  /* linear solver context */
 10:   PC          pc;   /* preconditioner */
 11:   PetscScalar v = 0.0;
 12:   PetscInt    Ii, Istart, Iend, m = 11;
 13:   PetscBool   consistent = PETSC_TRUE;

 15:   PetscFunctionBeginUser;
 16:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 17:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 18:   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-vv", &v, NULL));
 19:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-consistent", &consistent, NULL));

 21:   /* Create parallel diagonal matrix */
 22:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 23:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, m));
 24:   PetscCall(MatSetFromOptions(A));
 25:   PetscCall(MatMPIAIJSetPreallocation(A, 1, NULL, 1, NULL));
 26:   PetscCall(MatSeqAIJSetPreallocation(A, 1, NULL));
 27:   PetscCall(MatSetUp(A));
 28:   PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));

 30:   for (Ii = Istart; Ii < Iend; Ii++) {
 31:     PetscScalar vv = (PetscReal)Ii + 1;
 32:     PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &vv, INSERT_VALUES));
 33:   }

 35:   /* Make A singular or indefinite */
 36:   Ii = m - 1; /* last diagonal entry */
 37:   PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
 38:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 39:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 41:   PetscCall(MatCreateVecs(A, &x, &b));
 42:   if (consistent) {
 43:     PetscCall(VecSet(x, 1.0));
 44:     PetscCall(MatMult(A, x, b));
 45:     PetscCall(VecSet(x, 0.0));
 46:   } else {
 47:     PetscCall(VecSet(b, 1.0));
 48:   }

 50:   /* Create linear solver context */
 51:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 52:   PetscCall(KSPSetOperators(ksp, A, A));
 53:   PetscCall(KSPSetType(ksp, KSPMINRES));
 54:   PetscCall(KSPGetPC(ksp, &pc));
 55:   PetscCall(PCSetType(pc, PCNONE));
 56:   PetscCall(KSPSetFromOptions(ksp));
 57:   PetscCall(KSPSolve(ksp, b, x));

 59:   /* Test reuse */
 60:   PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
 61:   PetscCall(VecSet(x, 0.0));
 62:   PetscCall(KSPSolve(ksp, b, x));

 64:   PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));

 66:   /* Free work space. */
 67:   PetscCall(KSPDestroy(&ksp));
 68:   PetscCall(VecDestroy(&x));
 69:   PetscCall(VecDestroy(&b));
 70:   PetscCall(MatDestroy(&A));

 72:   PetscCall(PetscFinalize());
 73:   return 0;
 74: }

 76: /*TEST

 78:    test:
 79:       args: -ksp_converged_reason

 81:    test:
 82:       suffix: 2
 83:       nsize: 3
 84:       args: -ksp_converged_reason

 86:    test:
 87:       suffix: minres_qlp
 88:       args: -ksp_converged_reason -ksp_minres_qlp -ksp_minres_monitor

 90:    test:
 91:       suffix: minres_qlp_nonconsistent
 92:       args: -ksp_converged_reason -ksp_minres_qlp -ksp_minres_monitor -consistent 0

 94:    test:
 95:       suffix: minres_neg_curve
 96:       args: -ksp_converged_neg_curve -vv -1 -ksp_converged_reason -ksp_minres_qlp {{0 1}}

 98:    test:
 99:       suffix: cg_neg_curve
100:       args: -ksp_converged_neg_curve -vv -1 -ksp_converged_reason -ksp_type {{cg stcg}}
101:       requires: !complex

103: TEST*/