Actual source code: ex16.c
1: static char help[] = "Tests MatDenseGetArray() and MatView()/MatLoad() with binary viewers.\n\n";
3: #include <petscmat.h>
4: #include <petscviewer.h>
6: static PetscErrorCode CheckValues(Mat A, PetscBool one)
7: {
8: const PetscScalar *array;
9: PetscInt M, N, rstart, rend, lda, i, j;
11: PetscFunctionBegin;
12: PetscCall(MatDenseGetArrayRead(A, &array));
13: PetscCall(MatDenseGetLDA(A, &lda));
14: PetscCall(MatGetSize(A, &M, &N));
15: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
16: for (i = rstart; i < rend; i++) {
17: for (j = 0; j < N; j++) {
18: PetscInt ii = i - rstart, jj = j;
19: PetscReal v = (PetscReal)(one ? 1 : (1 + i + j * M));
20: PetscReal w = PetscRealPart(array[ii + jj * lda]);
21: PetscCheck(PetscAbsReal(v - w) <= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Matrix entry (%" PetscInt_FMT ",%" PetscInt_FMT ") should be %g, got %g", i, j, (double)v, (double)w);
22: }
23: }
24: PetscCall(MatDenseRestoreArrayRead(A, &array));
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: #define CheckValuesIJ(A) CheckValues(A, PETSC_FALSE)
29: #define CheckValuesOne(A) CheckValues(A, PETSC_TRUE)
31: int main(int argc, char **args)
32: {
33: Mat A;
34: PetscInt i, j, M = 4, N = 3, rstart, rend;
35: PetscScalar *array;
36: char mattype[256];
37: PetscViewer view;
39: PetscFunctionBeginUser;
40: PetscCall(PetscInitialize(&argc, &args, NULL, help));
41: PetscCall(PetscStrncpy(mattype, MATMPIDENSE, sizeof(mattype)));
42: PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_type", mattype, sizeof(mattype), NULL));
43: /*
44: Create a parallel dense matrix shared by all processors
45: */
46: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, NULL, &A));
47: PetscCall(MatConvert(A, mattype, MAT_INPLACE_MATRIX, &A));
48: /*
49: Set values into the matrix
50: */
51: for (i = 0; i < M; i++) {
52: for (j = 0; j < N; j++) {
53: PetscScalar v = (PetscReal)(1 + i + j * M);
54: PetscCall(MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES));
55: }
56: }
57: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
58: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
59: PetscCall(MatScale(A, 2.0));
60: PetscCall(MatScale(A, 1.0 / 2.0));
62: /*
63: Store the binary matrix to a file
64: */
65: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_WRITE, &view));
66: for (i = 0; i < 2; i++) {
67: PetscCall(MatView(A, view));
68: PetscCall(PetscViewerPushFormat(view, PETSC_VIEWER_NATIVE));
69: PetscCall(MatView(A, view));
70: PetscCall(PetscViewerPopFormat(view));
71: }
72: PetscCall(PetscViewerDestroy(&view));
73: PetscCall(MatDestroy(&A));
75: /*
76: Now reload the matrix and check its values
77: */
78: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_READ, &view));
79: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
80: PetscCall(MatSetType(A, mattype));
81: for (i = 0; i < 4; i++) {
82: if (i > 0) PetscCall(MatZeroEntries(A));
83: PetscCall(MatLoad(A, view));
84: PetscCall(CheckValuesIJ(A));
85: }
86: PetscCall(PetscViewerDestroy(&view));
88: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
89: PetscCall(PetscMalloc1((rend - rstart) * N, &array));
90: for (i = 0; i < (rend - rstart) * N; i++) array[i] = (PetscReal)1;
91: PetscCall(MatDensePlaceArray(A, array));
92: PetscCall(MatScale(A, 2.0));
93: PetscCall(MatScale(A, 1.0 / 2.0));
94: PetscCall(CheckValuesOne(A));
95: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_WRITE, &view));
96: PetscCall(MatView(A, view));
97: PetscCall(MatDenseResetArray(A));
98: PetscCall(PetscFree(array));
99: PetscCall(CheckValuesIJ(A));
100: PetscCall(PetscViewerBinarySetSkipHeader(view, PETSC_TRUE));
101: PetscCall(MatView(A, view));
102: PetscCall(PetscViewerBinarySetSkipHeader(view, PETSC_FALSE));
103: PetscCall(PetscViewerDestroy(&view));
104: PetscCall(MatDestroy(&A));
106: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
107: PetscCall(MatSetType(A, mattype));
108: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_READ, &view));
109: PetscCall(MatLoad(A, view));
110: PetscCall(CheckValuesOne(A));
111: PetscCall(MatZeroEntries(A));
112: PetscCall(PetscViewerBinarySetSkipHeader(view, PETSC_TRUE));
113: PetscCall(MatLoad(A, view));
114: PetscCall(PetscViewerBinarySetSkipHeader(view, PETSC_FALSE));
115: PetscCall(CheckValuesIJ(A));
116: PetscCall(PetscViewerDestroy(&view));
117: PetscCall(MatDestroy(&A));
119: {
120: PetscInt m = PETSC_DECIDE, n = PETSC_DECIDE;
121: PetscCall(PetscSplitOwnership(PETSC_COMM_WORLD, &m, &M));
122: PetscCall(PetscSplitOwnership(PETSC_COMM_WORLD, &n, &N));
123: /* TODO: MatCreateDense requires data!=NULL at all processes! */
124: PetscCall(PetscMalloc1(m * N + 1, &array));
126: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.dat", FILE_MODE_READ, &view));
127: PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, n, M, N, array, &A));
128: PetscCall(MatLoad(A, view));
129: PetscCall(CheckValuesOne(A));
130: PetscCall(PetscViewerBinarySetSkipHeader(view, PETSC_TRUE));
131: PetscCall(MatLoad(A, view));
132: PetscCall(PetscViewerBinarySetSkipHeader(view, PETSC_FALSE));
133: PetscCall(CheckValuesIJ(A));
134: PetscCall(MatDestroy(&A));
135: PetscCall(PetscViewerDestroy(&view));
137: PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, n, M, N, array, &A));
138: PetscCall(CheckValuesIJ(A));
139: PetscCall(MatDestroy(&A));
141: PetscCall(PetscFree(array));
142: }
144: PetscCall(PetscFinalize());
145: return 0;
146: }
148: /*TEST
150: testset:
151: args: -viewer_binary_mpiio 0
152: output_file: output/ex16.out
153: test:
154: suffix: stdio_1
155: nsize: 1
156: args: -mat_type seqdense
157: test:
158: suffix: stdio_2
159: nsize: 2
160: test:
161: suffix: stdio_3
162: nsize: 3
163: test:
164: suffix: stdio_4
165: nsize: 4
166: test:
167: suffix: stdio_5
168: nsize: 5
169: test:
170: requires: cuda
171: args: -mat_type seqdensecuda
172: suffix: stdio_cuda_1
173: nsize: 1
174: test:
175: requires: cuda
176: args: -mat_type mpidensecuda
177: suffix: stdio_cuda_2
178: nsize: 2
179: test:
180: requires: cuda
181: args: -mat_type mpidensecuda
182: suffix: stdio_cuda_3
183: nsize: 3
184: test:
185: requires: cuda
186: args: -mat_type mpidensecuda
187: suffix: stdio_cuda_4
188: nsize: 4
189: test:
190: requires: cuda
191: args: -mat_type mpidensecuda
192: suffix: stdio_cuda_5
193: nsize: 5
195: testset:
196: requires: mpiio
197: args: -viewer_binary_mpiio 1
198: output_file: output/ex16.out
199: test:
200: suffix: mpiio_1
201: nsize: 1
202: test:
203: suffix: mpiio_2
204: nsize: 2
205: test:
206: suffix: mpiio_3
207: nsize: 3
208: test:
209: suffix: mpiio_4
210: nsize: 4
211: test:
212: suffix: mpiio_5
213: nsize: 5
214: test:
215: requires: cuda
216: args: -mat_type mpidensecuda
217: suffix: mpiio_cuda_1
218: nsize: 1
219: test:
220: requires: cuda
221: args: -mat_type mpidensecuda
222: suffix: mpiio_cuda_2
223: nsize: 2
224: test:
225: requires: cuda
226: args: -mat_type mpidensecuda
227: suffix: mpiio_cuda_3
228: nsize: 3
229: test:
230: requires: cuda
231: args: -mat_type mpidensecuda
232: suffix: mpiio_cuda_4
233: nsize: 4
234: test:
235: requires: cuda
236: args: -mat_type mpidensecuda
237: suffix: mpiio_cuda_5
238: nsize: 5
240: TEST*/