static char help[]="eps";

#include <petscmat.h>
#include <petscvec.h>
#include <../src/mat/impls/sbaij/seq/sbaij.h>
#include <petscblaslapack.h>
#include <petscksp.h>
#include <math.h>
#include <petsc-private/vecimpl.h>
#include <petsc-private/taoimpl.h>
#include "../src/tao/pde_constrained/impls/lcl/lcl.h"
#include <slepceps.h>
#include "mpi.h"

extern PetscErrorCode CkEigenSolutions(PetscInt*,Mat*,PetscReal*,Vec*,PetscInt*,PetscInt*,PetscReal*);
extern PetscErrorCode FormElementStiffness_idx0(PetscScalar,PetscScalar,PetscScalar,PetscScalar,PetscScalar*);
extern PetscErrorCode FormElementStiffness(PetscScalar,PetscScalar,PetscScalar,PetscScalar,PetscScalar*);
extern PetscErrorCode FormElementMass_idx0(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar*);
extern PetscErrorCode FormElementMass(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar*);
extern PetscErrorCode FormRandomElementStiffness(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar*);



#undef __FUNCT__
#define __FUNCT__ "main"
 int main(int argc,char **args)
{
  Mat            K,M,DM,DK,ME,KE,Phi,TPhi,P,disp,vel,acc,f,invM,Kcap,a,b,invKcap,U,mmtop,kktop,a2DM;
  Vec            Iv,*Cv,xr,xi,delpe;
  Vec            random,P0,P1,P2,Pans,dKcap,delp,delu,delv,dela,v1,v1a3,acc1a5,v2,*av1,bacc1,v1a4,acc1a6,a1delu,a2delu,disp1,disp2,acc1,acc2/*(intial condition vector of acceleration)*/;    
  char           file1[PETSC_MAX_PATH_LEN],file2[PETSC_MAX_PATH_LEN];           /* input file name */  
  PetscBool      flg,flgB=PETSC_FALSE;
  PetscErrorCode ierr;
  PetscRandom    rnd;
  PetscBool      preload=PETSC_TRUE,isSymmetric,evecs,isHermitian,isGeneralized,isPositive,*aK;
  PetscScalar    *arrayK,*arrayM,*evecs_array,kr,ki,re,im,*axr;
  PetscMPIInt    size,rank;
  PetscInt       i,j,z,k,nev,il,iu,r,trial=1,change;
  PetscInt       maxit,its,lits,nconv,nini=1,ncon=0,S=100000,Dim=100000;
  PetscLogStage  stages[2];
  PetscReal      vl,vu,abstol=1.e-8;
  PetscScalar    Matlab_E,freq[Dim],total[Dim],average[Dim],REV[trial],Xr[Dim];
  PetscInt       t,/*n=5,m=2*/id[1],idx[6],IDX[3],index[1],idxmcsn[Dim],idxmcsm[1],idxrandom[trial],idx1[Dim],idxnt[1000],rows[Dim],one[1],istart,iend,nrows;
  /*S = (n+1)*(m);*/
  ST		 st;	
  EPS		 eps;
  EPSType        type;	
  PetscScalar    Ke[16], Me[16],K0[4],M0[4];
  PetscScalar    width = 25.65/1000;   /*25.65 mm */        
  PetscScalar    depth = 3.25/1000;    /* 3.25 mm   */        
  PetscScalar    Area = width*depth, rho   = 2662;         /*2662 kg/m3*/   
  PetscScalar 	 Tlel =.31;            /*length of the beam*/   
  PetscScalar    lel  = Tlel/5;        /*length of each element*/  
  PetscScalar    I = width*depth*depth*depth/12;          /*moment of inertia*/
  PetscScalar    E = 50000000000,En;       /*Modulus of elasticity*/
  PetscViewer    fd1,fd2,fd3,viewer;
  PetscReal      tol;
  	/*parameters used in Newmarks method*/  	
  PetscScalar 	 value[1000],*aP0,*aDM,*aDK,*aa,*ab,*aKcap,*aav1,*abacc1,*aP2,*adelp,*adelu,*av1a3,*av1a4,*aacc1a5,*aacc1a6,*adelv,*adela,*aa1delu,*aa2delu,*adisp1,*adisp2,*arv1,*arv2,*aacc1,*aacc2;
  PetscScalar	 ainvKcap[S],*adKcap;
  PetscScalar	 gaama=1/2,beta=1/4,dt=.01,tt=1000,a1=200,a2= 40000,a3=400,a4= 2,a5= 2,a6=0;






SlepcInitialize(&argc,&args,(char*)0,help);
 
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
     
/*Setting up socket connection for Petsc-Matlab interfacing*/	
	/*fd1 = PETSC_VIEWER_SOCKET_WORLD;
	/*fd2 = PETSC_VIEWER_SOCKET_2;
	/*fd3 = PETSC_VIEWER_SOCKET_3;
	PetscViewerSocketOpen(PETSC_COMM_WORLD,"NULL",PETSC_DEFAULT,&fd1);*/
 	
	/*PetscViewerSocketOpen(PETSC_COMM_WORLD,"NULL",PETSC_DEFAULT,&fd3);

/*if matrices are to be read from the file*/
	/*ierr = PetscOptionsGetString(NULL,"-file1",file1,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
	ierr = PetscPrintf(PETSC_COMM_WORLD," Reading Real 'Mass' matrix from a binary file...\n");CHKERRQ(ierr);
	ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file1,FILE_MODE_READ,&viewer);CHKERRQ(ierr);*/
	/*ierr = MatCreate(PETSC_COMM_WORLD,&mmtop);CHKERRQ(ierr);
	ierr = MatLoad(mmtop,fd1);CHKERRQ(ierr);
	/*ierr = MatView(mmtop,fd1);CHKERRQ(ierr);
 	
	/*ierr = MatSetFromOptions(Mfile);CHKERRQ(ierr);
 	ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);*/

	/*ierr = PetscOptionsGetString(NULL,"-file2",file2,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
	ierr = PetscPrintf(PETSC_COMM_WORLD," Reading Real 'Stiffness' matrix from a binary file...\n");CHKERRQ(ierr);
	ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file2,FILE_MODE_READ,&fd);CHKERRQ(ierr);*/
	/*PetscViewerSocketOpen(PETSC_COMM_WORLD,"NULL",PETSC_DEFAULT,&fd2); */	
	/*ierr = MatCreate(PETSC_COMM_WORLD,&kktop);CHKERRQ(ierr);
	ierr = MatLoad(kktop,fd1);CHKERRQ(ierr);
	/*ierr = MatView(kktop,fd1);CHKERRQ(ierr);
	/*ierr = MatSetFromOptions(Kfile);CHKERRQ(ierr);
 	ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);*/





 
 
 PetscRandomCreate(PETSC_COMM_WORLD,&rnd);
 PetscRandomSetType(rnd,PETSCRAND48);
 PetscRandomSetFromOptions(rnd);
 
 for(i=0;i<trial;i++) { PetscRandomGetValue(rnd,&Matlab_E);	REV[i]=Matlab_E-0.5;	}
/* for(i=0;i<trial/5;i++) printf("from processor %d:   %f   %f   %f    %f   %f    \n ",rank,REV[i],REV[trial/5+i],REV[2*trial/5+i],REV[3*trial/5+i],REV[4*trial/5+i]);*/

 PetscRandomDestroy(&rnd);      

for(i=0;i<S;i++) {freq[i]=0;total[i]=0;average[i]=0;}

/*_____________________________________________________________________ Creating mass matrix_______________________________________________________*/
  ierr  = MatCreate(PETSC_COMM_WORLD,&M);CHKERRQ(ierr);
  ierr  = MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,S,S);CHKERRQ(ierr);
  ierr  = MatSetFromOptions(M);CHKERRQ(ierr);
  ierr  = MatSetUp(M);CHKERRQ(ierr);
  /*
     Assemble matrix
  */
 
  
  for (i=0; i<5; i++) 
  {
    
    /* node numbers for the four corners of element */
    
    if(i==0)
    { ierr = FormElementMass_idx0(rho,lel,width,depth,M0);
    IDX[0]=0;
    IDX[1]=1;
    IDX[2]=2;
    ierr   = MatSetValues(M,3,IDX,3,IDX,M0,ADD_VALUES);CHKERRQ(ierr);}

    else  { ierr = FormElementMass(rho,lel,width,depth,Me);
    idx[0] = 3*(i-1);
    idx[1] = idx[0]+1; idx[2] = idx[1] + 1; idx[3] = idx[2] + 1; idx[4] = idx[3] + 1; idx[5] = idx[4] + 1;
    ierr   = MatSetValues(M,6,idx,6,idx,Me,ADD_VALUES);CHKERRQ(ierr);}
  }
  ierr = MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  /*MatView(M,fd1);
/* printf("Mass Matrix is" "/n" "/n"); 
 */
