// How to compute RARt with A and R as distributed (MPI) matrices ?
//
// ~> g++ -o matRARt.exe matRARt.cpp -lpetsc -lm; mpirun -n 2 matRARt.exe

#include <iostream>
#include <string>
#include "petsc.h"

using namespace std;

int main(int argc,char **argv) {
  if (argc != 2) {cout << "error: need arg = seq or mpi" << endl; return 1;}
  string matType = argv[1]; if (matType != "seq" && matType != "mpi") {cout << "error: need arg = seq or mpi" << endl; return 1;}

  PetscInitialize(&argc, &argv, NULL, NULL);
  int size = 0; MPI_Comm_size(MPI_COMM_WORLD, &size); if (size != 2) {cout << "error: np != 2" << endl; return 1;}
  int rank = 0; MPI_Comm_rank(MPI_COMM_WORLD, &rank);

  PetscInt idx[2] = {0, 1}; PetscScalar val[4] = {1., 0., 0., 1.};

  Mat R;
  if(matType == "mpi") MatCreateAIJ(PETSC_COMM_WORLD, 1, 1, 2, 2, 1, NULL, 1, NULL, &R);
  if(matType == "seq") MatCreateSeqAIJ(PETSC_COMM_SELF, 2, 2, 2, NULL, &R);
  MatSetValues(R, 2, idx, 2, idx, val, ADD_VALUES);
  MatAssemblyBegin(R, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(R, MAT_FINAL_ASSEMBLY);

  Mat A;
  if(matType == "mpi") MatCreateAIJ(PETSC_COMM_WORLD, 1, 1, 2, 2, 1, NULL, 1, NULL, &A);
  if(matType == "seq") MatCreateSeqAIJ(PETSC_COMM_SELF, 2, 2, 2, NULL, &A);
  MatSetValues(A, 2, idx, 2, idx, val, ADD_VALUES);
  MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

  Mat C; MatRARt(A, R, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C);
  if(matType == "mpi")              MatView(C, PETSC_VIEWER_STDOUT_WORLD);
  if(matType == "seq" && rank == 0) MatView(C, PETSC_VIEWER_STDOUT_SELF);

  PetscFinalize();
  return 0;
}