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

/* 使用Gauss-Jordan消去法求n階實矩陣的逆矩陣 */ 
int ReversalOp(double *matrixa, int rank) 
{ 
	int *is,*js,i,j,k,l,u,v; 
	double d,p; 
	is=(int*)malloc(rank*sizeof(int)); 
	js=(int*)malloc(rank*sizeof(int)); 


	for(k=0;k<=rank-1;k++) 
	{ 
		//第一步，全選主元
		//（從第 k 行、第 k 列開始的右下角子陣中選取絕對值最大的元素，并記住次元素所在的行號和列號，在通過行交換和列交換將它交換到主元素位置上。）
		d=0.0; 
		for (i=k; i<=rank-1; i++) 
			for (j=k; j<=rank-1; j++) 
			{ 
				l=i*rank+j; 
				p=fabs(matrixa[l]); 
				if (p>d) 
				{ 
					d=p; 
					is[k]=i;
					js[k]=j;
				} 
			} 
			if(d+1.0==1.0) 
			{
				free(is); 
				free(js); 
				printf("err**not inv\n"); 
				return(0); 
			} 
			if(is[k]!=k) 
				for(j=0;j<=rank-1;j++) 
				{
					u=k*rank+j;
					v=is[k]*rank+j; 
					p=matrixa[u]; 
					matrixa[u]=matrixa[v];
					matrixa[v]=p; 
				} 
				if(js[k]!=k) 
					for(i=0; i<=rank-1; i++) 
					{
						u=i*rank+k;
						v=i*rank+js[k]; 
						p=matrixa[u];
						matrixa[u]=matrixa[v];
						matrixa[v]=p; 
					} 
					
					//機算逆矩陣

					//第二步（m(k, k) = 1 / m(k, k)）
					l=k*rank+k; 
					matrixa[l]=1.0/matrixa[l]; 

					//第三步（m(k, j) = m(k, j) * m(k, k)，j = 0, 1, ..., n-1；j != k ）
					for(j=0;j<=rank-1;j++) 
						if (j!=k) 
						{
							u=k*rank+j; 
							matrixa[u]=matrixa[u]*matrixa[l];
						} 

						//第四步（m(i, j) = m(i, j) - m(i, k) * m(k, j)，i, j = 0, 1, ..., n-1；i, j != k ）
						for(i=0;i<=rank-1;i++) 
							if(i!=k) 
								for(j=0;j<=rank-1;j++) 
									if(j!=k) 
									{ 
										u=i*rank+j; 
										matrixa[u]=matrixa[u]-matrixa[i*rank+k]*matrixa[k*rank+j]; 
									} 
									//第五步（m(i, k) = -m(i, k) * m(k, k)，i = 0, 1, ..., n-1；i != k ）
									for(i=0;i<=rank-1;i++) 
										if(i!=k) 
										{
											u=i*rank+k; 
											matrixa[u]=-matrixa[u]*matrixa[l];
										} 
	} 
	/*最后，根據在全選主元過程中所記錄的行、列交換的信息進行恢復，恢復的原則如下：在全選主元過程中，先交換的
	行（列）后進行恢復；原來的行（列）交換用列（行）交換來恢復。*/
	for(k=rank-1;k>=0;k--) 
	{ 
		if (js[k]!=k) 
			for(j=0;j<=rank-1;j++) 
			{ 
				u=k*rank+j; 
				v=js[k]*rank+j; 
				p=matrixa[u];
				matrixa[u]=matrixa[v];
				matrixa[v]=p; 
			} 
			if(is[k]!=k) 
				for(i=0;i<=rank-1;i++) 
				{
					u=i*rank+k;
					v=i*rank+is[k]; 
					p=matrixa[u]; 
					matrixa[u]=matrixa[v]; 
					matrixa[v]=p; 
				} 
	} 
	free(is);
	free(js); 
	return(1); 
}