线性方程组的迭代求解(C++)

1. 简单迭代法

　　就是将一个线性方程组Ax=b等价变形为x=Bx+g，由此建立迭代公式$x^{(k+1)}=Bx^{(k)}+g$，取初始向量$x^0$，进行迭代。编写程序需要写几个函数：矩阵乘向量，向量相加。以下程序打印矩阵B和向量g，并打印迭代过程中每一步的x值。

void printMatrix(double *matrix, int row, int column) {
    for (int i = 0; i < row; i++) {
        for (int j = 0; j < column; j++) {
            printf("%6.2lf    ", *(matrix + i * column + j));
        }
        printf("
");
    }
}

void printVector(double *matrix, int len) {
    for (int i = 0; i < len; i++) {
        printf("%6.2lf    ", *(matrix + i));
    }
    printf("
");
}
double* MatrixMultiVector(double* matrix, int row, int column, double* vector) {
    double* result = new double[row];
    for (int i = 0; i < row; i++) {
        *(result + i) = 0;
        for (int j = 0; j < column; j++) {
            *(result + i) -= *(matrix + i * row + j) * (*(vector + j));
        }
    }
    return result;
}

void vectorAdd(double* vector1, double* vector2, int len) {
    for (int i = 0; i < len; i++) {
        *(vector1 + i) += *(vector2 + i);
    }
}

void simpleIterativeMethod(double *matrix, double* vector, double* X0, int order) {
    double* X;

    double temp;
    for (int i = 0; i < order;i++) {
        temp = *(matrix + i * order + i);
        for (int j = 0; j < order; j++) {
            if (i == j) {
                continue;
            }
            *(matrix + i * order + j) /= temp;
        }
        *(matrix + i * order + i) = 0;
        *(vector + i) /= temp;
    }
    printMatrix(matrix, order, order);
    printVector(vector, 3);
    int MAX_ITER = 10;
    int ind = 1;
    X = MatrixMultiVector(matrix, order, order, X0);
    vectorAdd(X, vector, order);
    printf("----迭代----
");
    printVector(X, order);
    while (ind < MAX_ITER) {
        X0 = X;
        X = MatrixMultiVector(matrix, order, order, X0);
        printVector(X, 3);
        vectorAdd(X , vector, order);
        printf("iter %d: ", ind);
        printVector(X, order);
        ind++;
    }
    if (X) {
        delete[] X;
    }
}
int main() {
    double matrix[] = {
        10, -2,-1,
        -2,10,-1,
        -1,-2,5
    };
    double vector[] = { 3,15,10 };
    double X0[] = { 0,0,0 };
    simpleIterativeMethod(matrix, vector, X0, 3);
    return 0;
}