Actual source code: ex159.c

petsc-3.3-p6 2013-02-11
  1: static const char help[] = "Test MatGetLocalSubMatrix() with multiple levels of nesting.\n\n";

  3: #include <petscmat.h>

  5: int main(int argc, char *argv[])
  6: {
  8:   IS is0a,is0b,is0,is1,isl0a,isl0b,isl0,isl1;
  9:   Mat A,Aexplicit;
 10:   PetscBool usenest;
 11:   PetscMPIInt rank,size;
 12:   PetscInt i,j;

 14:   PetscInitialize(&argc,&argv,PETSC_NULL,help);
 15:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 16:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 18:   {
 19:     const PetscInt ix0a[] = {rank*2+0},ix0b[] = {rank*2+1},ix0[] = {rank*3+0,rank*3+1},ix1[] = {rank*3+2};
 20:     ISCreateGeneral(PETSC_COMM_WORLD,1,ix0a,PETSC_COPY_VALUES,&is0a);
 21:     ISCreateGeneral(PETSC_COMM_WORLD,1,ix0b,PETSC_COPY_VALUES,&is0b);
 22:     ISCreateGeneral(PETSC_COMM_WORLD,2,ix0,PETSC_COPY_VALUES,&is0);
 23:     ISCreateGeneral(PETSC_COMM_WORLD,1,ix1,PETSC_COPY_VALUES,&is1);
 24:   }
 25:   {
 26:     ISCreateStride(PETSC_COMM_SELF,6,0,1,&isl0);
 27:     ISCreateStride(PETSC_COMM_SELF,3,0,1,&isl0a);
 28:     ISCreateStride(PETSC_COMM_SELF,3,3,1,&isl0b);
 29:     ISCreateStride(PETSC_COMM_SELF,3,6,1,&isl1);
 30:   }

 32:   usenest = PETSC_FALSE;
 33:   PetscOptionsGetBool(PETSC_NULL,"-nest",&usenest,PETSC_NULL);
 34:   if (usenest) {
 35:     ISLocalToGlobalMapping l2g;
 36:     const PetscInt l2gind[3] = {(rank-1+size)%size,rank,(rank+1)%size};
 37:     Mat B[9];
 38:     ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,3,l2gind,PETSC_COPY_VALUES,&l2g);
 39:     for (i=0; i<9; i++) {
 40:       MatCreateAIJ(PETSC_COMM_WORLD,1,1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_NULL,PETSC_DECIDE,PETSC_NULL,&B[i]);
 41:       MatSetUp(B[i]);
 42:       MatSetLocalToGlobalMapping(B[i],l2g,l2g);
 43:     }
 44:     {
 45:       const IS isx[] = {is0a,is0b};
 46:       const Mat Bx00[] = {B[0],B[1],B[3],B[4]},Bx01[] = {B[2],B[5]},Bx10[] = {B[6],B[7]};
 47:       Mat B00,B01,B10;
 48:       MatCreateNest(PETSC_COMM_WORLD,2,isx,2,isx,Bx00,&B00);
 49:       MatSetUp(B00);
 50:       MatCreateNest(PETSC_COMM_WORLD,2,isx,1,PETSC_NULL,Bx01,&B01);
 51:       MatSetUp(B01);
 52:       MatCreateNest(PETSC_COMM_WORLD,1,PETSC_NULL,2,isx,Bx10,&B10);
 53:       MatSetUp(B10);
 54:       {
 55:         Mat By[] = {B00,B01,B10,B[8]};
 56:         IS isy[] = {is0,is1};
 57:         MatCreateNest(PETSC_COMM_WORLD,2,isy,2,isy,By,&A);
 58:         MatSetUp(A);
 59:       }
 60:       MatDestroy(&B00);
 61:       MatDestroy(&B01);
 62:       MatDestroy(&B10);
 63:     }
 64:     for (i=0; i<9; i++) {MatDestroy(&B[i]);}
 65:     ISLocalToGlobalMappingDestroy(&l2g);
 66:   } else {
 67:     ISLocalToGlobalMapping l2g;
 68:     PetscInt l2gind[9];
 69:     for (i=0; i<3; i++) for (j=0; j<3; j++) l2gind[3*i+j] = ((rank-1+j+size) % size)*3 + i;
 70:     ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,9,l2gind,PETSC_COPY_VALUES,&l2g);
 71:     MatCreateAIJ(PETSC_COMM_WORLD,3,3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_NULL,PETSC_DECIDE,PETSC_NULL,&A);
 72:     MatSetLocalToGlobalMapping(A,l2g,l2g);
 73:     ISLocalToGlobalMappingDestroy(&l2g);
 74:   }

 76:   {
 77:     Mat A00,A11,A0a0a,A0a0b;
 78:     MatGetLocalSubMatrix(A,isl0,isl0,&A00);
 79:     MatGetLocalSubMatrix(A,isl1,isl1,&A11);
 80:     MatGetLocalSubMatrix(A00,isl0a,isl0a,&A0a0a);
 81:     MatGetLocalSubMatrix(A00,isl0a,isl0b,&A0a0b);

 83:     MatSetValueLocal(A0a0a,0,0,100*rank+1,ADD_VALUES);
 84:     MatSetValueLocal(A0a0a,0,1,100*rank+2,ADD_VALUES);
 85:     MatSetValueLocal(A0a0a,2,2,100*rank+9,ADD_VALUES);

 87:     MatSetValueLocal(A0a0b,1,1,100*rank+50+5,ADD_VALUES);

 89:     MatSetValueLocal(A11,0,0,1000*(rank+1)+1,ADD_VALUES);
 90:     MatSetValueLocal(A11,1,2,1000*(rank+1)+6,ADD_VALUES);

 92:     MatRestoreLocalSubMatrix(A00,isl0a,isl0a,&A0a0a);
 93:     MatRestoreLocalSubMatrix(A00,isl0a,isl0b,&A0a0b);
 94:     MatRestoreLocalSubMatrix(A,isl0,isl0,&A00);
 95:     MatRestoreLocalSubMatrix(A,isl1,isl1,&A11);
 96:   }
 97:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 98:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

100:   MatComputeExplicitOperator(A,&Aexplicit);
101:   MatView(Aexplicit,PETSC_VIEWER_STDOUT_WORLD);

103:   MatDestroy(&A);
104:   MatDestroy(&Aexplicit);
105:   ISDestroy(&is0a);
106:   ISDestroy(&is0b);
107:   ISDestroy(&is0);
108:   ISDestroy(&is1);
109:   ISDestroy(&isl0a);
110:   ISDestroy(&isl0b);
111:   ISDestroy(&isl0);
112:   ISDestroy(&isl1);
113:   PetscFinalize();
114:   return 0;
115: }