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 CHKERRQ(MatLoad(data_mat, inp_viewer)); 12 CHKERRQ(MatAssemblyBegin(data_mat, MAT_FINAL_ASSEMBLY)); 13 CHKERRQ(MatAssemblyEnd(data_mat, MAT_FINAL_ASSEMBLY)); 14 CHKERRQ(MatViewFromOptions(data_mat, NULL, "-view_mat")); 15 16 CHKERRQ(MatGetSize(data_mat, &M, &N)); 17 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "Data matrix size: %" PetscInt_FMT " %" PetscInt_FMT "\n", M,N)); 18 19 /* compute matrix norms */ 20 CHKERRQ(MatNorm(data_mat, NORM_1, &norms[0])); 21 CHKERRQ(MatNorm(data_mat, NORM_INFINITY, &norms[1])); 22 CHKERRQ(MatNorm(data_mat, NORM_FROBENIUS, &norms[2])); 23 CHKERRQ(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 CHKERRQ(MatMatTransposeMult(data_mat, data_mat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &corr_mat)); 27 28 /* compute autocorrelation matrix norms */ 29 CHKERRQ(MatNorm(corr_mat, NORM_1, &norms[3])); 30 CHKERRQ(MatNorm(corr_mat, NORM_INFINITY, &norms[4])); 31 CHKERRQ(MatNorm(corr_mat, NORM_FROBENIUS, &norms[5])); 32 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "Autocorrelation matrix norms: %g %g %g\n", (double)norms[3],(double)norms[4],(double)norms[5])); 33 34 CHKERRQ(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 CHKERRQ(PetscOptionsGetViewer(PETSC_COMM_SELF, NULL, NULL, option, r, fmt, &flg)); 44 if (flg) { 45 PetscFileMode mode; 46 CHKERRQ(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 PetscErrorCode ierr; 56 PetscInt i; 57 PetscReal norms0[NNORMS], norms1[NNORMS]; 58 PetscViewer inp_viewer; 59 PetscViewerFormat fmt; 60 Mat data_mat; 61 char mat_name[PETSC_MAX_PATH_LEN]="dmatrix"; 62 63 ierr = PetscInitialize(&argc, &argv, NULL, NULL);if (ierr) return ierr; 64 CHKERRQ(PetscOptionsGetString(NULL,NULL,"-mat_name",mat_name,sizeof(mat_name),NULL)); 65 66 /* load matrix sequentially */ 67 CHKERRQ(MatCreate(PETSC_COMM_SELF, &data_mat)); 68 CHKERRQ(MatSetType(data_mat,MATDENSE)); 69 CHKERRQ(PetscObjectSetName((PetscObject)data_mat, mat_name)); 70 CHKERRQ(GetReader(PETSC_COMM_SELF, "-serial_reader", &inp_viewer, &fmt)); 71 CHKERRQ(PetscViewerPushFormat(inp_viewer, fmt)); 72 CHKERRQ(MatLoadComputeNorms(data_mat, inp_viewer, norms0)); 73 CHKERRQ(PetscViewerPopFormat(inp_viewer)); 74 CHKERRQ(PetscViewerDestroy(&inp_viewer)); 75 CHKERRQ(MatViewFromOptions(data_mat, NULL, "-view_serial_mat")); 76 CHKERRQ(MatDestroy(&data_mat)); 77 78 /* load matrix in parallel */ 79 CHKERRQ(MatCreate(PETSC_COMM_WORLD, &data_mat)); 80 CHKERRQ(MatSetType(data_mat,MATDENSE)); 81 CHKERRQ(PetscObjectSetName((PetscObject)data_mat, mat_name)); 82 CHKERRQ(GetReader(PETSC_COMM_WORLD, "-parallel_reader", &inp_viewer, &fmt)); 83 CHKERRQ(PetscViewerPushFormat(inp_viewer, fmt)); 84 CHKERRQ(MatLoadComputeNorms(data_mat, inp_viewer, norms1)); 85 CHKERRQ(PetscViewerPopFormat(inp_viewer)); 86 CHKERRQ(PetscViewerDestroy(&inp_viewer)); 87 CHKERRQ(MatViewFromOptions(data_mat, NULL, "-view_parallel_mat")); 88 CHKERRQ(MatDestroy(&data_mat)); 89 90 for (i=0; i<NNORMS; i++) { 91 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); 92 } 93 94 ierr = PetscFinalize(); 95 return ierr; 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