xi_rsm(int m,  int n, float a[4][4],   float b[4][4],  float b1[4][4], float e[4][4],
       float (*xb)[4][2],  float  ps[4][2],  float qs[4][2], float cb[2][4],
	   float pa[2][4],     float c[2][4],   float ch[2][4],  int  (*k) [4]   )
{
    int  i, j, j0, i0, k1;
	float ee, g, ps1;

	for(i=1; i<=m; i++)
	{
	    (*k)[i]=n+i;
		cb[1][i]=0;
		for(j=1; j<=m; j++)
		{
		   if(i==j)
		     b[i][j]=1;
		   else
			 b[i][j]=0;

		}
	}


loop: for(j=1; j<=m; j++)
		for(i=1; i<=1; i++)
			{
				pa[i][j]=0;
				for(k1=1; k1<=m; k1++)
					pa[i][j] = pa[i][j] + cb[i][k1] * b[k1][j];
			}
		for(j=i; j<=n; j++)
			for(i=1; i<=1; i++)
			{
				ch[i][j]=0;
				for(k1=1; k1<=m; k1++)
					ch[i][j] = ch[i][j] + pa[i][k1] * a[k1][j];

			}
		ee=0;
		for(j=i; j<=n; j++)
		{
			ch[1][j] = ch[1][j] - c[1][j];
			if(ch[1][j]-ee>0)
			{
				ee=ch[1][j];
				j0=j;
			}		
		}
		if (ee<=0)
		{
		for (i=1; i<=m; i++)
			printf("\nOPTIMAL--SOLUTION\nI=%d, XB[i]=%f\n", (*k)[i], (*xb)[i][1]);
			return (0);
		}
		for (i=1; i<=m; i++)
			ps[i][1]=a[i][j0];
		for(j=1; j<=1; j++)
			for (i=1; i<=m; i++)
			{
				qs[i][j]=0;
				for(k1=1; k1<=m; k1++)
					qs[i][j] = qs[i][j] + b[i][k1] * ps[k1][j];
			}
		for(i=1; i<=m; i++)
			ps[i][1] = qs[i][1];
		ee=1e+10;
		for(i=1; i<=m; i++)
		{
			ps1=ps[i][1];
			if (ps1>0)
			{
				g=(*xb)[i][1]/ps1;
				if (g-ee<=0)
				{
					ee=g;
					i0=i;
				}
			}
		}
		if(ee==1e10)
		{
			printf("\nNO--SOLUTION\n");
			return (-1);	
		}
		(*k)[i0]=j0;
	for(i=1; i<=m; i++)
	{
		for(j=1; j<=m; j++)
		{
			e[i][j]=0;
			if(i==j)
				e[i][j]=1;
		}
	  e[i][i0] = -ps[i][1]/ps[i0][1];
	}
	e[i0][i0]=1/ps[i0][1];
	for(j=1; j<=m; j++)
		for(i=1; i<=m; i++)
		{
			b1[i][j]=0;
			for(k1=1; k1<=m; k1++)
				b1[i][j] = b1[i][j] +e[i][k1] * b[k1][j];
		}
	for(j=1; j<=1; j++)
		for(i=1; i<=m; i++)
		{
			qs[i][j]=0;
			for(k1=1; k1<=m; k1++)
				qs[i][j] = qs[i][j] + e[i][k1] * (*xb)[k1][j];		
		}
	for(i=1; i<=m; i++)
	{
		a[i][j0]=0;
		(*xb)[i][1] = qs[i][1];
		for(j=1; j<=m; j++)
			b[i][j] = b1[i][j];
	}

  a[i0][j0]=1;
  cb[1][i0] = c[1][j0];
  c[1][j0]=0;
  goto loop;
}   /* End of Func*/


# include <stdio.h>
# include <math.h>

main()
{
	 
	float a[4][4] = {   {0, 0, 0, 0},
						{0, 3, 4, 2},
						{0, 2, 1, 2},
						{0, 1, 3, 2}	};
	float xb[4][2] = {  {0, 0},
						{0, 60},
						{0, 40},
						{0, 80}			};
	float c[2][4] = {   {0, 0, 0, 0},
						{0, -2, -4, -3} };
	float b[4][4],b1[4][4],e[4][4],ps[4][2],qs[4][2],cb[2][4],pa[2][4],ch[2][4];
	int k[4];
	xi_rsm(3, 3, a, b, b1, e, &xb, ps, qs, cb, pa, c, ch, &k);
} /* End of main */