Young87

当前位置:首页 >个人收藏

矩阵求解(列行列主元消元法)

计算扩展矩阵解,利用列出行列的主元及利用双精度浮点数提高计算精度:

//Three order Mat
typedef struct STRTHREEORDERMAT
{
	//parameter and y (result)
	double m_param[3][4];
	//x (unknown)
	double m_unk[3];
}thrmat_t, *thrmat_tt;

//Three order Mat
typedef struct STRFOURORDERMAT
{
	//parameter and y (result)
	double m_param[4][5];
	//x (unknown)
	double m_unk[4];
}format_t, *format_tt;

//Six order Mat
typedef struct STRSIXORDERMAT
{
	//parameter and y (result)
	double m_param[6][7];
	//x (unknown)
	double m_unk[6];
}sixmat_t, *sixmat_tt;

//EliMination with maximal col and row privoting
static char gelimsolve(double **prow, short *pcol, double *unk, short row)
{
	double *temprow = NULL;
	short tempcol = 0;
	double calctemp = 0;
	short it0 = 0, it1 = 0, it2 = 0;
	short maxrow = 0, maxcol = 0;
	//elimination calc
	for (it0 = 0; it0 < row; ++it0)
	{
		//calc max col and row
		calctemp = FABS(*(*(prow + it0) + it0));
		maxrow = maxcol = it0;
		for (it1 = it0; it1 < row; ++it1)
		{
			for (it2 = it0; it2 < row; ++it2)
			{
				if (FABS(*(*(prow + it1) + *(pcol + it2))) > calctemp)
				{
					calctemp = FABS(*(*(prow + it1) + *(pcol + it2)));
					maxrow = it1;
					maxcol = it2;
				}
			}
		}
		//check if there no slove
		if (*(*(prow + maxrow) + *(pcol + maxcol)) == 0.0)
			return SONOSOLUT;
		else
		{
			//exchange row
			if (maxrow != it0)
			{
				temprow = *(prow + maxrow);
				*(prow + maxrow) = *(prow + it0);
				*(prow + it0) = temprow;
			}
			//exchange col
			if (maxcol != it0)
			{
				tempcol = *(pcol + maxcol);
				*(pcol + maxcol) = *(pcol + it0);
				*(pcol + it0) = tempcol;
			}
		}
		//calc current row calc col
		for (it1 = 0; it1 < it0; ++it1)
		{
			//check current cols value is not zero
			if (*(*(prow + it0) + *(pcol + it1)) != 0.0)
			{
				calctemp = *(*(prow + it1) + *(pcol + it1)) / *(*(prow + it0) + *(pcol + it1));
				for (it2 = it1 + 1; it2 < row + 1; ++it2)
				{
					*(*(prow + it0) + *(pcol + it2)) = 
						*(*(prow + it0) + *(pcol + it2)) * 
						calctemp - 
						*(*(prow + it1) + *(pcol + it2));
				}
				*(*(prow + it0) + *(pcol + it1)) = 0.0;
			}
		}
		//check curr cols is zro
		if (*(*(prow + it0) + *(pcol + it0)) == 0.0)
			--it0;
	}
	//return and calc unk
	for (it0 = row - 1; it0 > -1; --it0)
	{
		for (it1 = row - 1; it1 > it0; --it1)
			*(*(prow + it0) + row) -= *(*(prow + it0) + *(pcol + it1)) * *(unk + *(pcol + it1));
		*(unk + *(pcol + it0)) = *(*(prow + it0) + row) / *(*(prow + it0) + *(pcol + it0));
	}
	return SOSUCCESS;
}

//Gauss elimination solve 3 matrix
_EX_SOLE_ char gaelimsolve3mat(thrmat_tt mat)
{
	short it = 0;
	double *prow[3] = { 0 };
	short pcol[4] = { 0 };
	//init rows link
	for (it = 0; it < 3; ++it)
		*(prow + it) = *(mat->m_param + it);
	//init cols link
	for (it = 0; it < 4; ++it)
		*(pcol + it) = it;
	return gelimsolve(prow, pcol, mat->m_unk, 3);
}

//Gauss elimination solve 4 matrix
_EX_SOLE_ char gaelimsolve4mat(format_tt mat)
{
	short it = 0;
	double *prow[4] = { 0 };
	short pcol[5] = { 0 };
	//init rows link
	for (it = 0; it < 4; ++it)
		*(prow + it) = *(mat->m_param + it);
	//init cols link
	for (it = 0; it < 5; ++it)
		*(pcol + it) = it;
	return gelimsolve(prow, pcol, mat->m_unk, 4);
}

//Gauss elimination solve 6 matrix
_EX_SOLE_ char gaelimsolve6mat(sixmat_tt mat)
{
	short it = 0;
	double *prow[6] = { 0 };
	short pcol[7] = { 0 };
	//init rows link
	for (it = 0; it < 6; ++it)
		*(prow + it) = *(mat->m_param + it);
	//init cols link
	for (it = 0; it < 7; ++it)
		*(pcol + it) = it;
	return gelimsolve(prow, pcol, mat->m_unk, 6);
}

 

除特别声明,本站所有文章均为原创,如需转载请以超级链接形式注明出处:SmartCat's Blog

上一篇: python3的安装与环境变量的配置

下一篇: TypeScript开发环境搭建

精华推荐