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 PetscCall(PetscInitialize(&argc, &argv, NULL, NULL)); 63 PetscCall(PetscOptionsGetString(NULL,NULL,"-mat_name",mat_name,sizeof(mat_name),NULL)); 64 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)); 76 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)); 88 89 for (i=0; i<NNORMS; i++) { 90 PetscCheckFalse(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 93 PetscCall(PetscFinalize()); 94 return 0; 95 } 96 97 /*TEST 98 99 test: 100 suffix: 1 101 requires: hdf5 datafilespath complex 102 args: -serial_reader hdf5:${DATAFILESPATH}/matrices/hdf5/sample_data.h5::read -parallel_reader hdf5:${DATAFILESPATH}/matrices/hdf5/sample_data.h5::read 103 nsize: {{1 2 4}} 104 105 test: 106 requires: hdf5 datafilespath 107 args: -serial_reader hdf5:${DATAFILESPATH}/matrices/hdf5/tiny_rectangular_mat.h5::read -parallel_reader hdf5:${DATAFILESPATH}/matrices/hdf5/tiny_rectangular_mat.h5::read 108 nsize: {{1 2}} 109 test: 110 suffix: 2-complex 111 requires: complex 112 args: -mat_name ComplexMat 113 test: 114 suffix: 2-real 115 requires: !complex 116 args: -mat_name RealMat 117 118 TEST*/ 119