Actual source code: ex21.c
1: static const char help[] = "Tests MatGetSchurComplement\n";
3: #include petscksp.h
8: PetscErrorCode Create(MPI_Comm comm,Mat *inA,IS *is0,IS *is1)
9: {
11: Mat A;
12: PetscInt r,rend,M;
13: PetscMPIInt rank;
16: *inA = 0;
17: MatCreate(comm,&A);
18: MatSetSizes(A,4,4,PETSC_DETERMINE,PETSC_DETERMINE);
19: MatSetFromOptions(A);
20: MatGetOwnershipRange(A,&r,&rend);
21: MatGetSize(A,&M,PETSC_NULL);
23: ISCreateStride(comm,2,r,1,is0);
24: ISCreateStride(comm,2,r+2,1,is1);
26: MPI_Comm_rank(comm,&rank);
28: {
29: PetscInt
30: rows[] = {r,r+1,r+2,r+3},
31: cols0[] = {r+0,r+1,r+3,(r+4)%M,(r+M-4)%M},
32: cols1[] = {r+1,r+2,(r+4+1)%M,(r+M-4+1)%M},
33: cols2[] = {r,r+2,(r+4+2)%M},
34: cols3[] = {r+1,r+3,(r+4+3)%M};
35: PetscScalar RR = 1000*rank,
36: vals0[] = {RR+1,RR+2,RR+3,RR+4,RR+5},
37: vals1[] = {RR+6,RR+7,RR+8,RR+9},
38: vals2[] = {RR+10,RR+11,RR+12},
39: vals3[] = {RR+13,RR+14,RR+15};
40: MatSetValues(A,1,&rows[0],5,cols0,vals0,INSERT_VALUES);
41: MatSetValues(A,1,&rows[1],4,cols1,vals1,INSERT_VALUES);
42: MatSetValues(A,1,&rows[2],3,cols2,vals2,INSERT_VALUES);
43: MatSetValues(A,1,&rows[3],3,cols3,vals3,INSERT_VALUES);
44: }
45: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
46: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
47: *inA = A;
48: return(0);
49: }
53: PetscErrorCode Destroy(Mat A,IS is0,IS is1)
54: {
58: MatDestroy(A);
59: ISDestroy(is0);
60: ISDestroy(is1);
61: return(0);
62: }
66: int main(int argc,char *argv[])
67: {
69: Mat A,S,Sexplicit;
70: IS is0,is1;
72: PetscInitialize(&argc,&argv,0,help);
74: /* Test the Schur complement one way */
75: Create(PETSC_COMM_WORLD,&A,&is0,&is1);
76: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
77: ISView(is0,PETSC_VIEWER_STDOUT_WORLD);
78: ISView(is1,PETSC_VIEWER_STDOUT_WORLD);
79: MatGetSchurComplement(A,is0,is0,is1,is1,MAT_INITIAL_MATRIX,&S,MAT_IGNORE_MATRIX,PETSC_NULL);
80: MatComputeExplicitOperator(S,&Sexplicit);
81: PetscPrintf(PETSC_COMM_WORLD,"\nExplicit Schur complement of (0,0) in (1,1)\n");
82: MatView(Sexplicit,PETSC_VIEWER_STDOUT_WORLD);
83: Destroy(A,is0,is1);
84: MatDestroy(S);
85: MatDestroy(Sexplicit);
87: /* And the other */
88: Create(PETSC_COMM_WORLD,&A,&is0,&is1);
89: MatGetSchurComplement(A,is1,is1,is0,is0,MAT_INITIAL_MATRIX,&S,MAT_IGNORE_MATRIX,PETSC_NULL);
90: MatComputeExplicitOperator(S,&Sexplicit);
91: PetscPrintf(PETSC_COMM_WORLD,"\nExplicit Schur complement of (1,1) in (0,0)\n");
92: MatView(Sexplicit,PETSC_VIEWER_STDOUT_WORLD);
93: Destroy(A,is0,is1);
94: MatDestroy(S);
95: MatDestroy(Sexplicit);
97: /* This time just the preconditioner */
98: Create(PETSC_COMM_WORLD,&A,&is0,&is1);
99: MatGetSchurComplement(A,is0,is0,is1,is1,MAT_IGNORE_MATRIX,PETSC_NULL,MAT_INITIAL_MATRIX,&S);
100: PetscPrintf(PETSC_COMM_WORLD,"\nPreconditioning Schur complement of (0,0) in (1,1)\n");
101: MatView(S,PETSC_VIEWER_STDOUT_WORLD);
102: Destroy(A,is0,is1);
103: MatDestroy(S);
105: PetscFinalize();
106: return 0;
107: }