1 static const char help[] = "Test MatGetLocalSubMatrix() with multiple levels of nesting.\n\n"; 2 3 #include <petscmat.h> 4 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; 12 13 CHKERRQ(PetscInitialize(&argc,&argv,NULL,help)); 14 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 15 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 16 17 { 18 PetscInt ix0a[1],ix0b[1],ix0[2],ix1[1]; 19 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 CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD,1,ix0a,PETSC_COPY_VALUES,&is0a)); 25 CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD,1,ix0b,PETSC_COPY_VALUES,&is0b)); 26 CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD,2,ix0,PETSC_COPY_VALUES,&is0)); 27 CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD,1,ix1,PETSC_COPY_VALUES,&is1)); 28 } 29 { 30 CHKERRQ(ISCreateStride(PETSC_COMM_SELF,6,0,1,&isl0)); 31 CHKERRQ(ISCreateStride(PETSC_COMM_SELF,3,0,1,&isl0a)); 32 CHKERRQ(ISCreateStride(PETSC_COMM_SELF,3,3,1,&isl0b)); 33 CHKERRQ(ISCreateStride(PETSC_COMM_SELF,3,6,1,&isl1)); 34 } 35 36 usenest = PETSC_FALSE; 37 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-nest",&usenest,NULL)); 38 if (usenest) { 39 ISLocalToGlobalMapping l2g; 40 PetscInt l2gind[3]; 41 Mat B[9]; 42 43 l2gind[0] = (rank-1+size)%size; l2gind[1] = rank; l2gind[2] = (rank+1)%size; 44 CHKERRQ(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,3,l2gind,PETSC_COPY_VALUES,&l2g)); 45 for (i=0; i<9; i++) { 46 CHKERRQ(MatCreateAIJ(PETSC_COMM_WORLD,1,1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,NULL,PETSC_DECIDE,NULL,&B[i])); 47 CHKERRQ(MatSetUp(B[i])); 48 CHKERRQ(MatSetLocalToGlobalMapping(B[i],l2g,l2g)); 49 } 50 { 51 IS isx[2]; 52 Mat Bx00[4],Bx01[2],Bx10[2]; 53 Mat B00,B01,B10; 54 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]; 59 60 CHKERRQ(MatCreateNest(PETSC_COMM_WORLD,2,isx,2,isx,Bx00,&B00)); 61 CHKERRQ(MatSetUp(B00)); 62 CHKERRQ(MatCreateNest(PETSC_COMM_WORLD,2,isx,1,NULL,Bx01,&B01)); 63 CHKERRQ(MatSetUp(B01)); 64 CHKERRQ(MatCreateNest(PETSC_COMM_WORLD,1,NULL,2,isx,Bx10,&B10)); 65 CHKERRQ(MatSetUp(B10)); 66 { 67 Mat By[4]; 68 IS isy[2]; 69 70 By[0] = B00; By[1] = B01; By[2] = B10; By[3] = B[8]; 71 isy[0] = is0; isy[1] = is1; 72 73 CHKERRQ(MatCreateNest(PETSC_COMM_WORLD,2,isy,2,isy,By,&A)); 74 CHKERRQ(MatSetUp(A)); 75 } 76 CHKERRQ(MatDestroy(&B00)); 77 CHKERRQ(MatDestroy(&B01)); 78 CHKERRQ(MatDestroy(&B10)); 79 } 80 for (i=0; i<9; i++) CHKERRQ(MatDestroy(&B[i])); 81 CHKERRQ(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 CHKERRQ(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,9,l2gind,PETSC_COPY_VALUES,&l2g)); 87 CHKERRQ(MatCreateAIJ(PETSC_COMM_WORLD,3,3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,NULL,PETSC_DECIDE,NULL,&A)); 88 CHKERRQ(MatSetLocalToGlobalMapping(A,l2g,l2g)); 89 CHKERRQ(ISLocalToGlobalMappingDestroy(&l2g)); 90 } 91 92 { 93 Mat A00,A11,A0a0a,A0a0b; 94 CHKERRQ(MatGetLocalSubMatrix(A,isl0,isl0,&A00)); 95 CHKERRQ(MatGetLocalSubMatrix(A,isl1,isl1,&A11)); 96 CHKERRQ(MatGetLocalSubMatrix(A00,isl0a,isl0a,&A0a0a)); 97 CHKERRQ(MatGetLocalSubMatrix(A00,isl0a,isl0b,&A0a0b)); 98 99 CHKERRQ(MatSetValueLocal(A0a0a,0,0,100*rank+1,ADD_VALUES)); 100 CHKERRQ(MatSetValueLocal(A0a0a,0,1,100*rank+2,ADD_VALUES)); 101 CHKERRQ(MatSetValueLocal(A0a0a,2,2,100*rank+9,ADD_VALUES)); 102 103 CHKERRQ(MatSetValueLocal(A0a0b,1,1,100*rank+50+5,ADD_VALUES)); 104 105 CHKERRQ(MatSetValueLocal(A11,0,0,1000*(rank+1)+1,ADD_VALUES)); 106 CHKERRQ(MatSetValueLocal(A11,1,2,1000*(rank+1)+6,ADD_VALUES)); 107 108 CHKERRQ(MatRestoreLocalSubMatrix(A00,isl0a,isl0a,&A0a0a)); 109 CHKERRQ(MatRestoreLocalSubMatrix(A00,isl0a,isl0b,&A0a0b)); 110 CHKERRQ(MatRestoreLocalSubMatrix(A,isl0,isl0,&A00)); 111 CHKERRQ(MatRestoreLocalSubMatrix(A,isl1,isl1,&A11)); 112 } 113 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 114 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 115 116 CHKERRQ(MatComputeOperator(A,MATAIJ,&Aexplicit)); 117 CHKERRQ(MatView(Aexplicit,PETSC_VIEWER_STDOUT_WORLD)); 118 119 CHKERRQ(MatDestroy(&A)); 120 CHKERRQ(MatDestroy(&Aexplicit)); 121 CHKERRQ(ISDestroy(&is0a)); 122 CHKERRQ(ISDestroy(&is0b)); 123 CHKERRQ(ISDestroy(&is0)); 124 CHKERRQ(ISDestroy(&is1)); 125 CHKERRQ(ISDestroy(&isl0a)); 126 CHKERRQ(ISDestroy(&isl0b)); 127 CHKERRQ(ISDestroy(&isl0)); 128 CHKERRQ(ISDestroy(&isl1)); 129 CHKERRQ(PetscFinalize()); 130 return 0; 131 } 132 133 /*TEST 134 135 test: 136 nsize: 3 137 138 test: 139 suffix: nest 140 nsize: 3 141 args: -nest 142 143 TEST*/ 144