1 static char help[] = "Tests MatCopy() and MatStore/RetrieveValues().\n\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 Mat C, A;
8 PetscInt i, n = 10, midx[3], bs = 1;
9 PetscScalar v[3];
10 PetscBool flg, isAIJ;
11 MatType type;
12 PetscMPIInt size;
13
14 PetscFunctionBeginUser;
15 PetscCall(PetscInitialize(&argc, &args, NULL, help));
16 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
17 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
18 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL));
19
20 PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
21 PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, n, n));
22 PetscCall(MatSetType(C, MATAIJ));
23 PetscCall(MatSetFromOptions(C));
24 PetscCall(PetscObjectSetName((PetscObject)C, "initial"));
25
26 PetscCall(MatGetType(C, &type));
27 if (size == 1) {
28 PetscCall(PetscObjectTypeCompare((PetscObject)C, MATSEQAIJ, &isAIJ));
29 } else {
30 PetscCall(PetscObjectTypeCompare((PetscObject)C, MATMPIAIJ, &isAIJ));
31 }
32 PetscCall(MatSeqAIJSetPreallocation(C, 3, NULL));
33 PetscCall(MatMPIAIJSetPreallocation(C, 3, NULL, 3, NULL));
34 PetscCall(MatSeqBAIJSetPreallocation(C, bs, 3, NULL));
35 PetscCall(MatMPIBAIJSetPreallocation(C, bs, 3, NULL, 3, NULL));
36 PetscCall(MatSeqSBAIJSetPreallocation(C, bs, 3, NULL));
37 PetscCall(MatMPISBAIJSetPreallocation(C, bs, 3, NULL, 3, NULL));
38
39 v[0] = -1.;
40 v[1] = 2.;
41 v[2] = -1.;
42 for (i = 1; i < n - 1; i++) {
43 midx[2] = i - 1;
44 midx[1] = i;
45 midx[0] = i + 1;
46 PetscCall(MatSetValues(C, 1, &i, 3, midx, v, INSERT_VALUES));
47 }
48 i = 0;
49 midx[0] = 0;
50 midx[1] = 1;
51 v[0] = 2.0;
52 v[1] = -1.;
53 PetscCall(MatSetValues(C, 1, &i, 2, midx, v, INSERT_VALUES));
54 i = n - 1;
55 midx[0] = n - 2;
56 midx[1] = n - 1;
57 v[0] = -1.0;
58 v[1] = 2.;
59 PetscCall(MatSetValues(C, 1, &i, 2, midx, v, INSERT_VALUES));
60
61 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
62 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
63 PetscCall(MatView(C, NULL));
64 PetscCall(MatViewFromOptions(C, NULL, "-view"));
65
66 /* test matduplicate */
67 PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &A));
68 PetscCall(PetscObjectSetName((PetscObject)A, "duplicate_copy"));
69 PetscCall(MatViewFromOptions(A, NULL, "-view"));
70 PetscCall(MatEqual(A, C, &flg));
71 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatDuplicate(C,MAT_COPY_VALUES,): Matrices are NOT equal");
72 PetscCall(MatDestroy(&A));
73
74 /* test matrices with different nonzero patterns - Note: A is created with different nonzero pattern of C! */
75 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
76 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
77 PetscCall(MatSetFromOptions(A));
78 PetscCall(MatSetUp(A));
79 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
80 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
81
82 PetscCall(MatCopy(C, A, DIFFERENT_NONZERO_PATTERN));
83 PetscCall(PetscObjectSetName((PetscObject)A, "copy_diffnnz"));
84 PetscCall(MatViewFromOptions(A, NULL, "-view"));
85 PetscCall(MatEqual(A, C, &flg));
86 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatCopy(C,A,DIFFERENT_NONZERO_PATTERN): Matrices are NOT equal");
87
88 /* test matrices with same nonzero pattern */
89 PetscCall(MatDestroy(&A));
90 PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &A));
91 PetscCall(MatCopy(C, A, SAME_NONZERO_PATTERN));
92 PetscCall(PetscObjectSetName((PetscObject)A, "copy_samennz"));
93 PetscCall(MatViewFromOptions(A, NULL, "-view"));
94 PetscCall(MatEqual(A, C, &flg));
95 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatCopy(C,A,SAME_NONZERO_PATTERN): Matrices are NOT equal");
96
97 /* test subset nonzero pattern */
98 PetscCall(MatCopy(C, A, SUBSET_NONZERO_PATTERN));
99 PetscCall(PetscObjectSetName((PetscObject)A, "copy_subnnz"));
100 PetscCall(MatViewFromOptions(A, NULL, "-view"));
101 PetscCall(MatEqual(A, C, &flg));
102 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatCopy(C,A,SUBSET_NONZERO_PATTERN): Matrices are NOT equal");
103
104 /* Test MatCopy on a matrix obtained after MatConvert from AIJ
105 see https://lists.mcs.anl.gov/pipermail/petsc-dev/2019-April/024289.html */
106 PetscCall(MatHasCongruentLayouts(C, &flg));
107 if (flg) {
108 Mat Cs, Cse;
109 MatType Ctype, Cstype;
110
111 PetscCall(MatGetType(C, &Ctype));
112 PetscCall(MatTranspose(C, MAT_INITIAL_MATRIX, &Cs));
113 PetscCall(MatAXPY(Cs, 1.0, C, DIFFERENT_NONZERO_PATTERN));
114 PetscCall(MatConvert(Cs, MATAIJ, MAT_INPLACE_MATRIX, &Cs));
115 PetscCall(MatSetOption(Cs, MAT_SYMMETRIC, PETSC_TRUE));
116 PetscCall(MatGetType(Cs, &Cstype));
117
118 PetscCall(PetscObjectSetName((PetscObject)Cs, "symm_initial"));
119 PetscCall(MatViewFromOptions(Cs, NULL, "-view"));
120
121 PetscCall(MatConvert(Cs, Ctype, MAT_INITIAL_MATRIX, &Cse));
122 PetscCall(PetscObjectSetName((PetscObject)Cse, "symm_conv_init"));
123 PetscCall(MatViewFromOptions(Cse, NULL, "-view"));
124 PetscCall(MatMultEqual(Cs, Cse, 5, &flg));
125 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatConvert MAT_INITIAL_MATRIX %s -> %s: Matrices are NOT multequal", Ctype, Cstype);
126
127 PetscCall(MatConvert(Cs, Ctype, MAT_REUSE_MATRIX, &Cse));
128 PetscCall(PetscObjectSetName((PetscObject)Cse, "symm_conv_reuse"));
129 PetscCall(MatViewFromOptions(Cse, NULL, "-view"));
130 PetscCall(MatMultEqual(Cs, Cse, 5, &flg));
131 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatConvert MAT_REUSE_MATRIX %s -> %s: Matrices are NOT multequal", Ctype, Cstype);
132
133 PetscCall(MatCopy(Cs, Cse, SAME_NONZERO_PATTERN));
134 PetscCall(PetscObjectSetName((PetscObject)Cse, "symm_conv_copy_samennz"));
135 PetscCall(MatViewFromOptions(Cse, NULL, "-view"));
136 PetscCall(MatMultEqual(Cs, Cse, 5, &flg));
137 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatCopy(...SAME_NONZERO_PATTERN) %s -> %s: Matrices are NOT multequal", Ctype, Cstype);
138
139 PetscCall(MatCopy(Cs, Cse, SUBSET_NONZERO_PATTERN));
140 PetscCall(PetscObjectSetName((PetscObject)Cse, "symm_conv_copy_subnnz"));
141 PetscCall(MatViewFromOptions(Cse, NULL, "-view"));
142 PetscCall(MatMultEqual(Cs, Cse, 5, &flg));
143 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatCopy(...SUBSET_NONZERO_PATTERN) %s -> %s: Matrices are NOT multequal", Ctype, Cstype);
144
145 PetscCall(MatCopy(Cs, Cse, DIFFERENT_NONZERO_PATTERN));
146 PetscCall(PetscObjectSetName((PetscObject)Cse, "symm_conv_copy_diffnnz"));
147 PetscCall(MatViewFromOptions(Cse, NULL, "-view"));
148 PetscCall(MatMultEqual(Cs, Cse, 5, &flg));
149 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatCopy(...DIFFERENT_NONZERO_PATTERN) %s -> %s: Matrices are NOT multequal", Ctype, Cstype);
150
151 PetscCall(MatDestroy(&Cse));
152 PetscCall(MatDestroy(&Cs));
153 }
154
155 /* test MatStore/RetrieveValues() */
156 if (isAIJ) {
157 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE));
158 PetscCall(MatStoreValues(A));
159 PetscCall(MatZeroEntries(A));
160 PetscCall(MatRetrieveValues(A));
161 }
162
163 PetscCall(MatDestroy(&C));
164 PetscCall(MatDestroy(&A));
165 PetscCall(PetscFinalize());
166 return 0;
167 }
168
169 /*TEST
170
171 testset:
172 nsize: {{1 2}separate output}
173 args: -view ::ascii_info -mat_type {{aij baij sbaij mpiaij mpibaij mpisbaij}separate output} -mat_block_size {{1 2}separate output}
174
175 TEST*/
176