/*  Implementation of the Buslowicz's algorithm for inverting a polynomial matrix.
	ref: M Buslowicz "Inversion of Polynomial Matrices", Int. J. Control, 
	1980, vol. 33, No 5, 977-984.  
	Alina Solovyova, Fall 2001 */


# include <sys/types.h>
# include <sys/shm.h>
# include <sys/time.h>
# include <stdlib.h>
# include <sys/ipc.h>
# include <unistd.h>
# include <stdio.h>
# include <iostream.h>
# include <math.h>


#define SHM_MODE (SHM_R | SHM_W)

int rsquared;
int rn;
int rnplus1;
int rnplus1timesrsquared;

int n, r;
int & IndexFunction1( int *matrix, int i, int j);
int & IndexFunction2 (int *matrix, int i, int j, int k);
int & IndexFunction3 (int *matrix, int i, int k, int x, int y);


int main(int arqc, char *argv[])
{
	
	struct timeval starttime;
	struct timeval endtime;
	double totaltime;

 	    
        int *a;
//        int aid;
        int *H;
//	int Hid;
	int *R;
//	int Rid;
	int *G;
//	int Gid;

//	cout<< "\n Enter n (degree of a polynomial)  and r (size of matrix) ";
	cin>>n>>r;


       // alpha and Q are the outputed results

	int *alpha = new int[r*n+1];
 
	
	//Q(k) is r x r matrix;  k=0,..n(r-1)

	int *Q;
	Q=new int[(n*(r-1)+1)*r*r];



	int i,k,ll,kr, min,tr,l,mult;
	int x,y;     //x,y indices for matrix elements

        rsquared = r*r;
        rn = r*n;
        rnplus1 = r*n+1;
        rnplus1timesrsquared = (r*n+1)*r*r;
	
	//Identity Matrix
	int Ident[50][50];     //assume r<50
	for(x=0;x<r;x++)
		for(y=0;y<r;y++)
		{   
			if(x==y)
				Ident[x][y]=1;
   			else 
				Ident[x][y]=0;
		}


	
	// n -  degree of the polynomial (number of H matrices) 
	// r - size of the square matrix H (r x r)



   //********************* Setting up shared memory **************************

	//Each H is (r x r) matrix. First index H(i), x and y - for the matrix

       H= new int[(n+1)*r*r];

//  	Hid=shmget(IPC_PRIVATE, ((n+1)*r*r)*sizeof(int), SHM_MODE);
//	if(Hid < 0){
//		fprintf(stderr, "H: shmget error\n");
//      		exit(1);
//	}
     
       // a(i,k). i=1,..r-1;   k=0,..in

        a=new int [(r+1)*(r*n+1)];

//  	aid=shmget(IPC_PRIVATE, ((r+1)*(r*n+1)*sizeof(int)), SHM_MODE);
//	if(aid < 0){
//		fprintf(stderr, "a: shmget error\n");
//		exit(1);
//	}



	//Each G(i,k) is (r x r) matrix.  

  
	G= new int [(r+1)*(r*n+1)*r*r];

  //	Gid=shmget(IPC_PRIVATE,((r+1)*(r*n+1)*r*r*sizeof(int)), SHM_MODE);
//	if(Gid < 0){
//		fprintf(stderr, "G: shmget error\n");
//		exit(1);
//	}


	//Each R(i,k) is (r x r) matrix.  
	// i=1,..r-1; k=0,..in

  
	R= new int[(r+1)*(r*n+1)*r*r];

//  	Rid=shmget(IPC_PRIVATE,((r+1)*(r*n+1)*r*r*sizeof(int)), SHM_MODE);
//	if(Rid < 0){
//		fprintf(stderr, "R: shmget error\n");
//		exit(1);
//	}



//	a=(int *)shmat(aid,0,0);
//	H=(int *)shmat(Hid,0,0);
//	G=(int *)shmat(Gid,0,0);
//	R=(int *)shmat(Rid,0,0);

//       cout<<"\nEnter the elements of H matrices\n ";
	for(i=0;i<n+1;i++)
	{
//		cout<<"\n H("<<i<<"): ";
		for(x=0;x<r;x++)
			for(y=0;y<r;y++)
				cin>>IndexFunction2(H,i,x,y);
	}
	

	gettimeofday(&starttime, NULL);

    // **************** Calculations **************************
    
    //Everything for i=1

    for (k=0;k<n+1;k++)		//independent 
	{	
		for(x=0;x<r;x++)
			for(y=0;y<r;y++){
				IndexFunction3(G, 1,k,x,y)=IndexFunction2(H,k,x,y); //G{1,k}=H{1,k}
			}
           
			 tr=0;
   
			for(l=0; l<r;l++)
			   tr+=IndexFunction3(G, 1,k,l,l);
			IndexFunction1(a, 1, k)=-tr;
    		for(x=0;x<r;x++)
				for(y=0;y<r;y++)
				{
					IndexFunction3(R, 1,k,x,y)=IndexFunction3(G,1,k,x,y)
									+IndexFunction1(a, 1, k)*Ident[x][y];  
		   
				
				}


	}        
  
	for(i=2;i<r+1;i++)
		for(k=0;k<n*i+1;k++)
			for(x=0;x<r;x++)
				for(y=0;y<r;y++)
				{
					IndexFunction3(G,i,k,x,y)=0;
				}



	for(i=2;i<r+1; i++)  		//depends on i-1
	{

	for(k=0;k<n*i+1;k++)    	//independent
		{
			min=k;
			if(k>n) 
				min=n;
			for (ll=0;ll<min+1;ll++)
        
			{
				kr=k-ll;   
				if(kr<=(i-1)*n)  
				{ 

					for(x=0;x<r;x++)
						for(y=0;y<r;y++)
						{
						mult=0;
						for(l=0;l<r;l++)
						{
				             mult+=IndexFunction2(H,ll,x,l)*IndexFunction3(R,i-1,kr,l,y);
						}
					
 
							IndexFunction3(G,i,k,x,y)+=mult;
					
						}
				}
			}
		

			tr=0;
   
			for(l=0; l<r;l++)
			   tr+=IndexFunction3(G, i,k,l,l);
			IndexFunction1(a,i,k)=-tr/i;
			for(x=0;x<r;x++)
				for(y=0;y<r;y++)
				{

				IndexFunction3(R,i,k,x,y)=IndexFunction3(G,i,k,x,y)
						+IndexFunction1(a, i, k)*Ident[x][y];  

				}
		}
	}


//********************Done with calculations***************************
	




	gettimeofday(&endtime,NULL);
	totaltime=endtime.tv_sec - starttime.tv_sec;
	totaltime+=((endtime.tv_usec-starttime.tv_usec)/1000000.0);
	cout<<"   Total time = "<<totaltime<<endl;
      




//Calculate the output:

// alpha coeff.

	for(k=0;k<r*n+1;k++)
	{
		tr=0;
   
			for(l=0; l<r;l++)
			   tr+=IndexFunction3(G, r,k,l,l);
		alpha[k]=pow(-1,r+1)*tr/r;
	//	cout<<endl<<"k= "<<k<<' '<<"    alpha="<<alpha[k];
	//	cout.flush();
	}
    
//Q matrices
       for(k=0;k<n*(r-1)+1;k++)
       {
  //         cout<<endl<<"k= "<<k;

         for(x=0;x<r;x++)
		 {

           for(y=0;y<r;y++)
		   {
              IndexFunction2(Q,k,x,y)=pow(-1,r+1)*IndexFunction3(R,r-1,k,x,y);
//	        cout<<" "<<IndexFunction2(Q,k,x,y);	
		   }
            
		
		 }
		}


    //to check
		cout<<"Check\n";
	for(k=0;k<r*n+1;k++)
		for(x=0;x<r;x++)
			for(y=0;y<r;y++)
			{

			if((IndexFunction3(G,r,k,x,y)+pow(-1,r)*alpha[k]*Ident[x][y] )!= 0)
			
				cout<<"\n ***** didn't work*******\n";
			}


//	shmdt((char *)a);	
//	shmdt((char *)H);	
//	shmdt((char *)G);	
//	shmdt((char *)R);	



   // *************     deallocations  ****************



	delete [] Q;
	
	delete [] alpha;


//	if(shmctl(aid, IPC_RMID,0))
//		fprintf(stderr, "shmctl -- IPC_RMID error\n");
//	if(shmctl(Hid, IPC_RMID,0))
//		fprintf(stderr, "shmctl -- IPC_RMID error\n");
//	if(shmctl(Rid, IPC_RMID,0))
//		fprintf(stderr, "shmctl -- IPC_RMID error\n");
//	if(shmctl(Gid, IPC_RMID,0))
//		fprintf(stderr, "shmctl -- IPC_RMID error\n");
//	return 0;





}

// function 1 is for a 

int & IndexFunction1( int *matrix, int i, int j)
{
        int index;
        //index=i*(r*n+1)+j;
        index=i*rnplus1+j;
	return matrix[index];

}


//for H and Q i=0,..n+1 ; j=0, ..r; k= 0,..r   


int & IndexFunction2 (int *matrix, int i, int j, int k)
{
   int index;
   //index= i*r*r+j*r+k;
   index= i*rsquared+j*r+k;
   return matrix[index];


}
//function3 is for G and R
int & IndexFunction3(int *matrix, int i, int k, int x, int y)
{
	int index;
	//index=i*(r*n+1)*r*r + k*r*r+x*r+y;
        index=i*rnplus1timesrsquared + k*rsquared+x*r+y;
	return matrix[index];
}

