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