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