1 2 static char help[] = "Tests MatCreateSubMatrices() for SBAIJ matrices\n\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc,char **args) 7 { 8 Mat BAIJ,SBAIJ,*subBAIJ,*subSBAIJ; 9 PetscViewer viewer; 10 char file[PETSC_MAX_PATH_LEN]; 11 PetscBool flg; 12 PetscInt n = 2,issize,M,N; 13 PetscMPIInt rank; 14 IS isrow,iscol,irow[n],icol[n]; 15 16 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 17 PetscCall(PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg)); 18 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&viewer)); 19 20 PetscCall(MatCreate(PETSC_COMM_WORLD,&BAIJ)); 21 PetscCall(MatSetType(BAIJ,MATMPIBAIJ)); 22 PetscCall(MatLoad(BAIJ,viewer)); 23 PetscCall(PetscViewerDestroy(&viewer)); 24 25 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&viewer)); 26 PetscCall(MatCreate(PETSC_COMM_WORLD,&SBAIJ)); 27 PetscCall(MatSetType(SBAIJ,MATMPISBAIJ)); 28 PetscCall(MatLoad(SBAIJ,viewer)); 29 PetscCall(PetscViewerDestroy(&viewer)); 30 31 PetscCall(MatGetSize(BAIJ,&M,&N)); 32 issize = M/4; 33 PetscCall(ISCreateStride(PETSC_COMM_SELF,issize,0,1,&isrow)); 34 irow[0] = irow[1] = isrow; 35 issize = N/2; 36 PetscCall(ISCreateStride(PETSC_COMM_SELF,issize,0,1,&iscol)); 37 icol[0] = icol[1] = iscol; 38 PetscCall(MatCreateSubMatrices(BAIJ,n,irow,icol,MAT_INITIAL_MATRIX,&subBAIJ)); 39 PetscCall(MatCreateSubMatrices(BAIJ,n,irow,icol,MAT_REUSE_MATRIX,&subBAIJ)); 40 41 /* irow and icol must be same for SBAIJ matrices! */ 42 icol[0] = icol[1] = isrow; 43 PetscCall(MatCreateSubMatrices(SBAIJ,n,irow,icol,MAT_INITIAL_MATRIX,&subSBAIJ)); 44 PetscCall(MatCreateSubMatrices(SBAIJ,n,irow,icol,MAT_REUSE_MATRIX,&subSBAIJ)); 45 46 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 47 if (rank == 0) { 48 PetscCall(MatView(subBAIJ[0],PETSC_VIEWER_STDOUT_SELF)); 49 PetscCall(MatView(subSBAIJ[0],PETSC_VIEWER_STDOUT_SELF)); 50 } 51 52 /* Free data structures */ 53 PetscCall(ISDestroy(&isrow)); 54 PetscCall(ISDestroy(&iscol)); 55 PetscCall(MatDestroySubMatrices(n,&subBAIJ)); 56 PetscCall(MatDestroySubMatrices(n,&subSBAIJ)); 57 PetscCall(MatDestroy(&BAIJ)); 58 PetscCall(MatDestroy(&SBAIJ)); 59 60 PetscCall(PetscFinalize()); 61 return 0; 62 } 63