Actual source code: ex21.c

  1: static const char help[] = "Tests MatGetSchurComplement\n";

  3: #include <petscksp.h>

  5: PetscErrorCode Create(MPI_Comm comm,Mat *inA,IS *is0,IS *is1)
  6: {
  7:   Mat            A;
  8:   PetscInt       r,rend,M;
  9:   PetscMPIInt    rank;

 12:   *inA = 0;
 13:   MatCreate(comm,&A);
 14:   MatSetSizes(A,4,4,PETSC_DETERMINE,PETSC_DETERMINE);
 15:   MatSetFromOptions(A);
 16:   MatSetUp(A);
 17:   MatGetOwnershipRange(A,&r,&rend);
 18:   MatGetSize(A,&M,NULL);

 20:   ISCreateStride(comm,2,r,1,is0);
 21:   ISCreateStride(comm,2,r+2,1,is1);

 23:   MPI_Comm_rank(comm,&rank);

 25:   {
 26:     PetscInt    rows[4],cols0[5],cols1[5],cols2[3],cols3[3];
 27:     PetscScalar RR = 1000.*rank,vals0[5],vals1[4],vals2[3],vals3[3];

 29:     rows[0]            = r;
 30:     rows[1]            = r+1;
 31:     rows[2]            = r+2;
 32:     rows[3]            = r+3;

 34:     cols0[0]           = r+0;
 35:     cols0[1]           = r+1;
 36:     cols0[2]           = r+3;
 37:     cols0[3]           = (r+4)%M;
 38:     cols0[4]           = (r+M-4)%M;

 40:     cols1[0]           = r+1;
 41:     cols1[1]           = r+2;
 42:     cols1[2]           = (r+4+1)%M;
 43:     cols1[3]           = (r+M-4+1)%M;

 45:     cols2[0]           = r;
 46:     cols2[1]           = r+2;
 47:     cols2[2]           = (r+4+2)%M;

 49:     cols3[0]           = r+1;
 50:     cols3[1]           = r+3;
 51:     cols3[2]           = (r+4+3)%M;

 53:     vals0[0] = RR+1.;
 54:     vals0[1] = RR+2.;
 55:     vals0[2] = RR+3.;
 56:     vals0[3] = RR+4.;
 57:     vals0[4] = RR+5.;

 59:     vals1[0] = RR+6.;
 60:     vals1[1] = RR+7.;
 61:     vals1[2] = RR+8.;
 62:     vals1[3] = RR+9.;

 64:     vals2[0] = RR+10.;
 65:     vals2[1] = RR+11.;
 66:     vals2[2] = RR+12.;

 68:     vals3[0] = RR+13.;
 69:     vals3[1] = RR+14.;
 70:     vals3[2] = RR+15.;
 71:     MatSetValues(A,1,&rows[0],5,cols0,vals0,INSERT_VALUES);
 72:     MatSetValues(A,1,&rows[1],4,cols1,vals1,INSERT_VALUES);
 73:     MatSetValues(A,1,&rows[2],3,cols2,vals2,INSERT_VALUES);
 74:     MatSetValues(A,1,&rows[3],3,cols3,vals3,INSERT_VALUES);
 75:   }
 76:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 77:   MatAssemblyEnd  (A,MAT_FINAL_ASSEMBLY);
 78:   *inA = A;
 79:   return 0;
 80: }

 82: PetscErrorCode Destroy(Mat *A,IS *is0,IS *is1)
 83: {
 85:   MatDestroy(A);
 86:   ISDestroy(is0);
 87:   ISDestroy(is1);
 88:   return 0;
 89: }

 91: int main(int argc,char *argv[])
 92: {
 94:   Mat                        A,S = NULL,Sexplicit = NULL;
 95:   MatSchurComplementAinvType ainv_type = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
 96:   IS                         is0,is1;

 98:   PetscInitialize(&argc,&argv,0,help);
 99:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex21","KSP");
100:   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);
101:   PetscOptionsEnd();

