//高斯若当主函数 
Matrix* Gauss_Jordan_Elimination(Matrix* mat, int eliminateMode)
{
	int eliminateModeBak =  eliminateMode;
	
	if(eliminateMode == ELIMINATE_INVERSE2)
	{
		eliminateMode = ELIMINATE_INVERSE;
	}

	if(eliminateMode == ELIMINATE_GET_DET)
	{
		eliminateMode = ELIMINATE_TRIANGLE;
	}
	
	//获取第一个矩阵列数 
	int n1 = ( *(mat + 0) );
	//获取第一个矩阵行数 
	int m1 = ( *(mat + 1) );


	//低于二阶的矩阵没法阶梯化 
	if(n1 <2 || m1 < 2)
	{
		PRINTF_ERROR("The 1-order matrix can't run this.\n", "");
		exit(1);
		return real_matrix(1.0f / (*(mat + 2)));
	}

	//用传入矩阵的备份做 高斯若当处理 
	Matrix* mat1 = (Matrix *)calloc(n1 * m1 + MAT_PRELEN, sizeof(double));
	memcpy(mat1, mat, sizeof(double) * (n1 * m1 + MAT_PRELEN) ); 

	//行交换存储器 
	double exChangeLineTmp[MAX_LEN];
	
	Matrix* matInverse = NULL; 
	
	if(eliminateMode == ELIMINATE_INVERSE)
	{
		//动态分配一个存放逆矩阵的空间 
		matInverse = (Matrix *)calloc(n1 * m1 + MAT_PRELEN, sizeof(double));
		//写入列行数 
		*(matInverse + 0) =  (double)n1;
		*(matInverse + 1) =  (double)m1;
		// 先置一同阶单位矩阵 
		for(int i = 0; i < GET_MIN(n1, m1); i++) 
		{
			*(matInverse + MAT_PRELEN + i + i * n1) = 1.0f;
		}
	}
	
	//行交换的次数 
	int exChangeTimes = 0;
	//高斯若当消去法置换成上三角型矩阵 
	int i = 0, j = 0, j_start = 0;
	for(i = 0; i < GET_MIN(n1, m1); i++) 
	{	
		for(j = j_start; (j < m1) && (i < n1); j++) 
		{
			//找到主元直接中断内循环 
			if( (*(mat1 + MAT_PRELEN  + i + j_start * n1)) != 0)
			{
				break;
			}
			//找到主元为0，则通过行交换使得主元不为0 
			if( (*(mat1 + MAT_PRELEN  + i + j*n1) != 0) && (j != j_start))
			{
				
				memcpy(exChangeLineTmp, mat1 + MAT_PRELEN + j*n1, n1 * sizeof(double));			
				// i行与j行互换 
				memcpy(mat1 + MAT_PRELEN + j * n1,    mat1 + MAT_PRELEN + i * n1,            n1 * sizeof(double));
				memcpy(mat1 + MAT_PRELEN + i * n1,    exChangeLineTmp,                       n1 * sizeof(double));
				
				// 如果是求逆矩阵模式 
				if(eliminateMode == ELIMINATE_INVERSE)
				{ 
					memcpy(exChangeLineTmp, matInverse + MAT_PRELEN + j*n1, n1 * sizeof(double));			
					// i行与j行互换 
					memcpy(matInverse + MAT_PRELEN + j * n1,    matInverse + MAT_PRELEN + i * n1,            n1 * sizeof(double));
					memcpy(matInverse + MAT_PRELEN + i * n1,    exChangeLineTmp,                             n1 * sizeof(double));
				}
						
				exChangeTimes ++;
				break;
			}
		}

		//如果该列主元都为0，直接跳过，处理下一列 
		if((j == m1) && (*(mat1 + MAT_PRELEN  + i + j_start * n1) == 0))
		{
			// 这里根据开关选择是否采用阶梯消去法，阶梯消去用于非方阵的m*n型矩阵秩计算 
			if((eliminateModeBak != ELIMINATE_LADDER))
			{
				j_start ++;
			}
			continue;
		}

		//如果该列主元不为0，则本行除以主元 
		if( (*(mat1 + MAT_PRELEN  + i + j_start * n1)) != 0)
		{
			// 这里根据开关选择是否采用主元归一，
			if((eliminateMode == ELIMINATE_INVERSE))
			{
				// 主元归一, 避免误差 
				double mainPara = (*(mat1 + MAT_PRELEN  + i + j_start*n1));
				for(int k = i; k < n1; k++) 
				{
					(*(mat1 + MAT_PRELEN  + k + j_start*n1)) /= mainPara;
					
					// 精度处理
					_ZERO_( *(mat1 + MAT_PRELEN  + k + j_start*n1) );
				}
				
				if(eliminateMode == ELIMINATE_INVERSE)
				{
					for(int k = 0; k < n1; k++) 
					{
						(*(matInverse + MAT_PRELEN  + k + j_start*n1)) /= mainPara;

						// 精度处理
						_ZERO_( *(matInverse + MAT_PRELEN  + k + j_start*n1) );	
					}
				}
			}
		}

		//开始消元 
		for(j = ( (eliminateMode==ELIMINATE_DIAG) || (eliminateMode==ELIMINATE_INVERSE) )?0 :(j_start + 1); (j < m1) && (i < n1); j++) 
		{
			// 主元不消去自己 
			if(j == j_start)
			{
				continue;
			}
			
			//计算消去系数
			double disLinePara =  - (*(mat1 + MAT_PRELEN  + i + j*n1)) / (*(mat1 + MAT_PRELEN  + i + j_start*n1));
				
			//用主元行消去其他行 
			for(int k = i; k < n1; k++) 
			{
				(*(mat1 + MAT_PRELEN  + k + j*n1))  +=  disLinePara * (*(mat1 + MAT_PRELEN  + k + j_start*n1));
					
				// 精度处理
				_ZERO_(*(mat1 + MAT_PRELEN  + k + j*n1));
			}
			
			// 如果是求逆矩阵模式 
			if(eliminateMode == ELIMINATE_INVERSE)
			{
				//计算消去系数 (这里沿用左边的消去系数，因为 AE -> EA^(-1),故两边乘以的消去系数一致
				double disLineParaInverse = disLinePara ;
					
				//用主元行消去其他行 
				for(int k = ((eliminateMode==ELIMINATE_DIAG) || (eliminateMode==ELIMINATE_INVERSE) ?0 :i); k < n1; k++) 
				{
					(*(matInverse + MAT_PRELEN  + k + j*n1))  += disLineParaInverse * (*(matInverse + MAT_PRELEN  + k + j_start*n1));
					// 精度控制 
					_ZERO_( *(matInverse + MAT_PRELEN  + k + j*n1) );
				}
			} 
		}

		//下沉一行
		j_start ++;	
	}
	
	// 如果是det求值模式 
	if(eliminateModeBak == ELIMINATE_GET_DET)
	{	
		double detValue = 1.0f;	
		for(int j = 0; j < m1; j++) 
		{
			 detValue  *=  (*(mat1 + MAT_PRELEN  + j + j*m1));
		}
		
		// 精度控制 
		_ZERO_(detValue);
				
		//释放备份矩阵内存 
		free(mat1);
			
		//返回det值 1阶矩阵 
		return  real_matrix(detValue);
	}
	
	//如果是求逆矩阵模式 
	if(eliminateModeBak == ELIMINATE_INVERSE)
	{
		//获取首非零行个数 
		for(int j = 0; j < n1; j++) 
		{
			double mPara = (*(mat1 + MAT_PRELEN  + j + j*n1));
			if( mPara == 0.0f )
			{
				PRINTF_ERROR("The matrix has no inverse mat.\n", "");
				exit(1);
			}
		}
		
		//释放备份矩阵内存 
		free(mat1);	
		return  matInverse;
	}

	//如果是解方程组模式 
	if(eliminateModeBak == ELIMINATE_INVERSE2)
	{
		free(matInverse);
		return mat1;
	}

	//返回行交换次数 
	return mat1;
} 
复制代码