1 static const char help[] = "Tests MatGetSchurComplement\n";
2 
3 #include <petscksp.h>
4 
5 /*T
6     Concepts: Mat, Schur Complement
7 T*/
8 
Create(MPI_Comm comm,Mat * inA,IS * is0,IS * is1)9 PetscErrorCode Create(MPI_Comm comm,Mat *inA,IS *is0,IS *is1)
10 {
11   PetscErrorCode ierr;
12   Mat            A;
13   PetscInt       r,rend,M;
14   PetscMPIInt    rank;
15 
16   PetscFunctionBeginUser;
17   *inA = 0;
18   ierr = MatCreate(comm,&A);CHKERRQ(ierr);
19   ierr = MatSetSizes(A,4,4,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
20   ierr = MatSetFromOptions(A);CHKERRQ(ierr);
21   ierr = MatSetUp(A);CHKERRQ(ierr);
22   ierr = MatGetOwnershipRange(A,&r,&rend);CHKERRQ(ierr);
23   ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr);
24 
25   ierr = ISCreateStride(comm,2,r,1,is0);CHKERRQ(ierr);
26   ierr = ISCreateStride(comm,2,r+2,1,is1);CHKERRQ(ierr);
27 
28   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
29 
30   {
31     PetscInt    rows[4],cols0[5],cols1[5],cols2[3],cols3[3];
32     PetscScalar RR = 1000.*rank,vals0[5],vals1[4],vals2[3],vals3[3];
33 
34     rows[0]            = r;
35     rows[1]            = r+1;
36     rows[2]            = r+2;
37     rows[3]            = r+3;
38 
39     cols0[0]           = r+0;
40     cols0[1]           = r+1;
41     cols0[2]           = r+3;
42     cols0[3]           = (r+4)%M;
43     cols0[4]           = (r+M-4)%M;
44 
45     cols1[0]           = r+1;
46     cols1[1]           = r+2;
47     cols1[2]           = (r+4+1)%M;
48     cols1[3]           = (r+M-4+1)%M;
49 
50     cols2[0]           = r;
51     cols2[1]           = r+2;
52     cols2[2]           = (r+4+2)%M;
53 
54     cols3[0]           = r+1;
55     cols3[1]           = r+3;
56     cols3[2]           = (r+4+3)%M;
57 
58     vals0[0] = RR+1.;
59     vals0[1] = RR+2.;
60     vals0[2] = RR+3.;
61     vals0[3] = RR+4.;
62     vals0[4] = RR+5.;
63 
64     vals1[0] = RR+6.;
65     vals1[1] = RR+7.;
66     vals1[2] = RR+8.;
67     vals1[3] = RR+9.;
68 
69     vals2[0] = RR+10.;
70     vals2[1] = RR+11.;
71     vals2[2] = RR+12.;
72 
73     vals3[0] = RR+13.;
74     vals3[1] = RR+14.;
75     vals3[2] = RR+15.;
76     ierr = MatSetValues(A,1,&rows[0],5,cols0,vals0,INSERT_VALUES);CHKERRQ(ierr);
77     ierr = MatSetValues(A,1,&rows[1],4,cols1,vals1,INSERT_VALUES);CHKERRQ(ierr);
78     ierr = MatSetValues(A,1,&rows[2],3,cols2,vals2,INSERT_VALUES);CHKERRQ(ierr);
79     ierr = MatSetValues(A,1,&rows[3],3,cols3,vals3,INSERT_VALUES);CHKERRQ(ierr);
80   }
81   ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
82   ierr = MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
83   *inA = A;
84   PetscFunctionReturn(0);
85 }
86 
Destroy(Mat * A,IS * is0,IS * is1)87 PetscErrorCode Destroy(Mat *A,IS *is0,IS *is1)
88 {
89   PetscErrorCode ierr;
90 
91   PetscFunctionBeginUser;
92   ierr = MatDestroy(A);CHKERRQ(ierr);
93   ierr = ISDestroy(is0);CHKERRQ(ierr);
94   ierr = ISDestroy(is1);CHKERRQ(ierr);
95   PetscFunctionReturn(0);
96 }
97 
main(int argc,char * argv[])98 int main(int argc,char *argv[])
99 {
100   PetscErrorCode ierr;
101   Mat                        A,S = NULL,Sexplicit = NULL;
102   MatSchurComplementAinvType ainv_type = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
103   IS                         is0,is1;
104 
105   ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
106   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex21","KSP");CHKERRQ(ierr);
107   ierr = PetscOptionsEnum("-mat_schur_complement_ainv_type","Type of approximation for inv(A00) used when assembling Sp = A11 - A10 inv(A00) A01","MatSchurComplementAinvType",MatSchurComplementAinvTypes,(PetscEnum)ainv_type,(PetscEnum*)&ainv_type,NULL);CHKERRQ(ierr);
108   ierr = PetscOptionsEnd();CHKERRQ(ierr);
109 
110   /* Test the Schur complement one way */
111   ierr = Create(PETSC_COMM_WORLD,&A,&is0,&is1);CHKERRQ(ierr);
112   ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
113   ierr = ISView(is0,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
114   ierr = ISView(is1,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
115   ierr = MatGetSchurComplement(A,is0,is0,is1,is1,MAT_INITIAL_MATRIX,&S,ainv_type,MAT_IGNORE_MATRIX,NULL);CHKERRQ(ierr);
116   ierr = MatSetFromOptions(S);CHKERRQ(ierr);
117   ierr = MatComputeOperator(S,MATAIJ,&Sexplicit);CHKERRQ(ierr);
118   ierr = PetscPrintf(PETSC_COMM_WORLD,"\nExplicit Schur complement of (0,0) in (1,1)\n");CHKERRQ(ierr);
119   ierr = MatView(Sexplicit,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
120   ierr = Destroy(&A,&is0,&is1);CHKERRQ(ierr);
121   ierr = MatDestroy(&S);CHKERRQ(ierr);
122   ierr = MatDestroy(&Sexplicit);CHKERRQ(ierr);
123 
124   /* And the other */
125   ierr = Create(PETSC_COMM_WORLD,&A,&is0,&is1);CHKERRQ(ierr);
126   ierr = MatGetSchurComplement(A,is1,is1,is0,is0,MAT_INITIAL_MATRIX,&S,ainv_type,MAT_IGNORE_MATRIX,NULL);CHKERRQ(ierr);
127   ierr = MatSetFromOptions(S);CHKERRQ(ierr);
128   ierr = MatComputeOperator(S,MATAIJ,&Sexplicit);CHKERRQ(ierr);
129   ierr = PetscPrintf(PETSC_COMM_WORLD,"\nExplicit Schur complement of (1,1) in (0,0)\n");CHKERRQ(ierr);
130   ierr = MatView(Sexplicit,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
131   ierr = Destroy(&A,&is0,&is1);CHKERRQ(ierr);
132   ierr = MatDestroy(&S);CHKERRQ(ierr);
133   ierr = MatDestroy(&Sexplicit);CHKERRQ(ierr);
134 
135   /* This time just the preconditioning matrix. */
136   ierr = Create(PETSC_COMM_WORLD,&A,&is0,&is1);CHKERRQ(ierr);
137   ierr = MatGetSchurComplement(A,is0,is0,is1,is1,MAT_IGNORE_MATRIX,NULL,ainv_type,MAT_INITIAL_MATRIX,&S);CHKERRQ(ierr);
138   ierr = MatSetFromOptions(S);CHKERRQ(ierr);
139   ierr = PetscPrintf(PETSC_COMM_WORLD,"\nPreconditioning Schur complement of (0,0) in (1,1)\n");CHKERRQ(ierr);
140   ierr = MatView(S,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
141   /* Modify and refresh */
142   ierr = MatShift(A,1.);CHKERRQ(ierr);
143   ierr = MatGetSchurComplement(A,is0,is0,is1,is1,MAT_IGNORE_MATRIX,NULL,ainv_type,MAT_REUSE_MATRIX,&S);CHKERRQ(ierr);
144   ierr = PetscPrintf(PETSC_COMM_WORLD,"\nAfter update\n");CHKERRQ(ierr);
145   ierr = MatView(S,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
146   ierr = Destroy(&A,&is0,&is1);CHKERRQ(ierr);
147   ierr = MatDestroy(&S);CHKERRQ(ierr);
148 
149   ierr = PetscFinalize();
150   return ierr;
151 }
152 
153 /*TEST
154   test:
155     suffix: diag_1
156     args: -mat_schur_complement_ainv_type diag
157     nsize: 1
158   test:
159     suffix: blockdiag_1
160     args: -mat_schur_complement_ainv_type blockdiag
161     nsize: 1
162   test:
163     suffix: diag_2
164     args: -mat_schur_complement_ainv_type diag
165     nsize: 2
166   test:
167     suffix: blockdiag_2
168     args: -mat_schur_complement_ainv_type blockdiag
169     nsize: 2
170   test:
171     # does not work with single because residual norm computed by GMRES recurrence formula becomes invalid
172     requires: !single
173     suffix: diag_3
174     args: -mat_schur_complement_ainv_type diag
175     nsize: 3
176   test:
177     # does not work with single because residual norm computed by GMRES recurrence formula becomes invalid
178     requires: !single
179     suffix: blockdiag_3
180     args: -mat_schur_complement_ainv_type blockdiag
181     nsize: 3
182 TEST*/
183