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