Actual source code: ex53.c
2: static char help[] = "Tests various routines in MatMPIBAIJ format.\n";
4: #include <petscmat.h>
5: #define IMAX 15
6: int main(int argc, char **args)
7: {
8: Mat A, B, C, At, Bt;
9: PetscViewer fd;
10: char file[PETSC_MAX_PATH_LEN];
11: PetscRandom rand;
12: Vec xx, yy, s1, s2;
13: PetscReal s1norm, s2norm, rnorm, tol = 1.e-10;
14: PetscInt rstart, rend, rows[2], cols[2], m, n, i, j, M, N, ct, row, ncols1, ncols2, bs;
15: PetscMPIInt rank, size;
16: const PetscInt *cols1, *cols2;
17: PetscScalar vals1[4], vals2[4], v;
18: const PetscScalar *v1, *v2;
19: PetscBool flg;
21: PetscFunctionBeginUser;
22: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
23: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
24: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
26: /* Check out if MatLoad() works */
27: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg));
28: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Input file not specified");
29: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
30: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
31: PetscCall(MatSetType(A, MATBAIJ));
32: PetscCall(MatLoad(A, fd));
34: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B));
35: PetscCall(PetscViewerDestroy(&fd));
37: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
38: PetscCall(PetscRandomSetFromOptions(rand));
39: PetscCall(MatGetLocalSize(A, &m, &n));
40: PetscCall(VecCreate(PETSC_COMM_WORLD, &xx));
41: PetscCall(VecSetSizes(xx, m, PETSC_DECIDE));
42: PetscCall(VecSetFromOptions(xx));
43: PetscCall(VecDuplicate(xx, &s1));
44: PetscCall(VecDuplicate(xx, &s2));
45: PetscCall(VecDuplicate(xx, &yy));
46: PetscCall(MatGetBlockSize(A, &bs));
48: /* Test MatNorm() */
49: PetscCall(MatNorm(A, NORM_FROBENIUS, &s1norm));
50: PetscCall(MatNorm(B, NORM_FROBENIUS, &s2norm));
51: rnorm = PetscAbsScalar(s2norm - s1norm) / s2norm;
52: 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));
53: PetscCall(MatNorm(A, NORM_INFINITY, &s1norm));
54: PetscCall(MatNorm(B, NORM_INFINITY, &s2norm));
55: rnorm = PetscAbsScalar(s2norm - s1norm) / s2norm;
56: 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));
57: PetscCall(MatNorm(A, NORM_1, &s1norm));
58: PetscCall(MatNorm(B, NORM_1, &s2norm));
59: rnorm = PetscAbsScalar(s2norm - s1norm) / s2norm;
60: 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));
62: /* Test MatMult() */
63: for (i = 0; i < IMAX; i++) {
64: PetscCall(VecSetRandom(xx, rand));
65: PetscCall(MatMult(A, xx, s1));
66: PetscCall(MatMult(B, xx, s2));
67: PetscCall(VecAXPY(s2, -1.0, s1));
68: PetscCall(VecNorm(s2, NORM_2, &rnorm));
69: if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatMult - Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)rnorm, bs));
70: }
72: /* test MatMultAdd() */
73: for (i = 0; i < IMAX; i++) {
74: PetscCall(VecSetRandom(xx, rand));
75: PetscCall(VecSetRandom(yy, rand));
76: PetscCall(MatMultAdd(A, xx, yy, s1));
77: PetscCall(MatMultAdd(B, xx, yy, s2));
78: PetscCall(VecAXPY(s2, -1.0, s1));
79: PetscCall(VecNorm(s2, NORM_2, &rnorm));
80: if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatMultAdd - Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)rnorm, bs));
81: }
83: /* Test MatMultTranspose() */
84: for (i = 0; i < IMAX; i++) {
85: PetscCall(VecSetRandom(xx, rand));
86: PetscCall(MatMultTranspose(A, xx, s1));
87: PetscCall(MatMultTranspose(B, xx, s2));
88: PetscCall(VecNorm(s1, NORM_2, &s1norm));
89: PetscCall(VecNorm(s2, NORM_2, &s2norm));
90: rnorm = s2norm - s1norm;
91: 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));
92: }
93: /* Test MatMultTransposeAdd() */
94: for (i = 0; i < IMAX; i++) {
95: PetscCall(VecSetRandom(xx, rand));
96: PetscCall(VecSetRandom(yy, rand));
97: PetscCall(MatMultTransposeAdd(A, xx, yy, s1));
98: PetscCall(MatMultTransposeAdd(B, xx, yy, s2));
99: PetscCall(VecNorm(s1, NORM_2, &s1norm));
100: PetscCall(VecNorm(s2, NORM_2, &s2norm));
101: rnorm = s2norm - s1norm;
102: 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));
103: }
105: /* Check MatGetValues() */
106: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
107: PetscCall(MatGetSize(A, &M, &N));
109: for (i = 0; i < IMAX; i++) {
110: /* Create random row numbers ad col numbers */
111: PetscCall(PetscRandomGetValue(rand, &v));
112: cols[0] = (int)(PetscRealPart(v) * N);
113: PetscCall(PetscRandomGetValue(rand, &v));
114: cols[1] = (int)(PetscRealPart(v) * N);
115: PetscCall(PetscRandomGetValue(rand, &v));
116: rows[0] = rstart + (int)(PetscRealPart(v) * m);
117: PetscCall(PetscRandomGetValue(rand, &v));
118: rows[1] = rstart + (int)(PetscRealPart(v) * m);
120: PetscCall(MatGetValues(A, 2, rows, 2, cols, vals1));
121: PetscCall(MatGetValues(B, 2, rows, 2, cols, vals2));
123: for (j = 0; j < 4; j++) {
124: if (vals1[j] != vals2[j]) {
125: 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));
126: }
127: }
128: }
130: /* Test MatGetRow()/ MatRestoreRow() */
131: for (ct = 0; ct < 100; ct++) {
132: PetscCall(PetscRandomGetValue(rand, &v));
133: row = rstart + (PetscInt)(PetscRealPart(v) * m);
134: PetscCall(MatGetRow(A, row, &ncols1, &cols1, &v1));
135: PetscCall(MatGetRow(B, row, &ncols2, &cols2, &v2));
137: for (i = 0, j = 0; i < ncols1 && j < ncols2; j++) {
138: while (cols2[j] != cols1[i]) i++;
139: PetscCheck(v1[i] == v2[j], PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetRow() failed - vals incorrect.");
140: }
141: PetscCheck(j >= ncols2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetRow() failed - cols incorrect");
143: PetscCall(MatRestoreRow(A, row, &ncols1, &cols1, &v1));
144: PetscCall(MatRestoreRow(B, row, &ncols2, &cols2, &v2));
145: }
147: /* Test MatConvert() */
148: PetscCall(MatConvert(A, MATSAME, MAT_INITIAL_MATRIX, &C));
150: /* See if MatMult Says both are same */
151: for (i = 0; i < IMAX; i++) {
152: PetscCall(VecSetRandom(xx, rand));
153: PetscCall(MatMult(A, xx, s1));
154: PetscCall(MatMult(C, xx, s2));
155: PetscCall(VecNorm(s1, NORM_2, &s1norm));
156: PetscCall(VecNorm(s2, NORM_2, &s2norm));
157: rnorm = s2norm - s1norm;
158: 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));
159: }
160: PetscCall(MatDestroy(&C));
162: /* Test MatTranspose() */
163: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &At));
164: PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &Bt));
165: for (i = 0; i < IMAX; i++) {
166: PetscCall(VecSetRandom(xx, rand));
167: PetscCall(MatMult(At, xx, s1));
168: PetscCall(MatMult(Bt, xx, s2));
169: PetscCall(VecNorm(s1, NORM_2, &s1norm));
170: PetscCall(VecNorm(s2, NORM_2, &s2norm));
171: rnorm = s2norm - s1norm;
172: 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));
173: }
174: PetscCall(MatDestroy(&At));
175: PetscCall(MatDestroy(&Bt));
177: PetscCall(MatDestroy(&A));
178: PetscCall(MatDestroy(&B));
179: PetscCall(VecDestroy(&xx));
180: PetscCall(VecDestroy(&yy));
181: PetscCall(VecDestroy(&s1));
182: PetscCall(VecDestroy(&s2));
183: PetscCall(PetscRandomDestroy(&rand));
184: PetscCall(PetscFinalize());
185: return 0;
186: }
188: /*TEST
190: build:
191: requires: !complex
193: test:
194: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
195: nsize: 3
196: args: -matload_block_size 1 -f ${DATAFILESPATH}/matrices/small
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/ex53_1.out
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/ex53_1.out
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/ex53_1.out
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/ex53_1.out
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/ex53_1.out
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/ex53_1.out
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/ex53_1.out
250: TEST*/