xref: /petsc/src/ksp/ksp/tests/ex51.c (revision bd35522dd00d2ebcbc9f656e2902fae240dc8904)
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,(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));
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; value[1] = 2.0; value[2] = -1.0;
37   for (i=1; i<n-1; i++) {
38     col[0] = i-1; col[1] = i; col[2] = i+1;
39     PetscCall(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES));
40   }
41   i    = n - 1; col[0] = n - 2; col[1] = n - 1;
42   PetscCall(MatSetValues(A,1,&i,2,col,value,INSERT_VALUES));
43   i    = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0;
44   PetscCall(MatSetValues(A,1,&i,2,col,value,INSERT_VALUES));
45   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
46   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
47 
48   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49                 Create the linear solver and set various options
50      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
51   PetscCall(KSPCreate(PETSC_COMM_WORLD,&ksp));
52   PetscCall(KSPSetOperators(ksp,A,A));
53   PetscCall(MatShift(A,reverse?-alpha:-beta));
54   PetscCall(KSPGetPC(ksp,&pc));
55   PetscCall(PCSetType(pc,PCLU));
56   PetscCall(KSPSetFromOptions(ksp));
57 
58   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59                       Factorize first matrix
60      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
61   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"First matrix\n"));
62   PetscCall(KSPSetUp(ksp));
63   PetscCall(KSPGetConvergedReason(ksp,&reason));
64   if (reason < 0) {
65     PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp)),PETSC_VIEWER_DEFAULT));
66     PetscCall(KSPConvergedReasonView(ksp,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp))));
67     PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp))));
68   } else {
69     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Success!\n"));
70   }
71 
72   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73                       Factorize second matrix
74      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75   PetscCall(MatShift(A,reverse?alpha-beta:beta-alpha));
76   PetscCall(KSPSetOperators(ksp,A,A));
77 
78   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Second matrix\n"));
79   PetscCall(KSPSetUp(ksp));
80   PetscCall(KSPGetConvergedReason(ksp,&reason));
81   if (reason < 0) {
82     PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp)),PETSC_VIEWER_DEFAULT));
83     PetscCall(KSPConvergedReasonView(ksp,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp))));
84     PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp))));
85   } else {
86     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Success!\n"));
87     PetscCall(PCGetFailedReason(pc,&pcreason));
88     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"PC failed reason is %s\n",PCFailedReasons[pcreason]));
89   }
90 
91   /*
92      Free work space.
93   */
94   PetscCall(MatDestroy(&A));
95   PetscCall(KSPDestroy(&ksp));
96 
97   PetscCall(PetscFinalize());
98   return 0;
99 }
100 
101 /*TEST
102 
103    test:
104       args: -reverse
105 
106    test:
107       suffix: 2
108       args: -reverse -pc_type cholesky
109 
110 TEST*/
111