1 static char help[] = "Tests MatCreateSubMatrices() for SBAIJ matrices\n\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 Mat BAIJ, SBAIJ, *subBAIJ, *subSBAIJ;
8 PetscViewer viewer;
9 char file[PETSC_MAX_PATH_LEN];
10 PetscBool flg;
11 PetscInt n = 2, issize, M, N;
12 PetscMPIInt rank;
13 IS isrow, iscol, irow[n], icol[n];
14
15 PetscFunctionBeginUser;
16 PetscCall(PetscInitialize(&argc, &args, NULL, 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