# include "builderh.h"

#undef __FUNCT__  
#define __FUNCT__ "CreateMainMatPar"
PetscErrorCode PETSCKSP_DLLEXPORT CreateMainMatPar(ProblemType problem,PetscInt k,PetscInt v,PetscInt n,PetscInt s,\
		PetscReal h_x,PetscReal h_t,MainMatPar** mmtype) {
	MainMatPar *newmmt;
    PetscErrorCode ierr;
    PetscFunctionBegin;
    PetscValidPointer(mmtype,8);
    
	ierr = PetscNew(MainMatPar,&newmmt);CHKERRQ(ierr);	
  	newmmt->problem = problem;
  	newmmt->k = k;
  	newmmt->v = v;
  	newmmt->n_x = n;
  	newmmt->s_t = s;
  	newmmt->h_x = h_x;
  	newmmt->h_t = h_t;
  	//printf("%d ",newmmt->k);
  	*mmtype       = newmmt;  	
	PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "DestroyMainMatPar"
PetscErrorCode PETSCKSP_DLLEXPORT DestroyMainMatPar(MainMatPar *mmt) {
	PetscErrorCode ierr;
	ierr = PetscFree(mmt);CHKERRQ(ierr);
	PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "BuildMainMatrix"

PetscErrorCode PETSCKSP_DLLEXPORT BuildMainMatrix(Mat *matrix, MainMatPar *mmt) {
	Mat 		A;
	PetscErrorCode ierr;
	
	PetscFunctionBegin;
	PetscValidPointer(matrix,1);
	*matrix = PETSC_NULL;	
	
	if (mmt->problem==RC_1){
		ierr = Build11(&A,mmt);CHKERRQ(ierr);
	}else {
		
	}
	*matrix = A;
	
	PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "Build11"
PetscErrorCode PETSCKSP_DLLEXPORT Build11(Mat *matrix,MainMatPar* mmtype) {
	PetscInt k,v,n_x,s_t;
	PetscReal h_x,h_t;
	
	PetscScalar value[3],jacobi_value[3],alpha_v[2][3],beta_v[2][3];
	PetscInt i[3];
	PetscInt col[k],col_tmp,row;
	PetscErrorCode ierr;
	Mat 		A;
	// Initial 
	k= mmtype->k;
	v= mmtype->v;
	n_x =mmtype->n_x;
	s_t =mmtype->s_t;
	h_x =mmtype->h_x;
	h_t = mmtype->h_t;	
	*matrix = PETSC_NULL;	
		
  ierr = MatCreateSeqAIJ(PETSC_COMM_WORLD,n_x*s_t,n_x*s_t,k*k,PETSC_NULL,&A);CHKERRQ(ierr);
  /* 
     Assemble matrix
		alpha_v = [-1 1 0; 0 -1 1];
		beta_v = [5 8 -1; -1 8 5 ]/12;
  */
  jacobi_value[0] = -h_t*(PetscReal)((n_x+1)*(n_x+1)/PETSC_PI/PETSC_PI); 
  jacobi_value[1] = -2.0*jacobi_value[0]; 
  jacobi_value[2] = jacobi_value[0];
    
  alpha_v[0][0]=-1.0; alpha_v[0][1]=1.0; alpha_v[0][2]=0.0;
  alpha_v[1][0]=0.0; alpha_v[1][1]=-1.0; alpha_v[1][2]=1.0;
  
  beta_v[0][0]=5.0/12.0; beta_v[0][1]=8.0/12.0; beta_v[0][2]=-1.0/12.0;
  beta_v[1][0]=-1.0/12.0; beta_v[1][1]=8.0/12.0; beta_v[1][2]=5.0/12.0;
  
  // The first block line: 0
  value[0] = 1.0;
  for (i[0]=0; i[0]<n_x; i[0]++) {
  	row= i[0];
  	ierr = MatSetValues(A,1,&row,1,&row,value,INSERT_VALUES);CHKERRQ(ierr);
  }
  // The rest blocks
  /////////////////// 1: lines form 1 to v-1 ///////////////////
  for (i[0]=1; i[0]<v; i[0]++) { //block
  	// the first line of the block
  	row++;  	
  	for (i[2]=0;i[2]<k;i[2]++) {
  		// A*I-hB*J => [ -h*b*j[0] a-h*b*j[1] -h*b*j[2] ]
  		col[0] = i[2]*n_x; col[1] = i[2]*n_x+1;   
    	value[0]=alpha_v[i[0]-1][i[2]] + beta_v[i[0]-1][i[2]]*jacobi_value[1];
    	value[1]=beta_v[i[0]-1][i[2]]*jacobi_value[2];
  		//printf("1:row=%d  col[0]=%d  col[1]=%d",row,col[0],col[1]);
    	ierr = MatSetValues(A,1,&row,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
  	}  			
  	// middle lines
  	for (i[1]=1; i[1]<n_x-1; i[1]++) { //each line
  		row++;
  		for (i[2]=0;i[2]<k;i[2]++) {
  			// A*I-hB*J => [ -h*b*j[0] a-h*b*j[1] -h*b*j[2] ]
  			col[0] = i[2]*n_x+i[1]-1; col[1] = i[2]*n_x+i[1]; col[2] = i[2]*n_x+i[1]+1;    		
  			value[0]=beta_v[i[0]-1][i[2]]*jacobi_value[0];
  			value[1]=alpha_v[i[0]-1][i[2]] + beta_v[i[0]-1][i[2]]*jacobi_value[1];
  			value[2]=beta_v[i[0]-1][i[2]]*jacobi_value[2];
    		ierr = MatSetValues(A,1,&row,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
  		}  		
	}
	// the last line of the block
	row++;
	for (i[2]=0;i[2]<k;i[2]++) {
  		// A*I-hB*J => [ -h*b*j[0] a-h*b*j[1] -h*b*j[2] ]
  		col[0] = i[2]*n_x+n_x-2; col[1] = i[2]*n_x+n_x-1; 
    	value[0]=beta_v[i[0]-1][i[2]]*jacobi_value[0];
    	value[1]=alpha_v[i[0]-1][i[2]] + beta_v[i[0]-1][i[2]]*(-jacobi_value[0]);
  		ierr = MatSetValues(A,1,&row,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
  	}
  }
  /////////////////////////2: lines from v to s_t-k+v , total:s_t-k+1 //////////
  for (i[0]=0; i[0]<s_t-k+1; i[0]++) { //block
  	// the first line of the block
  	col_tmp = i[0]*n_x;
  	row++; 	  	
  	for (i[2]=0;i[2]<k;i[2]++) {
  		// A*I-hB*J => [ -h*b*j[0] a-h*b*j[1] -h*b*j[2] ]
  		col[0] = col_tmp+i[2]*n_x; col[1] = col_tmp+i[2]*n_x+1; 
    	value[0]=alpha_v[v-1][i[2]] + beta_v[v-1][i[2]]*jacobi_value[1];
    	value[1]=beta_v[v-1][i[2]]*jacobi_value[2];
  		//printf("2:row=%d  col[0]=%d  col[1]=%d",row,col[0],col[1]);
    	ierr = MatSetValues(A,1,&row,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
  	}  			
  	// middle lines
  	for (i[1]=1; i[1]<n_x-1; i[1]++) { //each line
  		row++;
  		for (i[2]=0;i[2]<k;i[2]++) {
  			// A*I-hB*J => [ -h*b*j[0] a-h*b*j[1] -h*b*j[2] ]
  			col[0] = col_tmp+i[2]*n_x+i[1]-1; 
  			col[1] = col_tmp+i[2]*n_x+i[1]; 
  			col[2] = col_tmp+i[2]*n_x+i[1]+1;    		
  			value[0]=beta_v[v-1][i[2]]*jacobi_value[0];
  			value[1]=alpha_v[v-1][i[2]] +beta_v[v-1][i[2]]*jacobi_value[1];
  			value[2]=beta_v[v-1][i[2]]*jacobi_value[2];
    		ierr = MatSetValues(A,1,&row,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
  		}  		
	}
	// the last line of the block
	row++;
	for (i[2]=0;i[2]<k;i[2]++) {
  		// A*I-hB*J => [ -h*b*j[0] a-h*b*j[1] -h*b*j[2] ]
  		col[0] = col_tmp+i[2]*n_x+n_x-2; col[1] = col_tmp+i[2]*n_x+n_x-1; 
    	value[0]=beta_v[v-1][i[2]]*jacobi_value[0];
    	value[1]=alpha_v[v-1][i[2]] + beta_v[v-1][i[2]]*(-jacobi_value[0]);
  		ierr = MatSetValues(A,1,&row,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
  	}
  }
  ////////////////////////////3: lines from s_t-k+v+1: s_t-1 ,total:k-v-1///////
  for (i[0]=0; i[0]<k-v-1; i[0]++) { //block
  	// the first line of the block
  	col_tmp = (s_t-k)*n_x;
  	row++;  	
  	for (i[2]=0;i[2]<k;i[2]++) {
  		// A*I-hB*J => [ -h*b*j[0] a-h*b*j[1] -h*b*j[2] ]
  		col[0] = col_tmp+i[2]*n_x; col[1] = col_tmp+i[2]*n_x+1; 
    	value[0]=alpha_v[i[0]+v][i[2]] + beta_v[i[0]+v][i[2]]*jacobi_value[1];
    	value[1]=beta_v[i[0]+v][i[2]]*jacobi_value[2];
  		//printf("3:row=%d  col[0]=%d  col[1]=%d",row,col[0],col[1]);
    	ierr = MatSetValues(A,1,&row,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
  	}  			
  	// middle lines
  	for (i[1]=1; i[1]<n_x-1; i[1]++) { //each line
  		row++;
  		for (i[2]=0;i[2]<k;i[2]++) {
  			// A*I-hB*J => [ -h*b*j[0] a-h*b*j[1] -h*b*j[2] ]
  			col[0] = col_tmp+i[2]*n_x+i[1]-1; col[1] = col_tmp+i[2]*n_x+i[1]; col[2] = col_tmp+i[2]*n_x+i[1]+1;    		
  			value[0]=beta_v[i[0]+v][i[2]]*jacobi_value[0];
  			value[1]=alpha_v[i[0]+v][i[2]] + beta_v[i[0]+v][i[2]]*jacobi_value[1];
  			value[2]=beta_v[i[0]+v][i[2]]*jacobi_value[2];
    		ierr = MatSetValues(A,1,&row,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
  		}  		
	}
	// the last line of the block
	row++;
	for (i[2]=0;i[2]<k;i[2]++) {
  		// A*I-hB*J => [ -h*b*j[0] a-h*b*j[1] -h*b*j[2] ]
  		col[0] = col_tmp+i[2]*n_x+n_x-2; col[1] = col_tmp+i[2]*n_x+n_x-1; 
    	value[0]= beta_v[i[0]+v][i[2]]*jacobi_value[0];
    	value[1]=alpha_v[i[0]+v][i[2]] + beta_v[i[0]+v][i[2]]*(-jacobi_value[0]);
  		ierr = MatSetValues(A,1,&row,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
  	}
  }
  //////////////////////////////////////////////////////////////////////////////
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  
  *matrix = A;
  
  PetscFunctionReturn(0);
}

//Mat            A,S; 
//PetscInt		its,S_indexes[k];
//////////////// Peconditioner Matirx ////////////////////////////////////////
  /*
  ierr = MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&S);CHKERRQ(ierr);
  // set S_indexes
  for (i[0]=0;i[0]<k;i[0]++) {
  	if (i[0]<v) {
  		S_indexes[i[0]]=i[0]-v+s_t;
  	}else {
  		S_indexes[i[0]]=i[0]-v;
  	}  	
  }
  // set Peconditioner Matrix
  row = -1;
  for (i[0]=0;i[0]<s_t;i[0]++) {
  	row++;
  	
  	// the first line of the block  	
  	for (i[2]=0;i[2]<k;i[2]++) {
  		// A*I-hB*J => [ -h*b*j[0] a-h*b*j[1] -h*b*j[2] ]
  		col_tmp=S_indexes[i[2]]*n_x;
  		col[0] = col_tmp; col[1] = col_tmp+1; 
    	value[0]=alpha_v[v-1][i[2]] + beta_v[v-1][i[2]]*jacobi_value[1];
    	value[1]=beta_v[v-1][i[2]]*jacobi_value[2];
  		//printf(" row=%d  col[0]=%d  col[1]=%d\n",row,col[0],col[1]);
    	ierr = MatSetValues(S,1,&row,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
  	}    		
  	// middle lines
  	for (i[1]=1; i[1]<n_x-1; i[1]++) { //each line
  		row++;
  		for (i[2]=0;i[2]<k;i[2]++) {
  			// A*I-hB*J => [ -h*b*j[0] a-h*b*j[1] -h*b*j[2] ]
  			col_tmp=S_indexes[i[2]]*n_x;
  			
  			col[0] = col_tmp+i[1]-1; 
  			col[1] = col_tmp+i[1]; 
  			col[2] = col_tmp+i[1]+1;    		
  			value[0]=beta_v[v-1][i[2]]*jacobi_value[0];
  			value[1]=alpha_v[v-1][i[2]] + beta_v[v-1][i[2]]*jacobi_value[1];
  			value[2]=beta_v[v-1][i[2]]*jacobi_value[2];
    		ierr = MatSetValues(S,1,&row,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
  		}  		
	}
	// the last line of the block
	row++;
	for (i[2]=0;i[2]<k;i[2]++) {
  		// A*I-hB*J => [ -h*b*j[0] a-h*b*j[1] -h*b*j[2] ]
  		col_tmp=S_indexes[i[2]]*n_x;
  		col[0] = col_tmp+n_x-2; col[1] = col_tmp+n_x-1; 
    	value[0]=beta_v[v-1][i[2]]*jacobi_value[0];
    	value[1]=alpha_v[v-1][i[2]] + beta_v[v-1][i[2]]*(-jacobi_value[0]);
  		ierr = MatSetValues(S,1,&row,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
  	}  	
  	// change indexes
  	for (i[1]=0;i[1]<k;i[1]++) {
  		if (S_indexes[i[1]]<s_t-1) {
  			S_indexes[i[1]]++;
  		}else{
  			S_indexes[i[1]]=0;
  		}
  	}	
  }
  ierr = MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  */