#include <stdio.h>
#include <stdlib.h>
#include <slepceps.h>
#include <petscviewerhdf5.h>

int main(int argc, char **argv){
  PetscErrorCode ierr;
  PetscReal dmat_norm, cmat_norm;
  PetscViewer h5_inp_viewer;
  Mat data_mat, corr_mat;
  const char* input_file = "sample_data.h5";

  ierr = PetscInitialize(&argc, &argv, NULL, NULL); CHKERRQ(ierr);

  // Set up data matrix.
  ierr = MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 882, 10, NULL, &data_mat); CHKERRQ(ierr);
  ierr = MatSetUp(data_mat); CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject)data_mat, "dmatrix"); CHKERRQ(ierr);
  ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, input_file, FILE_MODE_READ, &h5_inp_viewer); CHKERRQ(ierr);
  ierr = MatLoad(data_mat, h5_inp_viewer); CHKERRQ(ierr);
  ierr = PetscViewerDestroy(&h5_inp_viewer); CHKERRQ(ierr);
  ierr = MatAssemblyBegin(data_mat, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  ierr = MatAssemblyEnd(data_mat, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);

  // Get norm. This agrees between sequential and parallel runs.
  ierr = MatNorm(data_mat, NORM_FROBENIUS, &dmat_norm); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD, "Data matrix norm: %8.4e\n", dmat_norm); CHKERRQ(ierr);

  // Create autocorrelation matrix R = D.DT
  ierr = MatMatTransposeMult(data_mat, data_mat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &corr_mat); CHKERRQ(ierr);
  ierr = MatAssemblyBegin(corr_mat, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  ierr = MatAssemblyEnd(corr_mat, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);

  // Get autocorrelation matrix norm. This does NOT at all agree between sequential and parallel - depends on no. of processes.
  ierr = MatNorm(corr_mat, NORM_FROBENIUS, &cmat_norm); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD, "Autocorrelation matrix norm: %8.4e\n", cmat_norm); CHKERRQ(ierr);

  // Clean-up
  ierr = MatDestroy(&data_mat); CHKERRQ(ierr);
  ierr = MatDestroy(&corr_mat); CHKERRQ(ierr);

  ierr = PetscFinalize(); CHKERRQ(ierr);
}