Actual source code: ex21.c
petsc-3.5.4 2015-05-23
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: MatSetUp(A);
21: MatGetOwnershipRange(A,&r,&rend);
22: MatGetSize(A,&M,NULL);
24: ISCreateStride(comm,2,r,1,is0);
25: ISCreateStride(comm,2,r+2,1,is1);
27: MPI_Comm_rank(comm,&rank);
29: {
30: PetscInt
31: rows[] = {r,r+1,r+2,r+3},
32: cols0[] = {r+0,r+1,r+3,(r+4)%M,(r+M-4)%M},
33: cols1[] = {r+1,r+2,(r+4+1)%M,(r+M-4+1)%M},
34: cols2[] = {r,r+2,(r+4+2)%M},
35: cols3[] = {r+1,r+3,(r+4+3)%M};
36: PetscScalar RR = 1000*rank,
37: vals0[] = {RR+1,RR+2,RR+3,RR+4,RR+5},
38: vals1[] = {RR+6,RR+7,RR+8,RR+9},
39: vals2[] = {RR+10,RR+11,RR+12},
40: vals3[] = {RR+13,RR+14,RR+15};
41: MatSetValues(A,1,&rows[0],5,cols0,vals0,INSERT_VALUES);
42: MatSetValues(A,1,&rows[1],4,cols1,vals1,INSERT_VALUES);
43: MatSetValues(A,1,&rows[2],3,cols2,vals2,INSERT_VALUES);
44: MatSetValues(A,1,&rows[3],3,cols3,vals3,INSERT_VALUES);
45: }
46: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
47: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
48: *inA = A;
49: return(0);
50: }
54: PetscErrorCode Destroy(Mat *A,IS *is0,IS *is1)
55: {
59: MatDestroy(A);
60: ISDestroy(is0);
61: ISDestroy(is1);
62: return(0);
63: }
67: int main(int argc,char *argv[])
68: {
70: Mat A,S,Sexplicit;
71: IS is0,is1;
73: PetscInitialize(&argc,&argv,0,help);
75: /* Test the Schur complement one way */
76: Create(PETSC_COMM_WORLD,&A,&is0,&is1);
77: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
78: ISView(is0,PETSC_VIEWER_STDOUT_WORLD);
79: ISView(is1,PETSC_VIEWER_STDOUT_WORLD);
80: MatGetSchurComplement(A,is0,is0,is1,is1,MAT_INITIAL_MATRIX,&S,PETSC_FALSE,MAT_IGNORE_MATRIX,NULL);
81: MatComputeExplicitOperator(S,&Sexplicit);
82: PetscPrintf(PETSC_COMM_WORLD,"\nExplicit Schur complement of (0,0) in (1,1)\n");
83: MatView(Sexplicit,PETSC_VIEWER_STDOUT_WORLD);
84: Destroy(&A,&is0,&is1);
85: MatDestroy(&S);
86: MatDestroy(&Sexplicit);
88: /* And the other */
89: Create(PETSC_COMM_WORLD,&A,&is0,&is1);
90: MatGetSchurComplement(A,is1,is1,is0,is0,MAT_INITIAL_MATRIX,&S,PETSC_FALSE,MAT_IGNORE_MATRIX,NULL);
91: MatComputeExplicitOperator(S,&Sexplicit);
92: PetscPrintf(PETSC_COMM_WORLD,"\nExplicit Schur complement of (1,1) in (0,0)\n");
93: MatView(Sexplicit,PETSC_VIEWER_STDOUT_WORLD);
94: Destroy(&A,&is0,&is1);
95: MatDestroy(&S);
96: MatDestroy(&Sexplicit);
98: /* This time just the preconditioner */
99: Create(PETSC_COMM_WORLD,&A,&is0,&is1);
100: MatGetSchurComplement(A,is0,is0,is1,is1,MAT_IGNORE_MATRIX,NULL,PETSC_FALSE,MAT_INITIAL_MATRIX,&S);
101: PetscPrintf(PETSC_COMM_WORLD,"\nPreconditioning Schur complement of (0,0) in (1,1)\n");
102: MatView(S,PETSC_VIEWER_STDOUT_WORLD);
103: /* Modify and refresh */
104: MatShift(A,1.);
105: MatGetSchurComplement(A,is0,is0,is1,is1,MAT_IGNORE_MATRIX,NULL,PETSC_FALSE,MAT_REUSE_MATRIX,&S);
106: PetscPrintf(PETSC_COMM_WORLD,"\nAfter update\n");
107: MatView(S,PETSC_VIEWER_STDOUT_WORLD);
108: Destroy(&A,&is0,&is1);
109: MatDestroy(&S);
111: PetscFinalize();
112: return 0;
113: }