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