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: }