1 static char help[] = "Show MatShift BUG happening after copying a matrix with no rows on a process";
2 /*
3 Contributed by: Eric Chamberland
4 */
5 #include <petscmat.h>
6
7 /* DEFINE this to turn on/off the bug: */
8 #define SET_2nd_PROC_TO_HAVE_NO_LOCAL_LINES
9
main(int argc,char ** args)10 int main(int argc, char **args)
11 {
12 Mat C;
13 PetscInt i, m = 3;
14 PetscMPIInt rank, size;
15 PetscScalar v;
16 Mat lMatA;
17 PetscInt locallines;
18 PetscInt d_nnz[3] = {0, 0, 0};
19 PetscInt o_nnz[3] = {0, 0, 0};
20
21 PetscFunctionBeginUser;
22 PetscCall(PetscInitialize(&argc, &args, NULL, help));
23 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
24 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
25
26 PetscCheck(2 == size, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "Relevant with 2 processes only");
27 PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
28
29 #ifdef SET_2nd_PROC_TO_HAVE_NO_LOCAL_LINES
30 if (0 == rank) {
31 locallines = m;
32 d_nnz[0] = 1;
33 d_nnz[1] = 1;
34 d_nnz[2] = 1;
35 } else {
36 locallines = 0;
37 }
38 #else
39 if (0 == rank) {
40 locallines = m - 1;
41 d_nnz[0] = 1;
42 d_nnz[1] = 1;
43 } else {
44 locallines = 1;
45 d_nnz[0] = 1;
46 }
47 #endif
48
49 PetscCall(MatSetSizes(C, locallines, locallines, m, m));
50 PetscCall(MatSetFromOptions(C));
51 PetscCall(MatXAIJSetPreallocation(C, 1, d_nnz, o_nnz, NULL, NULL));
52
53 v = 2;
54 /* Assembly on the diagonal: */
55 for (i = 0; i < m; i++) PetscCall(MatSetValues(C, 1, &i, 1, &i, &v, ADD_VALUES));
56 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
57 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
58 PetscCall(MatSetOption(C, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
59 PetscCall(MatSetOption(C, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE));
60 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
61 PetscCall(MatConvert(C, MATSAME, MAT_INITIAL_MATRIX, &lMatA));
62 PetscCall(MatView(lMatA, PETSC_VIEWER_STDOUT_WORLD));
63
64 PetscCall(MatShift(lMatA, -1.0));
65
66 PetscCall(MatDestroy(&lMatA));
67 PetscCall(MatDestroy(&C));
68 PetscCall(PetscFinalize());
69 return 0;
70 }
71
72 /*TEST
73
74 test:
75 nsize: 2
76
77 TEST*/
78