#include <iostream>
#include "petsc.h"
using namespace std;
//***************************************************************************
#include <time.h>

int main(int argc, char **args)
{ 	PetscInitialize(&argc, &args, (char*)0, (char*)0);
	double SizeMatrix=100;
	double NonZeros=2;
	int Iterations=10;
	for (int j=0; j<Iterations; j++)
	{	SizeMatrix=2*SizeMatrix;
		std::cout<<"Matrix nxn, n="<<SizeMatrix<<" ";
		
		Mat Jac1, Jac11;

		MatCreate(PETSC_COMM_WORLD, &Jac1);
 		MatSetSizes(Jac1, PETSC_DECIDE, PETSC_DECIDE, SizeMatrix, SizeMatrix);
		MatSetFromOptions(Jac1);
		MatSeqAIJSetPreallocation(Jac1, NonZeros, PETSC_NULL);
		
		MatCreate(PETSC_COMM_WORLD, &Jac11);
 		MatSetSizes(Jac11, PETSC_DECIDE, PETSC_DECIDE, SizeMatrix, SizeMatrix);
		MatSetFromOptions(Jac11);
		MatSeqAIJSetPreallocation(Jac11, NonZeros, PETSC_NULL);
		
		for (int k=0; k<1; k++)
		{	for (int i=0; i<SizeMatrix; i++)
			{	double value = 1.0;
				int j=i;
				MatSetValues(Jac1, 1, &i, 1, &j, &value, INSERT_VALUES);
			}
		}
		MatAssemblyBegin(Jac1, MAT_FINAL_ASSEMBLY);
		MatAssemblyEnd(Jac1, MAT_FINAL_ASSEMBLY);	
		MatDuplicate(Jac1,MAT_COPY_VALUES,&Jac11);
				
		clock_t start1 = clock();
		for (int k=0; k<1000; k++)
		{	for (int i=0; i<SizeMatrix; i++)
			{	double value = 2.0;
				MatSetValues(Jac11, 1, &i, 1, &i, &value, INSERT_VALUES);
			}
		}
		printf("%f", ((double)clock() - start1) / CLOCKS_PER_SEC);
//******************************************************
		std::cout<<" ";
		Mat Jac2, Jac22;
		MatCreate(PETSC_COMM_WORLD, &Jac2);
 		MatSetSizes(Jac2, PETSC_DECIDE, PETSC_DECIDE, SizeMatrix, SizeMatrix);
		MatSetFromOptions(Jac2);
		MatSeqAIJSetPreallocation(Jac2, NonZeros, PETSC_NULL);
		
		MatCreate(PETSC_COMM_WORLD, &Jac22);
 		MatSetSizes(Jac22, PETSC_DECIDE, PETSC_DECIDE, SizeMatrix, SizeMatrix);
		MatSetFromOptions(Jac22);
		MatSeqAIJSetPreallocation(Jac22, NonZeros, PETSC_NULL);
		
		for (int k=0; k<1; k++)
		{	for (int i=0; i<SizeMatrix-1; i++)
			{	double value = 1.0;
				int j=i-1;
				MatSetValues(Jac2, 1, &i, 1, &j, &value, INSERT_VALUES);
			}
		}
		MatAssemblyBegin(Jac2, MAT_FINAL_ASSEMBLY);
		MatAssemblyEnd(Jac2, MAT_FINAL_ASSEMBLY);	
		MatDuplicate(Jac2,MAT_COPY_VALUES,&Jac22);
				
		clock_t start2 = clock();
		for (int k=0; k<1000; k++)
		{	for (int i=0; i<SizeMatrix; i++)
			{	double value = 2.0;
				MatSetValues(Jac2, 1, &i, 1, &i, &value, INSERT_VALUES);
			}
		}
		printf("%f", ((double)clock() - start2) / CLOCKS_PER_SEC);
//******************************************************
		std::cout<<" ";
		Mat Jac3, Jac33;
		MatCreate(PETSC_COMM_WORLD, &Jac3);
 		MatSetSizes(Jac3, PETSC_DECIDE, PETSC_DECIDE, SizeMatrix, SizeMatrix);
		MatSetFromOptions(Jac3);
		//MatSeqAIJSetPreallocation(Jac2, NonZeros, PETSC_NULL);
		
		MatCreate(PETSC_COMM_WORLD, &Jac33);
 		MatSetSizes(Jac33, PETSC_DECIDE, PETSC_DECIDE, SizeMatrix, SizeMatrix);
		MatSetFromOptions(Jac33);
		//MatSeqAIJSetPreallocation(Jac22, NonZeros, PETSC_NULL);
		
		for (int k=0; k<1; k++)
		{	for (int i=0; i<SizeMatrix-1; i++)
			{	double value = 1.0;
				int j=i-1;
				MatSetValues(Jac3, 1, &i, 1, &j, &value, INSERT_VALUES);
			}
		}
		MatAssemblyBegin(Jac3, MAT_FINAL_ASSEMBLY);
		MatAssemblyEnd(Jac3, MAT_FINAL_ASSEMBLY);	
		MatDuplicate(Jac3,MAT_COPY_VALUES,&Jac33);
				
		clock_t start3 = clock();
		for (int k=0; k<1000; k++)
		{	for (int i=0; i<SizeMatrix; i++)
			{	double value = 2.0;
				MatSetValues(Jac3, 1, &i, 1, &i, &value, INSERT_VALUES);
			}
		}
		printf("%f", ((double)clock() - start3) / CLOCKS_PER_SEC);
//****************************************************************************
		std::cout<<"\n";
	}
	PetscFinalize();
	return 0;
}