1*f478cf82SStefano Zampini #include <petscmat.h> 2*f478cf82SStefano Zampini #if defined(PETSC_HAVE_HDF5) 3*f478cf82SStefano Zampini #include <petscviewerhdf5.h> 4*f478cf82SStefano Zampini #endif 5*f478cf82SStefano Zampini 6*f478cf82SStefano Zampini int main(int argc, char **argv) 7*f478cf82SStefano Zampini { 8*f478cf82SStefano Zampini PetscErrorCode ierr; 9*f478cf82SStefano Zampini PetscReal dmat_norm[3], cmat_norm[3]; 10*f478cf82SStefano Zampini PetscViewer inp_viewer; 11*f478cf82SStefano Zampini Mat data_mat, corr_mat; 12*f478cf82SStefano Zampini char file[PETSC_MAX_PATH_LEN],hdf5_name[PETSC_MAX_PATH_LEN]; 13*f478cf82SStefano Zampini PetscBool flg; 14*f478cf82SStefano Zampini 15*f478cf82SStefano Zampini ierr = PetscInitialize(&argc, &argv, NULL, NULL); CHKERRQ(ierr); 16*f478cf82SStefano Zampini ierr = PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg);CHKERRQ(ierr); 17*f478cf82SStefano Zampini if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option"); 18*f478cf82SStefano Zampini ierr = PetscOptionsGetString(NULL,NULL,"-hdf5_name",hdf5_name,sizeof(hdf5_name),&flg);CHKERRQ(ierr); 19*f478cf82SStefano Zampini /* Set up data matrix */ 20*f478cf82SStefano Zampini ierr = MatCreate(PETSC_COMM_WORLD, &data_mat); CHKERRQ(ierr); 21*f478cf82SStefano Zampini ierr = MatSetType(data_mat,MATDENSE); CHKERRQ(ierr); 22*f478cf82SStefano Zampini if (flg) { 23*f478cf82SStefano Zampini ierr = PetscObjectSetName((PetscObject)data_mat, hdf5_name); CHKERRQ(ierr); 24*f478cf82SStefano Zampini ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, file, FILE_MODE_READ, &inp_viewer); CHKERRQ(ierr); 25*f478cf82SStefano Zampini } else { 26*f478cf82SStefano Zampini ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &inp_viewer); CHKERRQ(ierr); 27*f478cf82SStefano Zampini } 28*f478cf82SStefano Zampini ierr = MatLoad(data_mat, inp_viewer); CHKERRQ(ierr); 29*f478cf82SStefano Zampini ierr = PetscViewerDestroy(&inp_viewer); CHKERRQ(ierr); 30*f478cf82SStefano Zampini ierr = MatAssemblyBegin(data_mat, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 31*f478cf82SStefano Zampini ierr = MatAssemblyEnd(data_mat, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 32*f478cf82SStefano Zampini ierr = MatViewFromOptions(data_mat, NULL, "-view_mat");CHKERRQ(ierr); 33*f478cf82SStefano Zampini 34*f478cf82SStefano Zampini ierr = MatNorm(data_mat, NORM_1, &dmat_norm[0]); CHKERRQ(ierr); 35*f478cf82SStefano Zampini ierr = MatNorm(data_mat, NORM_INFINITY, &dmat_norm[1]); CHKERRQ(ierr); 36*f478cf82SStefano Zampini ierr = MatNorm(data_mat, NORM_FROBENIUS, &dmat_norm[2]); CHKERRQ(ierr); 37*f478cf82SStefano Zampini ierr = PetscPrintf(PETSC_COMM_WORLD, "Data matrix norms: %g %g %g\n", (double)dmat_norm[0],(double)dmat_norm[1],(double)dmat_norm[2]); CHKERRQ(ierr); 38*f478cf82SStefano Zampini 39*f478cf82SStefano Zampini /* compute autocorrelation matrix */ 40*f478cf82SStefano Zampini ierr = MatMatTransposeMult(data_mat, data_mat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &corr_mat); CHKERRQ(ierr); 41*f478cf82SStefano Zampini 42*f478cf82SStefano Zampini ierr = MatNorm(corr_mat, NORM_1, &cmat_norm[0]); CHKERRQ(ierr); 43*f478cf82SStefano Zampini ierr = MatNorm(corr_mat, NORM_INFINITY, &cmat_norm[1]); CHKERRQ(ierr); 44*f478cf82SStefano Zampini ierr = MatNorm(corr_mat, NORM_FROBENIUS, &cmat_norm[2]); CHKERRQ(ierr); 45*f478cf82SStefano Zampini ierr = PetscPrintf(PETSC_COMM_WORLD, "Autocorrelation matrix norms: %g %g %g\n", (double)cmat_norm[0],(double)cmat_norm[1],(double)cmat_norm[2]); CHKERRQ(ierr); 46*f478cf82SStefano Zampini 47*f478cf82SStefano Zampini ierr = MatDestroy(&data_mat); CHKERRQ(ierr); 48*f478cf82SStefano Zampini ierr = MatDestroy(&corr_mat); CHKERRQ(ierr); 49*f478cf82SStefano Zampini ierr = PetscFinalize(); CHKERRQ(ierr); 50*f478cf82SStefano Zampini } 51