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*/