1 #include <petscmat.h> 2 3 #define NNORMS 6 4 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_DEFAULT, &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 38 static PetscErrorCode GetReader(MPI_Comm comm, const char option[], PetscViewer *r, PetscViewerFormat *fmt) 39 { 40 PetscBool flg; 41 42 PetscFunctionBegin; 43 PetscCall(PetscOptionsGetViewer(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 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