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: }