/*
  Include "petscksp.h" so that we can use KSP solvers.  Note that this file
  automatically includes:
     petscsys.h    - base PETSc routines   petscvec.h - vectors
     petscmat.h    - matrices              petscpc.h  - preconditioners
     petscis.h     - index sets
     petscviewer.h - viewers

  Note:  The corresponding parallel example is ex23.c
*/
#include <petscksp.h>

int main(int argc, char **args)
{
  PetscInt    nlines = 8;   // lines
  PetscInt    ncolumns = 3; // columns
  PetscInt    random_size = 12, m;
  PetscMPIInt rank;
  PetscMPIInt size;
  Mat         R_full, A, S, R_part;

  // Initialize PETSc
  PetscCall(PetscInitialize(&argc, &args, NULL, NULL));

  PetscCallMPI(MPI_Comm_rank(MPI_COMM_WORLD, &rank));
  PetscCallMPI(MPI_Comm_size(MPI_COMM_WORLD, &size));

  // R_full with all values to zero
  PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, nlines, ncolumns, NULL, &R_full));
  PetscCall(MatView(R_full, PETSC_VIEWER_STDOUT_WORLD));

  // Creating and setting A and S to rand values
  PetscCall(MatDenseGetSubMatrix(R_full, PETSC_DECIDE, nlines / 2, PETSC_DECIDE, PETSC_DECIDE, &R_part));
  PetscCall(MatGetLocalSize(R_part, &m, NULL));
  PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, nlines / 2, random_size, NULL, &A));
  PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, random_size, ncolumns, NULL, &S));
  PetscCall(MatSetRandom(A, NULL));
  PetscCall(MatSetRandom(S, NULL));

  // Computing R_part
  PetscCall(MatMatMult(A, S, MAT_REUSE_MATRIX, PETSC_DECIDE, &R_part));

  // Visualizing R_full
  PetscCall(MatDenseRestoreSubMatrix(R_full, &R_part));
  PetscCall(MatView(R_full, PETSC_VIEWER_STDOUT_WORLD));

  // Destroying matrices
  PetscCall(MatDestroy(&R_part));
  PetscCall(MatDestroy(&R_full));

  PetscCall(PetscFinalize());
  return 0;
}