xref: /petsc/src/mat/tests/ex249.c (revision 58d68138c660dfb4e9f5b03334792cd4f2ffd7cc)
1 static char help[] = "Test MatCreateSubMatrices\n\n";
2 
3 #include <petscis.h>
4 #include <petscmat.h>
5 
6 int main(int argc, char **args) {
7   Mat         A, *submats, *submats2;
8   IS         *irow, *icol;
9   PetscInt    i, n;
10   PetscMPIInt rank;
11   PetscViewer matfd, rowfd, colfd;
12   PetscBool   same;
13   char        matfile[PETSC_MAX_PATH_LEN], rowfile[PETSC_MAX_PATH_LEN], colfile[PETSC_MAX_PATH_LEN];
14   char        rankstr[16] = {0};
15 
16   PetscFunctionBeginUser;
17   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
18   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
19 
20   PetscCall(PetscOptionsGetString(NULL, NULL, "-A", matfile, sizeof(matfile), NULL));
21   PetscCall(PetscOptionsGetString(NULL, NULL, "-row", rowfile, sizeof(rowfile), NULL));
22   PetscCall(PetscOptionsGetString(NULL, NULL, "-col", colfile, sizeof(colfile), NULL));
23 
24   /* Each rank has its own files for row/col ISes */
25   PetscCall(PetscSNPrintf(rankstr, 16, "-%d", rank));
26   PetscCall(PetscStrlcat(rowfile, rankstr, PETSC_MAX_PATH_LEN));
27   PetscCall(PetscStrlcat(colfile, rankstr, PETSC_MAX_PATH_LEN));
28 
29   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, matfile, FILE_MODE_READ, &matfd));
30   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, rowfile, FILE_MODE_READ, &rowfd));
31   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, colfile, FILE_MODE_READ, &colfd));
32 
33   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
34   PetscCall(MatSetFromOptions(A));
35   PetscCall(MatLoad(A, matfd));
36 
37   /* We stored the number of ISes at the beginning of rowfd */
38   PetscCall(PetscViewerBinaryRead(rowfd, &n, 1, NULL, PETSC_INT));
39   PetscCall(PetscMalloc2(n, &irow, n, &icol));
40   for (i = 0; i < n; i++) {
41     PetscCall(ISCreate(PETSC_COMM_SELF, &irow[i]));
42     PetscCall(ISCreate(PETSC_COMM_SELF, &icol[i]));
43     PetscCall(ISLoad(irow[i], rowfd));
44     PetscCall(ISLoad(icol[i], colfd));
45   }
46 
47   PetscCall(PetscViewerDestroy(&matfd));
48   PetscCall(PetscViewerDestroy(&rowfd));
49   PetscCall(PetscViewerDestroy(&colfd));
50 
51   /* Create submats for the first time */
52   PetscCall(MatCreateSubMatrices(A, n, irow, icol, MAT_INITIAL_MATRIX, &submats));
53 
54   /* Dup submats to submats2 for later comparison */
55   PetscCall(PetscMalloc1(n, &submats2));
56   for (i = 0; i < n; i++) { PetscCall(MatDuplicate(submats[i], MAT_COPY_VALUES, &submats2[i])); }
57 
58   /* Create submats again */
59   PetscCall(MatCreateSubMatrices(A, n, irow, icol, MAT_REUSE_MATRIX, &submats));
60 
61   /* Compare submats and submats2 */
62   for (i = 0; i < n; i++) {
63     PetscCall(MatEqual(submats[i], submats2[i], &same));
64     PetscCheck(same, PETSC_COMM_SELF, PETSC_ERR_PLIB, "submatrix %" PetscInt_FMT " is not same", i);
65   }
66 
67   PetscCall(MatDestroy(&A));
68   for (i = 0; i < n; i++) {
69     PetscCall(ISDestroy(&irow[i]));
70     PetscCall(ISDestroy(&icol[i]));
71   }
72   PetscCall(MatDestroySubMatrices(n, &submats));
73   PetscCall(MatDestroyMatrices(n, &submats2));
74   PetscCall(PetscFree2(irow, icol));
75   PetscCall(PetscFinalize());
76   return 0;
77 }
78 
79 /*TEST
80 
81    test:
82      suffix: 1
83      nsize: 2
84      requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
85      args: -mat_type {{aij baij}} -A ${DATAFILESPATH}/matrices/CreateSubMatrices/A -row ${DATAFILESPATH}/matrices/CreateSubMatrices/row -col ${DATAFILESPATH}/matrices/CreateSubMatrices/col
86 
87 TEST*/
88