Actual source code: ex53.c
1: static char help[] = "Tests various routines in MatMPIBAIJ format.\n";
3: #include <petscmat.h>
4: #define IMAX 15
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;
20: PetscFunctionBeginUser;
21: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
22: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
23: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
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));
33: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B));
34: PetscCall(PetscViewerDestroy(&fd));
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));
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));
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: }
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: }
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: }
104: /* Check MatGetValues() */
105: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
106: PetscCall(MatGetSize(A, &M, &N));
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);
119: PetscCall(MatGetValues(A, 2, rows, 2, cols, vals1));
120: PetscCall(MatGetValues(B, 2, rows, 2, cols, vals2));
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: }
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));
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");
142: PetscCall(MatRestoreRow(A, row, &ncols1, &cols1, &v1));
143: PetscCall(MatRestoreRow(B, row, &ncols2, &cols2, &v2));
144: }
146: /* Test MatConvert() */
147: PetscCall(MatConvert(A, MATSAME, MAT_INITIAL_MATRIX, &C));
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));
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));
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: }
187: /*TEST
189: build:
190: requires: !complex
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
197: test:
198: suffix: 2
199: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
200: nsize: 3
201: args: -matload_block_size 2 -f ${DATAFILESPATH}/matrices/small
202: output_file: output/ex53_1.out
204: test:
205: suffix: 3
206: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
207: nsize: 3
208: args: -matload_block_size 4 -f ${DATAFILESPATH}/matrices/small
209: output_file: output/ex53_1.out
211: test:
212: TODO: Matrix row/column sizes are not compatible with block size
213: suffix: 4
214: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
215: nsize: 3
216: args: -matload_block_size 5 -f ${DATAFILESPATH}/matrices/small
217: output_file: output/ex53_1.out
219: test:
220: suffix: 5
221: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
222: nsize: 3
223: args: -matload_block_size 6 -f ${DATAFILESPATH}/matrices/small
224: output_file: output/ex53_1.out
226: test:
227: TODO: Matrix row/column sizes are not compatible with block size
228: suffix: 6
229: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
230: nsize: 3
231: args: -matload_block_size 7 -f ${DATAFILESPATH}/matrices/small
232: output_file: output/ex53_1.out
234: test:
235: TODO: Matrix row/column sizes are not compatible with block size
236: suffix: 7
237: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
238: nsize: 3
239: args: -matload_block_size 8 -f ${DATAFILESPATH}/matrices/small
240: output_file: output/ex53_1.out
242: test:
243: suffix: 8
244: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
245: nsize: 4
246: args: -matload_block_size 3 -f ${DATAFILESPATH}/matrices/small
247: output_file: output/ex53_1.out
249: TEST*/