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