1 #include <petscmat.h>
2
3 #define NNORMS 6
4
MatLoadComputeNorms(Mat data_mat,PetscViewer inp_viewer,PetscReal norms[])5 static PetscErrorCode MatLoadComputeNorms(Mat data_mat, PetscViewer inp_viewer, PetscReal norms[])
6 {
7 Mat corr_mat;
8 PetscInt M, N;
9
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"));
15
16 PetscCall(MatGetSize(data_mat, &M, &N));
17 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Data matrix size: %" PetscInt_FMT " %" PetscInt_FMT "\n", M, N));
18
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]));
24
25 /* compute autocorrelation matrix */
26 PetscCall(MatMatTransposeMult(data_mat, data_mat, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &corr_mat));
27
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]));
33
34 PetscCall(MatDestroy(&corr_mat));
35 PetscFunctionReturn(PETSC_SUCCESS);
36 }
37
GetReader(MPI_Comm comm,const char option[],PetscViewer * r,PetscViewerFormat * fmt)38 static PetscErrorCode GetReader(MPI_Comm comm, const char option[], PetscViewer *r, PetscViewerFormat *fmt)
39 {
40 PetscBool flg;
41
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 }
52
main(int argc,char ** argv)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";
61
62 PetscFunctionBeginUser;
63 PetscCall(PetscInitialize(&argc, &argv, NULL, NULL));
64 PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_name", mat_name, sizeof(mat_name), NULL));
65
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));
77
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));
89
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);
91
92 PetscCall(PetscFinalize());
93 return 0;
94 }
95
96 /*TEST
97
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}}
103
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
116
117 TEST*/
118