Actual source code: ex51.c
1: static char help[] = "Test PCFailedReason.\n\n";
3: #include <petscksp.h>
5: int main(int argc, char **args)
6: {
7: Mat A; /* linear system matrix */
8: KSP ksp; /* linear solver context */
9: PC pc; /* preconditioner context */
10: PetscInt i, n = 10, col[3];
11: PetscMPIInt size;
12: PetscScalar value[3], alpha, beta, sx;
13: PetscBool reverse = PETSC_FALSE;
14: KSPConvergedReason reason;
15: PCFailedReason pcreason;
17: PetscFunctionBeginUser;
18: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
19: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
20: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
21: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
22: PetscCall(PetscOptionsGetBool(NULL, NULL, "-reverse", &reverse, NULL));
24: sx = PetscSinReal(n * PETSC_PI / 2 / (n + 1));
25: alpha = 4.0 * sx * sx; /* alpha is the largest eigenvalue of the matrix */
26: beta = 4.0;
28: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
29: Create the matrix
30: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
31: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
32: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
33: PetscCall(MatSetFromOptions(A));
34: PetscCall(MatSetUp(A));
36: value[0] = -1.0;
37: value[1] = 2.0;
38: value[2] = -1.0;
39: for (i = 1; i < n - 1; i++) {
40: col[0] = i - 1;
41: col[1] = i;
42: col[2] = i + 1;
43: PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
44: }
45: i = n - 1;
46: col[0] = n - 2;
47: col[1] = n - 1;
48: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
49: i = 0;
50: col[0] = 0;
51: col[1] = 1;
52: value[0] = 2.0;
53: value[1] = -1.0;
54: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
55: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
56: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
58: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59: Create the linear solver and set various options
60: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
61: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
62: PetscCall(KSPSetOperators(ksp, A, A));
63: PetscCall(MatShift(A, reverse ? -alpha : -beta));
64: PetscCall(KSPGetPC(ksp, &pc));
65: PetscCall(PCSetType(pc, PCLU));
66: PetscCall(KSPSetFromOptions(ksp));
68: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: Factorize first matrix
70: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
71: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "First matrix\n"));
72: PetscCall(KSPSetUp(ksp));
73: PetscCall(KSPGetConvergedReason(ksp, &reason));
74: if (reason < 0) {
75: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp)), PETSC_VIEWER_DEFAULT));
76: PetscCall(KSPConvergedReasonView(ksp, PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp))));
77: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp))));
78: } else {
79: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Success!\n"));
80: }
82: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: Factorize second matrix
84: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85: PetscCall(MatShift(A, reverse ? alpha - beta : beta - alpha));
86: PetscCall(KSPSetOperators(ksp, A, A));
88: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Second matrix\n"));
89: PetscCall(KSPSetUp(ksp));
90: PetscCall(KSPGetConvergedReason(ksp, &reason));
91: if (reason < 0) {
92: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp)), PETSC_VIEWER_DEFAULT));
93: PetscCall(KSPConvergedReasonView(ksp, PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp))));
94: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp))));
95: } else {
96: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Success!\n"));
97: PetscCall(PCGetFailedReason(pc, &pcreason));
98: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "PC failed reason is %s\n", PCFailedReasons[pcreason]));
99: }
101: /*
102: Free work space.
103: */
104: PetscCall(MatDestroy(&A));
105: PetscCall(KSPDestroy(&ksp));
107: PetscCall(PetscFinalize());
108: return 0;
109: }
111: /*TEST
113: test:
114: args: -reverse
116: test:
117: suffix: 2
118: args: -reverse -pc_type cholesky
120: TEST*/