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,NULL,help);if (ierr) return ierr;
 15:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 16:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 18:   {
 19:     PetscInt ix0a[1],ix0b[1],ix0[2],ix1[1];

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

 37:   usenest = PETSC_FALSE;
 38:   PetscOptionsGetBool(NULL,NULL,"-nest",&usenest,NULL);
 39:   if (usenest) {
 40:     ISLocalToGlobalMapping l2g;
 41:     PetscInt               l2gind[3];
 42:     Mat                    B[9];

 44:     l2gind[0] = (rank-1+size)%size; l2gind[1] = rank; l2gind[2] = (rank+1)%size;
 45:     ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,3,l2gind,PETSC_COPY_VALUES,&l2g);
 46:     for (i=0; i<9; i++) {
 47:       MatCreateAIJ(PETSC_COMM_WORLD,1,1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,NULL,PETSC_DECIDE,NULL,&B[i]);
 48:       MatSetUp(B[i]);
 49:       MatSetLocalToGlobalMapping(B[i],l2g,l2g);
 50:     }
 51:     {
 52:       IS  isx[2];
 53:       Mat Bx00[4],Bx01[2],Bx10[2];
 54:       Mat B00,B01,B10;

 56:       isx[0]  = is0a; isx[1] = is0b;
 57:       Bx00[0] = B[0]; Bx00[1] = B[1]; Bx00[2] = B[3]; Bx00[3] = B[4];
 58:       Bx01[0] = B[2]; Bx01[1] = B[5];
 59:       Bx10[0] = B[6]; Bx10[1] = B[7];

 61:       MatCreateNest(PETSC_COMM_WORLD,2,isx,2,isx,Bx00,&B00);
 62:       MatSetUp(B00);
 63:       MatCreateNest(PETSC_COMM_WORLD,2,isx,1,NULL,Bx01,&B01);
 64:       MatSetUp(B01);
 65:       MatCreateNest(PETSC_COMM_WORLD,1,NULL,2,isx,Bx10,&B10);
 66:       MatSetUp(B10);
 67:       {
 68:         Mat By[4];
 69:         IS  isy[2];

 71:         By[0]  = B00; By[1] = B01; By[2] = B10; By[3] = B[8];
 72:         isy[0] = is0; isy[1] = is1;

 74:         MatCreateNest(PETSC_COMM_WORLD,2,isy,2,isy,By,&A);
 75:         MatSetUp(A);
 76:       }
 77:       MatDestroy(&B00);
 78:       MatDestroy(&B01);
 79:       MatDestroy(&B10);
 80:     }
 81:     for (i=0; i<9; i++) {MatDestroy(&B[i]);}
 82:     ISLocalToGlobalMappingDestroy(&l2g);
 83:   } else {
 84:     ISLocalToGlobalMapping l2g;
 85:     PetscInt               l2gind[9];
 86:     for (i=0; i<3; i++) for (j=0; j<3; j++) l2gind[3*i+j] = ((rank-1+j+size) % size)*3 + i;
 87:     ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,9,l2gind,PETSC_COPY_VALUES,&l2g);
 88:     MatCreateAIJ(PETSC_COMM_WORLD,3,3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,NULL,PETSC_DECIDE,NULL,&A);
 89:     MatSetLocalToGlobalMapping(A,l2g,l2g);
 90:     ISLocalToGlobalMappingDestroy(&l2g);
 91:   }

 93:   {
 94:     Mat A00,A11,A0a0a,A0a0b;
 95:     MatGetLocalSubMatrix(A,isl0,isl0,&A00);
 96:     MatGetLocalSubMatrix(A,isl1,isl1,&A11);
 97:     MatGetLocalSubMatrix(A00,isl0a,isl0a,&A0a0a);
 98:     MatGetLocalSubMatrix(A00,isl0a,isl0b,&A0a0b);

100:     MatSetValueLocal(A0a0a,0,0,100*rank+1,ADD_VALUES);
101:     MatSetValueLocal(A0a0a,0,1,100*rank+2,ADD_VALUES);
102:     MatSetValueLocal(A0a0a,2,2,100*rank+9,ADD_VALUES);

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

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

109:     MatRestoreLocalSubMatrix(A00,isl0a,isl0a,&A0a0a);
110:     MatRestoreLocalSubMatrix(A00,isl0a,isl0b,&A0a0b);
111:     MatRestoreLocalSubMatrix(A,isl0,isl0,&A00);
112:     MatRestoreLocalSubMatrix(A,isl1,isl1,&A11);
113:   }
114:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
115:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

117:   MatComputeExplicitOperator(A,&Aexplicit);
118:   MatView(Aexplicit,PETSC_VIEWER_STDOUT_WORLD);

120:   MatDestroy(&A);
121:   MatDestroy(&Aexplicit);
122:   ISDestroy(&is0a);
123:   ISDestroy(&is0b);
124:   ISDestroy(&is0);
125:   ISDestroy(&is1);
126:   ISDestroy(&isl0a);
127:   ISDestroy(&isl0b);
128:   ISDestroy(&isl0);
129:   ISDestroy(&isl1);
130:   PetscFinalize();
131:   return ierr;
132: }