1 static char help[] = "Tests various routines in MatMPIBAIJ format.\n";
2
3 #include <petscmat.h>
4 #define IMAX 15
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 Mat A, B, C, At, Bt;
8 PetscViewer fd;
9 char file[PETSC_MAX_PATH_LEN];
10 PetscRandom rand;
11 Vec xx, yy, s1, s2;
12 PetscReal s1norm, s2norm, rnorm, tol = 1.e-10;
13 PetscInt rstart, rend, rows[2], cols[2], m, n, i, j, M, N, ct, row, ncols1, ncols2, bs;
14 PetscMPIInt rank, size;
15 const PetscInt *cols1, *cols2;
16 PetscScalar vals1[4], vals2[4], v;
17 const PetscScalar *v1, *v2;
18 PetscBool flg;
19
20 PetscFunctionBeginUser;
21 PetscCall(PetscInitialize(&argc, &args, NULL, help));
22 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
23 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
24
25 /* Check out if MatLoad() works */
26 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg));
27 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Input file not specified");
28 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
29 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
30 PetscCall(MatSetType(A, MATBAIJ));
31 PetscCall(MatLoad(A, fd));
32
33 PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B));
34 PetscCall(PetscViewerDestroy(&fd));
35
36 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
37 PetscCall(PetscRandomSetFromOptions(rand));
38 PetscCall(MatGetLocalSize(A, &m, &n));
39 PetscCall(VecCreate(PETSC_COMM_WORLD, &xx));
40 PetscCall(VecSetSizes(xx, m, PETSC_DECIDE));
41 PetscCall(VecSetFromOptions(xx));
42 PetscCall(VecDuplicate(xx, &s1));
43 PetscCall(VecDuplicate(xx, &s2));
44 PetscCall(VecDuplicate(xx, &yy));
45 PetscCall(MatGetBlockSize(A, &bs));
46
47 /* Test MatNorm() */
48 PetscCall(MatNorm(A, NORM_FROBENIUS, &s1norm));
49 PetscCall(MatNorm(B, NORM_FROBENIUS, &s2norm));
50 rnorm = PetscAbsScalar(s2norm - s1norm) / s2norm;
51 if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs));
52 PetscCall(MatNorm(A, NORM_INFINITY, &s1norm));
53 PetscCall(MatNorm(B, NORM_INFINITY, &s2norm));
54 rnorm = PetscAbsScalar(s2norm - s1norm) / s2norm;
55 if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs));
56 PetscCall(MatNorm(A, NORM_1, &s1norm));
57 PetscCall(MatNorm(B, NORM_1, &s2norm));
58 rnorm = PetscAbsScalar(s2norm - s1norm) / s2norm;
59 if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatNorm_NORM_1()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs));
60
61 /* Test MatMult() */
62 for (i = 0; i < IMAX; i++) {
63 PetscCall(VecSetRandom(xx, rand));
64 PetscCall(MatMult(A, xx, s1));
65 PetscCall(MatMult(B, xx, s2));
66 PetscCall(VecAXPY(s2, -1.0, s1));
67 PetscCall(VecNorm(s2, NORM_2, &rnorm));
68 if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatMult - Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)rnorm, bs));
69 }
70
71 /* test MatMultAdd() */
72 for (i = 0; i < IMAX; i++) {
73 PetscCall(VecSetRandom(xx, rand));
74 PetscCall(VecSetRandom(yy, rand));
75 PetscCall(MatMultAdd(A, xx, yy, s1));
76 PetscCall(MatMultAdd(B, xx, yy, s2));
77 PetscCall(VecAXPY(s2, -1.0, s1));
78 PetscCall(VecNorm(s2, NORM_2, &rnorm));
79 if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatMultAdd - Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)rnorm, bs));
80 }
81
82 /* Test MatMultTranspose() */
83 for (i = 0; i < IMAX; i++) {
84 PetscCall(VecSetRandom(xx, rand));
85 PetscCall(MatMultTranspose(A, xx, s1));
86 PetscCall(MatMultTranspose(B, xx, s2));
87 PetscCall(VecNorm(s1, NORM_2, &s1norm));
88 PetscCall(VecNorm(s2, NORM_2, &s2norm));
89 rnorm = s2norm - s1norm;
90 if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatMultTranspose - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs));
91 }
92 /* Test MatMultTransposeAdd() */
93 for (i = 0; i < IMAX; i++) {
94 PetscCall(VecSetRandom(xx, rand));
95 PetscCall(VecSetRandom(yy, rand));
96 PetscCall(MatMultTransposeAdd(A, xx, yy, s1));
97 PetscCall(MatMultTransposeAdd(B, xx, yy, s2));
98 PetscCall(VecNorm(s1, NORM_2, &s1norm));
99 PetscCall(VecNorm(s2, NORM_2, &s2norm));
100 rnorm = s2norm - s1norm;
101 if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatMultTransposeAdd - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs));
102 }
103
104 /* Check MatGetValues() */
105 PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
106 PetscCall(MatGetSize(A, &M, &N));
107
108 for (i = 0; i < IMAX; i++) {
109 /* Create random row numbers ad col numbers */
110 PetscCall(PetscRandomGetValue(rand, &v));
111 cols[0] = (int)(PetscRealPart(v) * N);
112 PetscCall(PetscRandomGetValue(rand, &v));
113 cols[1] = (int)(PetscRealPart(v) * N);
114 PetscCall(PetscRandomGetValue(rand, &v));
115 rows[0] = rstart + (int)(PetscRealPart(v) * m);
116 PetscCall(PetscRandomGetValue(rand, &v));
117 rows[1] = rstart + (int)(PetscRealPart(v) * m);
118
119 PetscCall(MatGetValues(A, 2, rows, 2, cols, vals1));
120 PetscCall(MatGetValues(B, 2, rows, 2, cols, vals2));
121
122 for (j = 0; j < 4; j++) {
123 if (vals1[j] != vals2[j]) {
124 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]: Error: MatGetValues rstart = %2" PetscInt_FMT " row = %2" PetscInt_FMT " col = %2" PetscInt_FMT " val1 = %e val2 = %e bs = %" PetscInt_FMT "\n", rank, rstart, rows[j / 2], cols[j % 2], (double)PetscRealPart(vals1[j]), (double)PetscRealPart(vals2[j]), bs));
125 }
126 }
127 }
128
129 /* Test MatGetRow()/ MatRestoreRow() */
130 for (ct = 0; ct < 100; ct++) {
131 PetscCall(PetscRandomGetValue(rand, &v));
132 row = rstart + (PetscInt)(PetscRealPart(v) * m);
133 PetscCall(MatGetRow(A, row, &ncols1, &cols1, &v1));
134 PetscCall(MatGetRow(B, row, &ncols2, &cols2, &v2));
135
136 for (i = 0, j = 0; i < ncols1 && j < ncols2; j++) {
137 while (cols2[j] != cols1[i]) i++;
138 PetscCheck(v1[i] == v2[j], PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetRow() failed - vals incorrect.");
139 }
140 PetscCheck(j >= ncols2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetRow() failed - cols incorrect");
141
142 PetscCall(MatRestoreRow(A, row, &ncols1, &cols1, &v1));
143 PetscCall(MatRestoreRow(B, row, &ncols2, &cols2, &v2));
144 }
145
146 /* Test MatConvert() */
147 PetscCall(MatConvert(A, MATSAME, MAT_INITIAL_MATRIX, &C));
148
149 /* See if MatMult Says both are same */
150 for (i = 0; i < IMAX; i++) {
151 PetscCall(VecSetRandom(xx, rand));
152 PetscCall(MatMult(A, xx, s1));
153 PetscCall(MatMult(C, xx, s2));
154 PetscCall(VecNorm(s1, NORM_2, &s1norm));
155 PetscCall(VecNorm(s2, NORM_2, &s2norm));
156 rnorm = s2norm - s1norm;
157 if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error in MatConvert: MatMult - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs));
158 }
159 PetscCall(MatDestroy(&C));
160
161 /* Test MatTranspose() */
162 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &At));
163 PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &Bt));
164 for (i = 0; i < IMAX; i++) {
165 PetscCall(VecSetRandom(xx, rand));
166 PetscCall(MatMult(At, xx, s1));
167 PetscCall(MatMult(Bt, xx, s2));
168 PetscCall(VecNorm(s1, NORM_2, &s1norm));
169 PetscCall(VecNorm(s2, NORM_2, &s2norm));
170 rnorm = s2norm - s1norm;
171 if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error in MatConvert:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs));
172 }
173 PetscCall(MatDestroy(&At));
174 PetscCall(MatDestroy(&Bt));
175
176 PetscCall(MatDestroy(&A));
177 PetscCall(MatDestroy(&B));
178 PetscCall(VecDestroy(&xx));
179 PetscCall(VecDestroy(&yy));
180 PetscCall(VecDestroy(&s1));
181 PetscCall(VecDestroy(&s2));
182 PetscCall(PetscRandomDestroy(&rand));
183 PetscCall(PetscFinalize());
184 return 0;
185 }
186
187 /*TEST
188
189 build:
190 requires: !complex
191
192 test:
193 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
194 nsize: 3
195 args: -matload_block_size 1 -f ${DATAFILESPATH}/matrices/small
196 output_file: output/empty.out
197
198 test:
199 suffix: 2
200 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
201 nsize: 3
202 args: -matload_block_size 2 -f ${DATAFILESPATH}/matrices/small
203 output_file: output/empty.out
204
205 test:
206 suffix: 3
207 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
208 nsize: 3
209 args: -matload_block_size 4 -f ${DATAFILESPATH}/matrices/small
210 output_file: output/empty.out
211
212 test:
213 TODO: Matrix row/column sizes are not compatible with block size
214 suffix: 4
215 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
216 nsize: 3
217 args: -matload_block_size 5 -f ${DATAFILESPATH}/matrices/small
218 output_file: output/empty.out
219
220 test:
221 suffix: 5
222 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
223 nsize: 3
224 args: -matload_block_size 6 -f ${DATAFILESPATH}/matrices/small
225 output_file: output/empty.out
226
227 test:
228 TODO: Matrix row/column sizes are not compatible with block size
229 suffix: 6
230 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
231 nsize: 3
232 args: -matload_block_size 7 -f ${DATAFILESPATH}/matrices/small
233 output_file: output/empty.out
234
235 test:
236 TODO: Matrix row/column sizes are not compatible with block size
237 suffix: 7
238 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
239 nsize: 3
240 args: -matload_block_size 8 -f ${DATAFILESPATH}/matrices/small
241 output_file: output/empty.out
242
243 test:
244 suffix: 8
245 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
246 nsize: 4
247 args: -matload_block_size 3 -f ${DATAFILESPATH}/matrices/small
248 output_file: output/empty.out
249
250 TEST*/
251