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;
13: PetscInitialize(&argc,&argv,NULL,help);
14: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
15: MPI_Comm_size(PETSC_COMM_WORLD,&size);
17: {
18: PetscInt ix0a[1],ix0b[1],ix0[2],ix1[1];
20: ix0a[0] = rank*2+0;
21: ix0b[0] = rank*2+1;
22: ix0[0] = rank*3+0; ix0[1] = rank*3+1;
23: ix1[0] = rank*3+2;
24: ISCreateGeneral(PETSC_COMM_WORLD,1,ix0a,PETSC_COPY_VALUES,&is0a);
25: ISCreateGeneral(PETSC_COMM_WORLD,1,ix0b,PETSC_COPY_VALUES,&is0b);
26: ISCreateGeneral(PETSC_COMM_WORLD,2,ix0,PETSC_COPY_VALUES,&is0);
27: ISCreateGeneral(PETSC_COMM_WORLD,1,ix1,PETSC_COPY_VALUES,&is1);
28: }
29: {
30: ISCreateStride(PETSC_COMM_SELF,6,0,1,&isl0);
31: ISCreateStride(PETSC_COMM_SELF,3,0,1,&isl0a);
32: ISCreateStride(PETSC_COMM_SELF,3,3,1,&isl0b);
33: ISCreateStride(PETSC_COMM_SELF,3,6,1,&isl1);
34: }
36: usenest = PETSC_FALSE;
37: PetscOptionsGetBool(NULL,NULL,"-nest",&usenest,NULL);
38: if (usenest) {
39: ISLocalToGlobalMapping l2g;
40: PetscInt l2gind[3];
41: Mat B[9];
43: l2gind[0] = (rank-1+size)%size; l2gind[1] = rank; l2gind[2] = (rank+1)%size;
44: ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,3,l2gind,PETSC_COPY_VALUES,&l2g);
45: for (i=0; i<9; i++) {
46: MatCreateAIJ(PETSC_COMM_WORLD,1,1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,NULL,PETSC_DECIDE,NULL,&B[i]);
47: MatSetUp(B[i]);
48: MatSetLocalToGlobalMapping(B[i],l2g,l2g);
49: }
50: {
51: IS isx[2];
52: Mat Bx00[4],Bx01[2],Bx10[2];
53: Mat B00,B01,B10;
55: isx[0] = is0a; isx[1] = is0b;
56: Bx00[0] = B[0]; Bx00[1] = B[1]; Bx00[2] = B[3]; Bx00[3] = B[4];
57: Bx01[0] = B[2]; Bx01[1] = B[5];
58: Bx10[0] = B[6]; Bx10[1] = B[7];
60: MatCreateNest(PETSC_COMM_WORLD,2,isx,2,isx,Bx00,&B00);
61: MatSetUp(B00);
62: MatCreateNest(PETSC_COMM_WORLD,2,isx,1,NULL,Bx01,&B01);
63: MatSetUp(B01);
64: MatCreateNest(PETSC_COMM_WORLD,1,NULL,2,isx,Bx10,&B10);
65: MatSetUp(B10);
66: {
67: Mat By[4];
68: IS isy[2];
70: By[0] = B00; By[1] = B01; By[2] = B10; By[3] = B[8];
71: isy[0] = is0; isy[1] = is1;
73: MatCreateNest(PETSC_COMM_WORLD,2,isy,2,isy,By,&A);
74: MatSetUp(A);
75: }
76: MatDestroy(&B00);
77: MatDestroy(&B01);
78: MatDestroy(&B10);
79: }
80: for (i=0; i<9; i++) MatDestroy(&B[i]);
81: ISLocalToGlobalMappingDestroy(&l2g);
82: } else {
83: ISLocalToGlobalMapping l2g;
84: PetscInt l2gind[9];
85: for (i=0; i<3; i++) for (j=0; j<3; j++) l2gind[3*i+j] = ((rank-1+j+size) % size)*3 + i;
86: ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,9,l2gind,PETSC_COPY_VALUES,&l2g);
87: MatCreateAIJ(PETSC_COMM_WORLD,3,3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,NULL,PETSC_DECIDE,NULL,&A);
88: MatSetLocalToGlobalMapping(A,l2g,l2g);
89: ISLocalToGlobalMappingDestroy(&l2g);
90: }
92: {
93: Mat A00,A11,A0a0a,A0a0b;
94: MatGetLocalSubMatrix(A,isl0,isl0,&A00);
95: MatGetLocalSubMatrix(A,isl1,isl1,&A11);
96: MatGetLocalSubMatrix(A00,isl0a,isl0a,&A0a0a);
97: MatGetLocalSubMatrix(A00,isl0a,isl0b,&A0a0b);
99: MatSetValueLocal(A0a0a,0,0,100*rank+1,ADD_VALUES);
100: MatSetValueLocal(A0a0a,0,1,100*rank+2,ADD_VALUES);
101: MatSetValueLocal(A0a0a,2,2,100*rank+9,ADD_VALUES);
103: MatSetValueLocal(A0a0b,1,1,100*rank+50+5,ADD_VALUES);
105: MatSetValueLocal(A11,0,0,1000*(rank+1)+1,ADD_VALUES);
106: MatSetValueLocal(A11,1,2,1000*(rank+1)+6,ADD_VALUES);
108: MatRestoreLocalSubMatrix(A00,isl0a,isl0a,&A0a0a);
109: MatRestoreLocalSubMatrix(A00,isl0a,isl0b,&A0a0b);
110: MatRestoreLocalSubMatrix(A,isl0,isl0,&A00);
111: MatRestoreLocalSubMatrix(A,isl1,isl1,&A11);
112: }
113: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
114: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
116: MatComputeOperator(A,MATAIJ,&Aexplicit);
117: MatView(Aexplicit,PETSC_VIEWER_STDOUT_WORLD);
119: MatDestroy(&A);
120: MatDestroy(&Aexplicit);
121: ISDestroy(&is0a);
122: ISDestroy(&is0b);
123: ISDestroy(&is0);
124: ISDestroy(&is1);
125: ISDestroy(&isl0a);
126: ISDestroy(&isl0b);
127: ISDestroy(&isl0);
128: ISDestroy(&isl1);
129: PetscFinalize();
130: return 0;
131: }
133: /*TEST
135: test:
136: nsize: 3
138: test:
139: suffix: nest
140: nsize: 3
141: args: -nest
143: TEST*/