矩阵的基本运算 ------ C 实现

最近写了一个C程序 ------------- 矩阵的基本运算。

包括四则运算， 求逆， rref（行最简型）， 转置。

并有几个实际应用例子， 求解方程组， 线性回归， 线性规划正在进行中。 

贴出来请大牛们指点指点。就贴最难的一个， rref吧。

/************************************
 * pTag is the pivots tag
 * if 1, the column is pivot
 * and 0, the column is not pivot
 * Return:
 *	return elimination matrix
*************************************/
matrix
rref(matrix A, int *pTag)         // reduced row echelon form
{
    // initialize Elimination matrix
    matrix E = Imat(A->rows);
    double **MA = A->data;
    double **ME = E->data;
    double e[A->rows];		// factor of reduced

    // echelon loop
    int k;
    for (int pivots = 0, j = 0; j < A->cols && pivots < A->rows; j++)
    {
	if (MA[j][pivots] == 0.0)
	{
	    // if this element is 0.0, find a nonzero element behind in j's column
	    for (k = pivots + 1; k < A->rows && MA[j][k] == 0.0; k++);

	    // row swap if find else operate next line
	    if (k != A->rows)
	    {
		RowSwap(A, pivots, k);
		RowSwap(E, pivots, k);
	    }
	    else
	    {
		pTag[j] = 0;
		continue;
	    }
	}

	// eliminate
	pTag[j] = 1;
	for (int i = 0; i < A->rows; i++)
	    e[i] = (i == pivots ? 0 : MA[j][i] / MA[j][pivots]);
	for (int col = 0; col < A->cols; col++)
	    for (int row = 0; row < A->rows; row++)
		MA[col][row] -= e[row] * MA[col][pivots];
	for (int col = 0; col < E->cols; col++)
	    for (int row = 0; row < E->rows; row++)
		ME[col][row] -= e[row] * ME[col][pivots];
	pivots++;
    }
    
    // normalise after echelon, maybe can layout it to 
    // echelon loop, but it will reduce the precision
    double fac;
    for (int i = 0, k = 0; i < A->rows; i++)
    {
	while (pTag[k++] == 0);
	fac = MA[k - 1][i];
	for (int j = 0; j < E->cols; j++)
	    ME[j][i] /= fac;
	for (int j = 0; j < A->cols; j++)
	    MA[j][i] /= fac;
    }
    return E;
}

上传半天不见有反应, CSDN, 有点怪. 完整历程及示例请到这里下载

这里