#define NDIMS 3

#include <stdio.h>
#include <stdlib.h>
#include <petscmat.h>

int bin_make(const char*);
int bin_read(const char*);

int main(int argc, char **argv){
  PetscErrorCode ierr;
  const char* bin_file = "./myfile.bin";
  
  ierr = PetscInitialize(&argc, &argv, NULL, NULL); CHKERRQ(ierr);

  ierr = PetscPrintf(PETSC_COMM_WORLD, "Creating binary %s\n", bin_file); CHKERRQ(ierr);
  ierr = bin_make(bin_file); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD, "Reading from binary %s\n", bin_file); CHKERRQ(ierr);
  ierr = bin_read(bin_file); CHKERRQ(ierr);
  ierr = PetscFinalize(); CHKERRQ(ierr);

  return(ierr);
}

int bin_make(const char* bin_file){
  PetscErrorCode ierr;
  PetscViewer bin_viewer;
  PetscInt my_arr[5] = {1, 1, 3, 5, 8};
  Vec my_vec;
  Mat my_mat;
  PetscReal vnorm, mnorm;

  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, bin_file, FILE_MODE_WRITE, &bin_viewer); CHKERRQ(ierr);

  ierr = PetscIntView(5, my_arr, bin_viewer); CHKERRQ(ierr);

  ierr = VecCreate(PETSC_COMM_WORLD, &my_vec); CHKERRQ(ierr);
  ierr = VecSetSizes(my_vec, PETSC_DECIDE, 6); CHKERRQ(ierr);
  ierr = VecSetType(my_vec, VECSTANDARD); CHKERRQ(ierr);
  ierr = VecSetRandom(my_vec, NULL); CHKERRQ(ierr);
  ierr = VecView(my_vec, bin_viewer); CHKERRQ(ierr);

  ierr = MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 5, 5, NULL, &my_mat); CHKERRQ(ierr);
  ierr = MatSetRandom(my_mat, NULL); CHKERRQ(ierr);
  ierr = MatView(my_mat, bin_viewer); CHKERRQ(ierr);

  ierr = VecNorm(my_vec, NORM_2, &vnorm); CHKERRQ(ierr);
  ierr = MatNorm(my_mat, NORM_FROBENIUS, &mnorm); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD, " Vector norm: %.2f. Matrix norm: %.2f.\n", vnorm, mnorm); CHKERRQ(ierr);

  ierr = VecDestroy(&my_vec); CHKERRQ(ierr);
  ierr = MatDestroy(&my_mat); CHKERRQ(ierr);
  ierr = PetscViewerDestroy(&bin_viewer); CHKERRQ(ierr);

  return(ierr);
}

int bin_read(const char* bin_file){
  PetscErrorCode ierr;
  PetscViewer bin_viewer;
  PetscInt my_arr[5];
  PetscReal vnorm, mnorm;
  Vec my_vec;
  Mat my_mat;
  int bin_fd;
  PetscBool skip_arr = PETSC_FALSE;
  PetscBool skip_vec = PETSC_FALSE;
  PetscBool skip_mat = PETSC_FALSE;

  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, bin_file, FILE_MODE_READ, &bin_viewer); CHKERRQ(ierr);

  if (skip_arr){
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Skipping PetscInt array\n"); CHKERRQ(ierr);
    ierr = PetscBinarySeek(bin_fd, PETSC_BINARY_INT_SIZE*5, PETSC_BINARY_SEEK_CUR, NULL); CHKERRQ(ierr);
  }
  else{
    ierr = PetscViewerBinaryGetDescriptor(bin_viewer, &bin_fd); CHKERRQ(ierr);
    ierr = PetscBinaryRead(bin_fd, my_arr, 5, NULL, PETSC_INT); CHKERRQ(ierr);
    ierr = PetscIntView(5, my_arr, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
  }

  if (skip_vec){
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Skipping Vec\n"); CHKERRQ(ierr);
    ierr = PetscBinarySeek(bin_fd, PETSC_BINARY_SCALAR_SIZE*6, PETSC_BINARY_SEEK_CUR, NULL); CHKERRQ(ierr);
  }
  else{
    ierr = VecCreate(PETSC_COMM_WORLD, &my_vec); CHKERRQ(ierr);
    ierr = VecLoad(my_vec, bin_viewer); CHKERRQ(ierr);
    ierr = VecNorm(my_vec, NORM_2, &vnorm); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, " Vector norm: %.2f\n", vnorm); CHKERRQ(ierr);
    ierr = VecDestroy(&my_vec); CHKERRQ(ierr);
  }

  if (skip_mat){
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Skipping Mat\n"); CHKERRQ(ierr);
    ierr = PetscBinarySeek(bin_fd, PETSC_BINARY_SCALAR_SIZE*5*5, PETSC_BINARY_SEEK_CUR, NULL); CHKERRQ(ierr);
  }
  else{
    ierr = MatCreate(PETSC_COMM_WORLD, &my_mat); CHKERRQ(ierr);
    ierr = MatLoad(my_mat, bin_viewer); CHKERRQ(ierr);
    ierr = MatNorm(my_mat, NORM_FROBENIUS, &mnorm); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Matrix norm: %.2f.\n", mnorm); CHKERRQ(ierr);
    ierr = MatDestroy(&my_mat); CHKERRQ(ierr);
  }

  ierr = PetscViewerDestroy(&bin_viewer); CHKERRQ(ierr);

  return(ierr);
}