Actual source code: ex84.c
1: #include <petscmat.h>
3: #define NNORMS 6
5: static PetscErrorCode MatLoadComputeNorms(Mat data_mat, PetscViewer inp_viewer, PetscReal norms[])
6: {
7: Mat corr_mat;
8: PetscInt M, N;
10: PetscFunctionBegin;
11: PetscCall(MatLoad(data_mat, inp_viewer));
12: PetscCall(MatAssemblyBegin(data_mat, MAT_FINAL_ASSEMBLY));
13: PetscCall(MatAssemblyEnd(data_mat, MAT_FINAL_ASSEMBLY));
14: PetscCall(MatViewFromOptions(data_mat, NULL, "-view_mat"));
16: PetscCall(MatGetSize(data_mat, &M, &N));
17: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Data matrix size: %" PetscInt_FMT " %" PetscInt_FMT "\n", M, N));
19: /* compute matrix norms */
20: PetscCall(MatNorm(data_mat, NORM_1, &norms[0]));
21: PetscCall(MatNorm(data_mat, NORM_INFINITY, &norms[1]));
22: PetscCall(MatNorm(data_mat, NORM_FROBENIUS, &norms[2]));
23: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Data matrix norms: %g %g %g\n", (double)norms[0], (double)norms[1], (double)norms[2]));
25: /* compute autocorrelation matrix */
26: PetscCall(MatMatTransposeMult(data_mat, data_mat, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &corr_mat));
28: /* compute autocorrelation matrix norms */
29: PetscCall(MatNorm(corr_mat, NORM_1, &norms[3]));
30: PetscCall(MatNorm(corr_mat, NORM_INFINITY, &norms[4]));
31: PetscCall(MatNorm(corr_mat, NORM_FROBENIUS, &norms[5]));
32: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Autocorrelation matrix norms: %g %g %g\n", (double)norms[3], (double)norms[4], (double)norms[5]));
34: PetscCall(MatDestroy(&corr_mat));
35: PetscFunctionReturn(PETSC_SUCCESS);
36: }
38: static PetscErrorCode GetReader(MPI_Comm comm, const char option[], PetscViewer *r, PetscViewerFormat *fmt)
39: {
40: PetscBool flg;
42: PetscFunctionBegin;
43: PetscCall(PetscOptionsCreateViewer(PETSC_COMM_SELF, NULL, NULL, option, r, fmt, &flg));
44: if (flg) {
45: PetscFileMode mode;
46: PetscCall(PetscViewerFileGetMode(*r, &mode));
47: flg = (PetscBool)(mode == FILE_MODE_READ);
48: }
49: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Need to specify %s viewer_type:file:format:read", option);
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: int main(int argc, char **argv)
54: {
55: PetscInt i;
56: PetscReal norms0[NNORMS], norms1[NNORMS];
57: PetscViewer inp_viewer;
58: PetscViewerFormat fmt;
59: Mat data_mat;
60: char mat_name[PETSC_MAX_PATH_LEN] = "dmatrix";
62: PetscFunctionBeginUser;
63: PetscCall(PetscInitialize(&argc, &argv, NULL, NULL));
64: PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_name", mat_name, sizeof(mat_name), NULL));
66: /* load matrix sequentially */
67: PetscCall(MatCreate(PETSC_COMM_SELF, &data_mat));
68: PetscCall(MatSetType(data_mat, MATDENSE));
69: PetscCall(PetscObjectSetName((PetscObject)data_mat, mat_name));
70: PetscCall(GetReader(PETSC_COMM_SELF, "-serial_reader", &inp_viewer, &fmt));
71: PetscCall(PetscViewerPushFormat(inp_viewer, fmt));
72: PetscCall(MatLoadComputeNorms(data_mat, inp_viewer, norms0));
73: PetscCall(PetscViewerPopFormat(inp_viewer));
74: PetscCall(PetscViewerDestroy(&inp_viewer));
75: PetscCall(MatViewFromOptions(data_mat, NULL, "-view_serial_mat"));
76: PetscCall(MatDestroy(&data_mat));
78: /* load matrix in parallel */
79: PetscCall(MatCreate(PETSC_COMM_WORLD, &data_mat));
80: PetscCall(MatSetType(data_mat, MATDENSE));
81: PetscCall(PetscObjectSetName((PetscObject)data_mat, mat_name));
82: PetscCall(GetReader(PETSC_COMM_WORLD, "-parallel_reader", &inp_viewer, &fmt));
83: PetscCall(PetscViewerPushFormat(inp_viewer, fmt));
84: PetscCall(MatLoadComputeNorms(data_mat, inp_viewer, norms1));
85: PetscCall(PetscViewerPopFormat(inp_viewer));
86: PetscCall(PetscViewerDestroy(&inp_viewer));
87: PetscCall(MatViewFromOptions(data_mat, NULL, "-view_parallel_mat"));
88: PetscCall(MatDestroy(&data_mat));
90: for (i = 0; i < NNORMS; i++) PetscCheck(PetscAbs(norms0[i] - norms1[i]) <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "norm0[%" PetscInt_FMT "] = %g != %g = norms1[%" PetscInt_FMT "]", i, (double)norms0[i], (double)norms1[i], i);
92: PetscCall(PetscFinalize());
93: return 0;
94: }
96: /*TEST
98: test:
99: suffix: 1
100: requires: hdf5 datafilespath complex
101: args: -serial_reader hdf5:${DATAFILESPATH}/matrices/hdf5/sample_data.h5::read -parallel_reader hdf5:${DATAFILESPATH}/matrices/hdf5/sample_data.h5::read
102: nsize: {{1 2 4}}
104: test:
105: requires: hdf5 datafilespath
106: args: -serial_reader hdf5:${DATAFILESPATH}/matrices/hdf5/tiny_rectangular_mat.h5::read -parallel_reader hdf5:${DATAFILESPATH}/matrices/hdf5/tiny_rectangular_mat.h5::read
107: nsize: {{1 2}}
108: test:
109: suffix: 2-complex
110: requires: complex
111: args: -mat_name ComplexMat
112: test:
113: suffix: 2-real
114: requires: !complex
115: args: -mat_name RealMat
117: TEST*/