Actual source code: ex159.c

  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: {
  7:   IS          is0a, is0b, is0, is1, isl0a, isl0b, isl0, isl1;
  8:   Mat         A, Aexplicit;
  9:   PetscBool   usenest;
 10:   PetscMPIInt rank, size;
 11:   PetscInt    i, j;

 14:   PetscInitialize(&argc, &argv, NULL, help);
 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;
 24:     ix0[1]  = rank * 3 + 1;
 25:     ix1[0]  = rank * 3 + 2;
 26:     ISCreateGeneral(PETSC_COMM_WORLD, 1, ix0a, PETSC_COPY_VALUES, &is0a);
 27:     ISCreateGeneral(PETSC_COMM_WORLD, 1, ix0b, PETSC_COPY_VALUES, &is0b);
 28:     ISCreateGeneral(PETSC_COMM_WORLD, 2, ix0, PETSC_COPY_VALUES, &is0);
 29:     ISCreateGeneral(PETSC_COMM_WORLD, 1, ix1, PETSC_COPY_VALUES, &is1);
 30:   }
 31:   {
 32:     ISCreateStride(PETSC_COMM_SELF, 6, 0, 1, &isl0);
 33:     ISCreateStride(PETSC_COMM_SELF, 3, 0, 1, &isl0a);
 34:     ISCreateStride(PETSC_COMM_SELF, 3, 3, 1, &isl0b);
 35:     ISCreateStride(PETSC_COMM_SELF, 3, 6, 1, &isl1);
 36:   }

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

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

 59:       isx[0]  = is0a;
 60:       isx[1]  = is0b;
 61:       Bx00[0] = B[0];
 62:       Bx00[1] = B[1];
 63:       Bx00[2] = B[3];
 64:       Bx00[3] = B[4];
 65:       Bx01[0] = B[2];
 66:       Bx01[1] = B[5];
 67:       Bx10[0] = B[6];
 68:       Bx10[1] = B[7];

 70:       MatCreateNest(PETSC_COMM_WORLD, 2, isx, 2, isx, Bx00, &B00);
 71:       MatSetUp(B00);
 72:       MatCreateNest(PETSC_COMM_WORLD, 2, isx, 1, NULL, Bx01, &B01);
 73:       MatSetUp(B01);
 74:       MatCreateNest(PETSC_COMM_WORLD, 1, NULL, 2, isx, Bx10, &B10);
 75:       MatSetUp(B10);
 76:       {
 77:         Mat By[4];
 78:         IS  isy[2];

 80:         By[0]  = B00;
 81:         By[1]  = B01;
 82:         By[2]  = B10;
 83:         By[3]  = B[8];
 84:         isy[0] = is0;
 85:         isy[1] = is1;

 87:         MatCreateNest(PETSC_COMM_WORLD, 2, isy, 2, isy, By, &A);
 88:         MatSetUp(A);
 89:       }
 90:       MatDestroy(&B00);
 91:       MatDestroy(&B01);
 92:       MatDestroy(&B10);
 93:     }
 94:     for (i = 0; i < 9; i++) MatDestroy(&B[i]);
 95:     ISLocalToGlobalMappingDestroy(&l2g);
 96:   } else {
 97:     ISLocalToGlobalMapping l2g;
 98:     PetscInt               l2gind[9];
 99:     for (i = 0; i < 3; i++)
100:       for (j = 0; j < 3; j++) l2gind[3 * i + j] = ((rank - 1 + j + size) % size) * 3 + i;
101:     ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, 9, l2gind, PETSC_COPY_VALUES, &l2g);
102:     MatCreateAIJ(PETSC_COMM_WORLD, 3, 3, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, NULL, PETSC_DECIDE, NULL, &A);
103:     MatSetLocalToGlobalMapping(A, l2g, l2g);
104:     ISLocalToGlobalMappingDestroy(&l2g);
105:   }

107:   {
108:     Mat A00, A11, A0a0a, A0a0b;
109:     MatGetLocalSubMatrix(A, isl0, isl0, &A00);
110:     MatGetLocalSubMatrix(A, isl1, isl1, &A11);
111:     MatGetLocalSubMatrix(A00, isl0a, isl0a, &A0a0a);
112:     MatGetLocalSubMatrix(A00, isl0a, isl0b, &A0a0b);

114:     MatSetValueLocal(A0a0a, 0, 0, 100 * rank + 1, ADD_VALUES);
115:     MatSetValueLocal(A0a0a, 0, 1, 100 * rank + 2, ADD_VALUES);
116:     MatSetValueLocal(A0a0a, 2, 2, 100 * rank + 9, ADD_VALUES);

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

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

123:     MatRestoreLocalSubMatrix(A00, isl0a, isl0a, &A0a0a);
124:     MatRestoreLocalSubMatrix(A00, isl0a, isl0b, &A0a0b);
125:     MatRestoreLocalSubMatrix(A, isl0, isl0, &A00);
126:     MatRestoreLocalSubMatrix(A, isl1, isl1, &A11);
127:   }
128:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
129:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

131:   MatComputeOperator(A, MATAIJ, &Aexplicit);
132:   MatView(Aexplicit, PETSC_VIEWER_STDOUT_WORLD);

134:   MatDestroy(&A);
135:   MatDestroy(&Aexplicit);
136:   ISDestroy(&is0a);
137:   ISDestroy(&is0b);
138:   ISDestroy(&is0);
139:   ISDestroy(&is1);
140:   ISDestroy(&isl0a);
141:   ISDestroy(&isl0b);
142:   ISDestroy(&isl0);
143:   ISDestroy(&isl1);
144:   PetscFinalize();
145:   return 0;
146: }

148: /*TEST

150:    test:
151:       nsize: 3

153:    test:
154:       suffix: nest
155:       nsize: 3
156:       args: -nest

158: TEST*/