1 static char help[] = "Test PCFailedReason.\n\n"; 2 3 #include <petscksp.h> 4 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; 16 17 PetscFunctionBeginUser; 18 PetscCall(PetscInitialize(&argc, &args, NULL, 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)); 23 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; 27 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)); 35 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)); 57 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)); 67 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 } 81 82 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 83 Factorize second matrix 84 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 85 PetscCall(MatShift(A, reverse ? alpha - beta : beta - alpha)); 86 PetscCall(KSPSetOperators(ksp, A, A)); 87 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 } 100 101 /* 102 Free work space. 103 */ 104 PetscCall(MatDestroy(&A)); 105 PetscCall(KSPDestroy(&ksp)); 106 107 PetscCall(PetscFinalize()); 108 return 0; 109 } 110 111 /*TEST 112 113 test: 114 args: -reverse 115 116 test: 117 suffix: 2 118 args: -reverse -pc_type cholesky 119 120 TEST*/ 121