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(0); 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(0); 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++) { 91 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); 92 } 93 94 PetscCall(PetscFinalize()); 95 return 0; 96 } 97 98 /*TEST 99 100 test: 101 suffix: 1 102 requires: hdf5 datafilespath complex 103 args: -serial_reader hdf5:${DATAFILESPATH}/matrices/hdf5/sample_data.h5::read -parallel_reader hdf5:${DATAFILESPATH}/matrices/hdf5/sample_data.h5::read 104 nsize: {{1 2 4}} 105 106 test: 107 requires: hdf5 datafilespath 108 args: -serial_reader hdf5:${DATAFILESPATH}/matrices/hdf5/tiny_rectangular_mat.h5::read -parallel_reader hdf5:${DATAFILESPATH}/matrices/hdf5/tiny_rectangular_mat.h5::read 109 nsize: {{1 2}} 110 test: 111 suffix: 2-complex 112 requires: complex 113 args: -mat_name ComplexMat 114 test: 115 suffix: 2-real 116 requires: !complex 117 args: -mat_name RealMat 118 119 TEST*/ 120