103:   /* Test the Schur complement one way */
104:   Create(PETSC_COMM_WORLD,&A,&is0,&is1);
105:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);
106:   ISView(is0,PETSC_VIEWER_STDOUT_WORLD);
107:   ISView(is1,PETSC_VIEWER_STDOUT_WORLD);
108:   MatGetSchurComplement(A,is0,is0,is1,is1,MAT_INITIAL_MATRIX,&S,ainv_type,MAT_IGNORE_MATRIX,NULL);
109:   MatSetFromOptions(S);
110:   MatComputeOperator(S,MATAIJ,&Sexplicit);
111:   PetscPrintf(PETSC_COMM_WORLD,"\nExplicit Schur complement of (0,0) in (1,1)\n");
112:   MatView(Sexplicit,PETSC_VIEWER_STDOUT_WORLD);
113:   Destroy(&A,&is0,&is1);
114:   MatDestroy(&S);
115:   MatDestroy(&Sexplicit);

117:   /* And the other */
118:   Create(PETSC_COMM_WORLD,&A,&is0,&is1);
119:   MatGetSchurComplement(A,is1,is1,is0,is0,MAT_INITIAL_MATRIX,&S,ainv_type,MAT_IGNORE_MATRIX,NULL);
120:   MatSetFromOptions(S);
121:   MatComputeOperator(S,MATAIJ,&Sexplicit);
122:   PetscPrintf(PETSC_COMM_WORLD,"\nExplicit Schur complement of (1,1) in (0,0)\n");
123:   MatView(Sexplicit,PETSC_VIEWER_STDOUT_WORLD);
124:   Destroy(&A,&is0,&is1);
125:   MatDestroy(&S);
126:   MatDestroy(&Sexplicit);

128:   /* This time just the preconditioning matrix. */
129:   Create(PETSC_COMM_WORLD,&A,&is0,&is1);
130:   MatGetSchurComplement(A,is0,is0,is1,is1,MAT_IGNORE_MATRIX,NULL,ainv_type,MAT_INITIAL_MATRIX,&S);
131:   MatSetFromOptions(S);
132:   PetscPrintf(PETSC_COMM_WORLD,"\nPreconditioning Schur complement of (0,0) in (1,1)\n");
133:   MatView(S,PETSC_VIEWER_STDOUT_WORLD);
134:   /* Modify and refresh */
135:   MatShift(A,1.);
136:   MatGetSchurComplement(A,is0,is0,is1,is1,MAT_IGNORE_MATRIX,NULL,ainv_type,MAT_REUSE_MATRIX,&S);
137:   PetscPrintf(PETSC_COMM_WORLD,"\nAfter update\n");
138:   MatView(S,PETSC_VIEWER_STDOUT_WORLD);
139:   Destroy(&A,&is0,&is1);
140:   MatDestroy(&S);

142:   PetscFinalize();
143:   return 0;
144: }

146: /*TEST
147:   test:
148:     suffix: diag_1
149:     args: -mat_schur_complement_ainv_type diag
150:     nsize: 1
151:   test:
152:     suffix: blockdiag_1
153:     args: -mat_schur_complement_ainv_type blockdiag
154:     nsize: 1
155:   test:
156:     suffix: diag_2
157:     args: -mat_schur_complement_ainv_type diag
158:     nsize: 2
159:   test:
160:     suffix: blockdiag_2
161:     args: -mat_schur_complement_ainv_type blockdiag
162:     nsize: 2
163:   test:
164:     # does not work with single because residual norm computed by GMRES recurrence formula becomes invalid
165:     requires: !single
166:     suffix: diag_3
167:     args: -mat_schur_complement_ainv_type diag
168:     nsize: 3
169:   test:
170:     # does not work with single because residual norm computed by GMRES recurrence formula becomes invalid
171:     requires: !single
172:     suffix: blockdiag_3
173:     args: -mat_schur_complement_ainv_type blockdiag
174:     nsize: 3
175: TEST*/