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