static char help[] = "Practice creating a matrix.\n";
#include "petscmat.h"
#undef __FUNCT__
#define __FUNCT__ "main"

int main(int argc,char **args)
{
	Mat A;
	PetscInt rows = 6, columns = 6;
	PetscErrorCode ierr;
	int myRank, numNodes, i, j;

	PetscInitialize(&argc,&args,(char *)0,help);
	ierr = MPI_Comm_size(PETSC_COMM_WORLD, &numNodes);CHKERRQ(ierr);
	ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &myRank);CHKERRQ(ierr);
	double sqrtNodes = sqrt(numNodes);
	if(myRank == 0)
	{
		printf("numNodes = %d\n", numNodes);
	}
	
	if(sqrtNodes / (int)sqrtNodes == 1)
	{
		if(rows % (int)sqrtNodes != 0 || columns % (int)sqrtNodes !=0)
		{
			printf("ERROR: Matrix cannot be divided into one square block per processor\nMy myRank = %d   rows = %d   columns = %d\n", myRank, rows, columns);
			ierr = PetscFinalize();CHKERRQ(ierr);
			return 0;
		}
	}
	else
	{
		printf("ERROR: Nodes cannot be arranged as a grid");
		ierr = PetscFinalize();CHKERRQ(ierr);
		return 0;
	}
	PetscInt rowsPerNode = rows / (int)sqrtNodes, columnsPerNode = columns / (int)sqrtNodes;
	PetscInt globalRowIDs[rowsPerNode], globalColumnIDs[columnsPerNode];
	PetscScalar myMat[rowsPerNode][columnsPerNode];
	
	for(i=0; i<rowsPerNode; i++)
	{
		for(j=0; j<columnsPerNode; j++)
		{
			myMat[i][j] = i*rowsPerNode + j + 1 + myRank*10;
		}
		globalRowIDs[i] = i + rowsPerNode*(myRank / (int)sqrtNodes);
		globalColumnIDs[i] = i + columnsPerNode*(myRank % (int)sqrtNodes);
	}
	printf("my rank = %d   my glob rows = %d %d %d   my glob cols = %d %d %d\n",myRank,globalRowIDs[0],globalRowIDs[1],globalRowIDs[2],globalColumnIDs[0],globalColumnIDs[1],globalColumnIDs[2]);

	ierr = MatCreate(PETSC_COMM_WORLD, &A);CHKERRQ(ierr); 
	ierr = MatSetSizes(A, rowsPerNode, columnsPerNode, rows, columns);CHKERRQ(ierr);
	ierr = MatSetType(A, MATAIJ);CHKERRQ(ierr);

    ierr = MatSetValues(A, rowsPerNode, globalRowIDs, columnsPerNode, globalColumnIDs, &(myMat[0][0]), INSERT_VALUES);CHKERRQ(ierr);
	MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
	//ierr = PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_DENSE);CHKERRQ(ierr);
	ierr = MatView(A, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

    ierr = MatDestroy(A);CHKERRQ(ierr);
	ierr = PetscFinalize();CHKERRQ(ierr);
	
	return 0;
}