/*PetscLogStagePush(stages[0]);*/





change=0;
/* ******************************************************************LOOP FOR MONTE CARLO SIMULATION STARTS******************************************************	*/



for(z=0;z<trial;z++)
{ 

 /*printf( " \n \n %d th trial on  process %d out of %d processes is working fine! \n", z,rank+1, size);


  /* _________________________________________________________________________________Stiffness matrices ___________________________________________________________________*/
 
  ierr  = MatCreate(PETSC_COMM_WORLD,&K);CHKERRQ(ierr);
  ierr  = MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,S,S);CHKERRQ(ierr); 
  ierr  = MatSetFromOptions(K);CHKERRQ(ierr);
  ierr  = MatSetUp(K);CHKERRQ(ierr);
  /*
     Assemble matrix
  */
  
  
  for (i=0; i<5; i++) 
  {     
    
 
    if(i==0)
    {
    ierr = FormElementStiffness_idx0(E,Area,lel,I,K0);
    IDX[0]=2*i;
    IDX[1]=2*i+1;
    IDX[2]=2*i+2;
    ierr   = MatSetValues(K,3,IDX,3,IDX,K0,ADD_VALUES);CHKERRQ(ierr);
    }   
   else if(i<4)
   { 
    ierr = FormElementStiffness(E,Area,lel,I,Ke);CHKERRQ(ierr);
    idx[0] = 3*(i-1);
    idx[1] = idx[0]+1; idx[2] = idx[1] + 1; idx[3] = idx[2] + 1; idx[4] = idx[3] + 1; idx[5] = idx[4] + 1;
    ierr   = MatSetValues(K,6,idx,6,idx,Ke,ADD_VALUES);CHKERRQ(ierr);}
  
   else
   { 
    En = E*(1 + REV[z]/10);
    ierr = FormRandomElementStiffness(En,Area,lel,I,Ke);
    idx[0] = 3*(i-1);
    idx[1] = idx[0]+1; idx[2] = idx[1] + 1; idx[3] = idx[2] + 1;  idx[4] = idx[3] + 1; idx[5] = idx[4] + 1;
    ierr   = MatSetValues(K,6,idx,6,idx,Ke,ADD_VALUES);CHKERRQ(ierr);
   }
  }
   
  ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 /*  printf("Mass Matrix is" "/n" "/n"); 
 MatView(K,viewer);*/
				/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      
         									 Create the eigensolver and set various options
  /*
     Create eigensolver context
  */
  
   ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRQ(ierr); 													  

  /*
     Set operators. In this case, it is a generalized eigenvalue problem
  */
  ierr = EPSSetOperators(eps,K,M);CHKERRQ(ierr); 	/*==========================Do necessary changes with Kfile and Mfile according to the location from which they are taken======*/
  MatGetVecs(K,NULL,&xr);
 /* VecSet(xr,100);*/
  /*MatGetVecs(K,NULL,&xi);	*/			/*=========Do necessary changes=======*/		
  /*
     If the user provided initial guesses or constraints, pass them here
  */
  /*VecLoad(Iv,fd1);*/
  ierr = EPSSetInitialSpace(eps,1,&xr);CHKERRQ(ierr);
 /* ierr = EPSSetDeflationSpace(eps,ncon,Cv);CHKERRQ(ierr);*/

  /*
     Set solver parameters at runtime
  */
  ierr = EPSSetFromOptions(eps);CHKERRQ(ierr);

 				 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                   									Solve the eigensystem   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 
   /*set problem as generalized hermitian*/
  										
  EPSSetProblemType(eps,EPS_GHEP); 


  ierr = EPSIsGeneralized(eps,&isGeneralized);CHKERRQ(ierr);
  ierr = EPSIsHermitian(eps,&isHermitian);CHKERRQ(ierr);
  ierr = EPSIsPositive(eps,&isPositive);CHKERRQ(ierr);
  ierr = EPSSolve(eps);CHKERRQ(ierr);

  



   /*ierr = EPSPrintSolution(eps,NULL);CHKERRQ(ierr);*/
 	 
  
  ierr = EPSGetConverged(eps,&nconv);CHKERRQ(ierr);

  /*printf(" converged eigenpairs are %d \n",nconv);*/
	 /* setting up the eigenvector matrix */
	 ierr  = MatCreate(PETSC_COMM_WORLD,&Phi);CHKERRQ(ierr);
 	 ierr  = MatSetSizes(Phi,PETSC_DECIDE,PETSC_DECIDE,Dim,nconv);CHKERRQ(ierr); 
 	 ierr  = MatSetFromOptions(Phi);CHKERRQ(ierr);
  	 ierr  = MatSetUp(Phi);CHKERRQ(ierr);
	
	
	
     /*Save eigenvectors, if requested in Phi matrix*/

	/*initializing Phi matrix*/
	
  for (i=0;i<nconv;i++)
 {    
	
     ierr = EPSGetEigenpair(eps,i,&kr,&ki,xr,xi);CHKERRQ(ierr);
    /*ierr = PetscPrintf(PETSC_COMM_WORLD,"eigenvalues are:   %12f   \n",sqrt((double)kr)/(2*3.1415926535));CHKERRQ(ierr);*/


	VecGetArray(xr,&axr);
	for(j=0;j<Dim;j++){idx1[j]=j;}
	MatSetValues(Phi,Dim,&idx1,1,&i,axr,ADD_VALUES);
	VecRestoreArray(xr,&axr);
	
/*ierr = padmatrix(xr,S,nconv,Phi,i);CHKERRQ(ierr);*/


      re = kr;
      im = ki;
  total[i]= total[i] + kr;
     
	
	
    }
	/*ierr = VecDestroy(&xr);CHKERRQ(ierr);*/
        ierr = MatAssemblyBegin(Phi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 	ierr = MatAssemblyEnd(Phi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
	/*ierr = MatView(Phi,fd1);CHKERRQ(ierr);
   /* ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr);*/
  



  
  
 
  ierr = PetscFree(xi);CHKERRQ(ierr);
  

   ierr = EPSDestroy(&eps);CHKERRQ(ierr);
    ierr = VecDestroy(&xr);CHKERRQ(ierr);

/*_________________________________________________________________________________________________________________________________
 					Solving system dynamics
_________________________________________________________________________________________________________________________________*/


  ierr = MatTranspose(Phi,MAT_INITIAL_MATRIX,&TPhi);CHKERRQ(ierr);      /* B = A^T */
  ierr = MatMatMult(M,Phi,MAT_INITIAL_MATRIX,1,&ME);CHKERRQ(ierr);       /*ME = M*Phi*/
  ierr = MatMatMult(TPhi,ME,MAT_INITIAL_MATRIX,1,&DM);CHKERRQ(ierr);    /*DM = TPhi*ME*/

   
  ierr = MatMatMult(K,Phi,MAT_INITIAL_MATRIX,1,&KE);CHKERRQ(ierr);       /*ME = K*Phi*/
  ierr = MatMatMult(TPhi,KE,MAT_INITIAL_MATRIX,1,&DK);CHKERRQ(ierr);    /*DM = TPhi*KE*/
  /*ierr = PetscPrintf(PETSC_COMM_WORLD,"Decoupled mass matrix is\n");CHKERRQ(ierr);*/
 /* MatView(DM,fd1);
  /*ierr = PetscPrintf(PETSC_COMM_WORLD,"Decoupled stiffness matrix is\n");CHKERRQ(ierr);
  MatView(DK,PETSC_VIEWER_STDOUT_WORLD);*/

    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
    
    /*ierr = MatDestroy(&Phi);CHKERRQ(ierr);*/
    
	ierr = MatDestroy(&KE);CHKERRQ(ierr);
	ierr = MatDestroy(&ME);CHKERRQ(ierr);

/* ************************************************************************** LOOP FOR MONTE CARLO ENDS******************************************************/
}
 for(i=0;i<nconv;i++) {
     average[i]=total[i]/trial;
     PetscPrintf(PETSC_COMM_WORLD,"average eigenvalues are :%f \n",sqrt(average[i])/(2*3.1415926535));
    }



   	


  /* Free work space. */



 

  ierr = MatDestroy(&M);CHKERRQ(ierr);
  

  /*ierr = PetscFinalize();*/
  ierr = SlepcFinalize();
  return 0; 
}