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: MatLoad(data_mat, inp_viewer);
11: MatAssemblyBegin(data_mat, MAT_FINAL_ASSEMBLY);
12: MatAssemblyEnd(data_mat, MAT_FINAL_ASSEMBLY);
13: MatViewFromOptions(data_mat, NULL, "-view_mat");
15: MatGetSize(data_mat, &M, &N);
16: PetscPrintf(PETSC_COMM_WORLD, "Data matrix size: %" PetscInt_FMT " %" PetscInt_FMT "\n", M, N);
18: /* compute matrix norms */
19: MatNorm(data_mat, NORM_1, &norms[0]);
20: MatNorm(data_mat, NORM_INFINITY, &norms[1]);
21: MatNorm(data_mat, NORM_FROBENIUS, &norms[2]);
22: PetscPrintf(PETSC_COMM_WORLD, "Data matrix norms: %g %g %g\n", (double)norms[0], (double)norms[1], (double)norms[2]);
24: /* compute autocorrelation matrix */
25: MatMatTransposeMult(data_mat, data_mat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &corr_mat);
27: /* compute autocorrelation matrix norms */
28: MatNorm(corr_mat, NORM_1, &norms[3]);
29: MatNorm(corr_mat, NORM_INFINITY, &norms[4]);
30: MatNorm(corr_mat, NORM_FROBENIUS, &norms[5]);
31: PetscPrintf(PETSC_COMM_WORLD, "Autocorrelation matrix norms: %g %g %g\n", (double)norms[3], (double)norms[4], (double)norms[5]);
33: MatDestroy(&corr_mat);
34: return 0;
35: }
37: static PetscErrorCode GetReader(MPI_Comm comm, const char option[], PetscViewer *r, PetscViewerFormat *fmt)
38: {
39: PetscBool flg;
41: PetscOptionsGetViewer(PETSC_COMM_SELF, NULL, NULL, option, r, fmt, &flg);
42: if (flg) {
43: PetscFileMode mode;
44: PetscViewerFileGetMode(*r, &mode);
45: flg = (PetscBool)(mode == FILE_MODE_READ);
46: }
48: return 0;
49: }
51: int main(int argc, char **argv)
52: {
53: PetscInt i;
54: PetscReal norms0[NNORMS], norms1[NNORMS];
55: PetscViewer inp_viewer;
56: PetscViewerFormat fmt;
57: Mat data_mat;
58: char mat_name[PETSC_MAX_PATH_LEN] = "dmatrix";
61: PetscInitialize(&argc, &argv, NULL, NULL);
62: PetscOptionsGetString(NULL, NULL, "-mat_name", mat_name, sizeof(mat_name), NULL);
64: /* load matrix sequentially */
65: MatCreate(PETSC_COMM_SELF, &data_mat);
66: MatSetType(data_mat, MATDENSE);
67: PetscObjectSetName((PetscObject)data_mat, mat_name);
68: GetReader(PETSC_COMM_SELF, "-serial_reader", &inp_viewer, &fmt);
69: PetscViewerPushFormat(inp_viewer, fmt);
70: MatLoadComputeNorms(data_mat, inp_viewer, norms0);
71: PetscViewerPopFormat(inp_viewer);
72: PetscViewerDestroy(&inp_viewer);
73: MatViewFromOptions(data_mat, NULL, "-view_serial_mat");
74: MatDestroy(&data_mat);
76: /* load matrix in parallel */
77: MatCreate(PETSC_COMM_WORLD, &data_mat);
78: MatSetType(data_mat, MATDENSE);
79: PetscObjectSetName((PetscObject)data_mat, mat_name);
80: GetReader(PETSC_COMM_WORLD, "-parallel_reader", &inp_viewer, &fmt);
81: PetscViewerPushFormat(inp_viewer, fmt);
82: MatLoadComputeNorms(data_mat, inp_viewer, norms1);
83: PetscViewerPopFormat(inp_viewer);
84: PetscViewerDestroy(&inp_viewer);
85: MatViewFromOptions(data_mat, NULL, "-view_parallel_mat");
86: MatDestroy(&data_mat);
90: PetscFinalize();
91: return 0;
92: }
94: /*TEST
96: test:
97: suffix: 1
98: requires: hdf5 datafilespath complex
99: args: -serial_reader hdf5:${DATAFILESPATH}/matrices/hdf5/sample_data.h5::read -parallel_reader hdf5:${DATAFILESPATH}/matrices/hdf5/sample_data.h5::read
100: nsize: {{1 2 4}}
102: test:
103: requires: hdf5 datafilespath
104: args: -serial_reader hdf5:${DATAFILESPATH}/matrices/hdf5/tiny_rectangular_mat.h5::read -parallel_reader hdf5:${DATAFILESPATH}/matrices/hdf5/tiny_rectangular_mat.h5::read
105: nsize: {{1 2}}
106: test:
107: suffix: 2-complex
108: requires: complex
109: args: -mat_name ComplexMat
110: test:
111: suffix: 2-real
112: requires: !complex
113: args: -mat_name RealMat
115: TEST*/