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