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:   PetscReal         norms0[NNORMS], norms1[NNORMS];
 56:   PetscViewer       inp_viewer;
 57:   PetscViewerFormat fmt;
 58:   Mat               data_mat;
 59:   char              mat_name[PETSC_MAX_PATH_LEN] = "dmatrix";

 61:   PetscFunctionBeginUser;
 62:   PetscCall(PetscInitialize(&argc, &argv, NULL, NULL));
 63:   PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_name", mat_name, sizeof(mat_name), NULL));

 65:   /* load matrix sequentially */
 66:   PetscCall(MatCreate(PETSC_COMM_SELF, &data_mat));
 67:   PetscCall(MatSetType(data_mat, MATDENSE));
 68:   PetscCall(PetscObjectSetName((PetscObject)data_mat, mat_name));
 69:   PetscCall(GetReader(PETSC_COMM_SELF, "-serial_reader", &inp_viewer, &fmt));
 70:   PetscCall(PetscViewerPushFormat(inp_viewer, fmt));
 71:   PetscCall(MatLoadComputeNorms(data_mat, inp_viewer, norms0));
 72:   PetscCall(PetscViewerPopFormat(inp_viewer));
 73:   PetscCall(PetscViewerDestroy(&inp_viewer));
 74:   PetscCall(MatViewFromOptions(data_mat, NULL, "-view_serial_mat"));
 75:   PetscCall(MatDestroy(&data_mat));

 77:   /* load matrix in parallel */
 78:   PetscCall(MatCreate(PETSC_COMM_WORLD, &data_mat));
 79:   PetscCall(MatSetType(data_mat, MATDENSE));
 80:   PetscCall(PetscObjectSetName((PetscObject)data_mat, mat_name));
 81:   PetscCall(GetReader(PETSC_COMM_WORLD, "-parallel_reader", &inp_viewer, &fmt));
 82:   PetscCall(PetscViewerPushFormat(inp_viewer, fmt));
 83:   PetscCall(MatLoadComputeNorms(data_mat, inp_viewer, norms1));
 84:   PetscCall(PetscViewerPopFormat(inp_viewer));
 85:   PetscCall(PetscViewerDestroy(&inp_viewer));
 86:   PetscCall(MatViewFromOptions(data_mat, NULL, "-view_parallel_mat"));
 87:   PetscCall(MatDestroy(&data_mat));

 89:   for (PetscInt 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:   PetscCall(PetscFinalize());
 92:   return 0;
 93: }

 95: /*TEST

 97:   test:
 98:     suffix: 1
 99:     requires: hdf5 datafilespath complex
100:     args: -serial_reader hdf5:${DATAFILESPATH}/matrices/hdf5/sample_data.h5::read -parallel_reader hdf5:${DATAFILESPATH}/matrices/hdf5/sample_data.h5::read
101:     nsize: {{1 2 4}}

103:   test:
104:     requires: hdf5 datafilespath
105:     args: -serial_reader hdf5:${DATAFILESPATH}/matrices/hdf5/tiny_rectangular_mat.h5::read -parallel_reader hdf5:${DATAFILESPATH}/matrices/hdf5/tiny_rectangular_mat.h5::read
106:     nsize: {{1 2}}
107:     test:
108:       suffix: 2-complex
109:       requires: complex
110:       args: -mat_name ComplexMat
111:     test:
112:       suffix: 2-real
113:       requires: !complex
114:       args: -mat_name RealMat

116: TEST*/