/*
  A sample code to test Petsc functions
*/
#include <petscksp.h>
#include <iostream>
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
  Mat            A;        /* linear system matrix */
  PetscInt       i,j,m = 3,n = 3;
  PetscErrorCode ierr;
  PetscScalar    v;


  PetscInitialize(&argc,&args,(char*)0,NULL);

  // create matrix A 
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = MatMPIAIJSetPreallocation(A,3,NULL,3,NULL);CHKERRQ(ierr);
  ierr = MatSeqAIJSetPreallocation(A,3,NULL);CHKERRQ(ierr);

  for (i=0; i<m; i++) {
    for (j=0; j<n; j++) {
      v = i*n + j + 1;
      std::cout << "inserting value " << v << " into position " << i << ", " << j << std::endl;
      MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES); CHKERRQ(ierr);
    }
  }

  // change one entry so it is non singular
  i = 2;
  j = 0;
  v = 8;

  std::cout << "inserting value " << v << " into position " << i << ", " << j << std::endl;
  MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES); CHKERRQ(ierr);

  // assemble matrix
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);


  MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  // Test MatGetLocalSubMatrix
  Mat SubA;
  IS isx;
  PetscInt isx_vals[2];
  isx_vals[0] = 0;
  isx_vals[1] = 1;
  CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD, 2, isx_vals, PETSC_COPY_VALUES, &isx));
  CHKERRQ(MatGetLocalSubMatrix(A, isx, isx, &SubA));  // error


  // cleanup
  ierr = MatDestroy(&A); CHKERRQ(ierr);
  ierr = PetscFinalize();

  return 0;
}