1 static char help[] = "Test MatCreateSubMatrix with -mat_type nest and block sizes.\n"; 2 3 #include <petscmat.h> 4 5 int main(int argc, char **argv) 6 { 7 Mat A, B, C, mats[6]; 8 IS rows[2]; 9 ISLocalToGlobalMapping cmap, rmap; 10 const PetscInt indices[3] = {0, 1, 2}; 11 PetscInt i; 12 PetscMPIInt size; 13 14 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 15 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 16 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only coded for one process"); 17 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 18 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2,1)); 19 PetscCall(MatSetBlockSizes(A,2,1)); 20 PetscCall(MatSetType(A,MATAIJ)); 21 PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,2,1,indices,PETSC_COPY_VALUES,&rmap)); 22 PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,1,indices,PETSC_COPY_VALUES,&cmap)); 23 PetscCall(MatSetLocalToGlobalMapping(A,rmap,cmap)); 24 PetscCall(ISLocalToGlobalMappingDestroy(&rmap)); 25 PetscCall(ISLocalToGlobalMappingDestroy(&cmap)); 26 PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,2,indices,PETSC_COPY_VALUES,&rmap)); 27 PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,1,3,indices,PETSC_COPY_VALUES,&cmap)); 28 PetscCall(MatCreate(PETSC_COMM_WORLD,&B)); 29 PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,2,3)); 30 PetscCall(MatSetBlockSizes(B,1,1)); 31 PetscCall(MatSetType(B,MATAIJ)); 32 PetscCall(MatSetLocalToGlobalMapping(B,rmap,cmap)); 33 PetscCall(ISLocalToGlobalMappingDestroy(&rmap)); 34 PetscCall(ISLocalToGlobalMappingDestroy(&cmap)); 35 PetscCall(MatSetUp(A)); 36 PetscCall(MatSetUp(B)); 37 mats[0] = A; 38 mats[1] = B; 39 mats[2] = A; 40 mats[3] = NULL; 41 mats[4] = B; 42 mats[5] = A; 43 PetscCall(MatCreateNest(PETSC_COMM_WORLD,2,NULL,3,NULL,mats,&C)); 44 PetscCall(MatSetUp(C)); 45 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 46 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 47 PetscCall(MatView(C,NULL)); 48 PetscCall(MatNestGetISs(C,rows,NULL)); 49 for (i=0; i<2; i++) { 50 Mat submat; 51 IS cols[3]; 52 PetscInt j; 53 PetscCall(ISView(rows[i],NULL)); 54 PetscCall(MatCreateSubMatrix(C,rows[i],NULL,MAT_INITIAL_MATRIX,&submat)); 55 PetscCall(MatView(submat,NULL)); 56 PetscCall(MatNestGetISs(submat,NULL,cols)); 57 for (j=0; j<3; j++) { 58 PetscCall(ISView(cols[j],NULL)); 59 } 60 PetscCall(MatDestroy(&submat)); 61 } 62 PetscCall(MatDestroy(&C)); 63 PetscCall(MatDestroy(&B)); 64 PetscCall(MatDestroy(&A)); 65 PetscCall(PetscFinalize()); 66 return 0; 67 } 68 69 /*TEST 70 71 test: 72 73 TEST*/ 74