static char help[] = "Template for Solution of system of equations  Ax=b";


#include <petscmat.h>
#include <petscvec.h>
#include <stdio.h>

#undef __FUNC__
#define __FUNC__ "main"

void constructA(Mat *A,PetscInt N);
void LU(Mat *A,PetscInt n);
void forward(Mat *A,Vec *b,PetscInt n);
void backward(Mat *A,Vec *b,PetscInt n);


int main(int argc,char **args)
{
   FILE * file;
   Mat A,L;
   Vec b,bCholesky;
   PetscInt M,N,n,i,j,k,fdisp;
   MatScalar *a;
   PetscScalar *vecp;
   PetscScalar h,f=1.0,v;

   PetscErrorCode ierr;
   N=10;
   M=N+2;
   n=N*N;
   h = 1/((PetscScalar)N+1);
   f = f*h*h;
   PetscInitialize(&argc,&args,(char *) 0,help);
    fdisp=1; //Only for surf plot display and png filename
   /*
       Create Matrix A,L
   */

   ierr = MatCreate(PETSC_COMM_WORLD,&A); CHKERRQ(ierr);
   ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
   ierr = MatSetFromOptions(A);CHKERRQ(ierr);
   ierr = MatSetUp(A);

   ierr = MatCreate(PETSC_COMM_WORLD,&L); CHKERRQ(ierr);
   ierr = MatSetSizes(L,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
   ierr = MatSetFromOptions(L);CHKERRQ(ierr);
   ierr = MatSetUp(L);

   /*
       Constructing A - a symetric N^2*N^2 matrix
   */
	MatSetType(A,MATSEQDENSE);
    MatSeqDenseSetPreallocation(A,PETSC_NULL);
    MatSetType(L,MATSEQDENSE);
    MatSeqDenseSetPreallocation(L,PETSC_NULL);

    ierr = MatAssemblyBegin(L,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(L,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    constructA(&A,N);
´

   /*
       Create Vector b
   */

   ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr);
   ierr = VecSetSizes(b,PETSC_DECIDE,n);CHKERRQ(ierr);
   ierr = VecSetFromOptions(b);CHKERRQ(ierr);

 

   for (i=0;i<n;i++)
   {
            VecSetValues(b,1,&i,&f,INSERT_VALUES);
          
   }

   VecAssemblyBegin(b);
   VecAssemblyEnd(b);
   


   /*
       Print out A to matlab file
   */


  file = fopen ("ass2.m","w");
  fprintf(file,"clc;close all;clear all;\n");



	MatGetArray(A,&a);
	fprintf(file,"A=[");
	for (i = 0;i<n;i++)
	{
	//	printf("Size: %d\n",sizeof(a[i]));
		for (j=0;j<n;j++)
		{
			fprintf(file,"%1.3f ",a[i*n+j]);

		}
		fprintf(file,";\n");
	}
	fprintf(file,"];\n");
	MatRestoreArray(A,&a);
     /*
       Print out b to matlab file
   */

   VecGetArray(b,&vecp);
   fprintf(file,"b=[");
   for (i=0;i<n;i++)
   {
       fprintf(file,"%1.5f;\n",vecp[i]);
   }
   fprintf(file,"];\n");
   VecRestoreArray(b,&vecp);

  

   /*
       LU-factorization: write function LU for this factorization
    */
    LU(&A,n);

    /*  Write functions forward and backward in order to get solution */

    forward(&A,&b,n);
    backward(&A,&b,n);


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