Actual source code: ex21.c

petsc-3.5.4 2015-05-23
Report Typos and Errors
  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: }