数值线性代数实验-共轭梯度法

一开始用c++的运算符重载程序总是莫名其妙的崩掉，然后以为是运算符重载的问题就写了个class对矩阵重新封装，结果还是崩，然后好久才发现是我把空间开的太大导致程序崩掉，无语，这样就浪费了我一个上午。。。。

课本上的例题：

$$
\begin{equation}
{\left[ \begin{array}{ccc}
10&1&2&3&4\\
1& 9& -1& 2& -3\\
2& -1& 7& 3& -5\\
3& 2& 3& 12& -1\\
4& -3& -5 &-1 &15
\end{array}
\right ]}\times {
\left[ \begin{array}{ccc}
x_1\\
x_2\\
x_3\\
x_4\\
x_5
\end{array}
\right ]}
={
\left[ \begin{array}{ccc}
12\\ -27\\ 14\\ -17\\ 12
\end{array}
\right ]}
\end{equation}
$$

输入数据

5 5
10 1 2 3 4
1 9 -1 2 -3
2 -1 7 3 -5
3 2 3 12 -1
4 -3 -5 -1 15
12 -27 14 -17 12

#include "cstdio"
#include "cstring"
#include "cstdlib"
#include "cmath"
#include "iostream"
using namespace std;

const double eps = 1e-6;
const int maxr = 10;
const int maxc = 10;

class Matrix {
private:
    double m[maxr][maxc];
    int row, col;
public:
    Matrix operator - (Matrix b) {
        Matrix c;
        c.col = this->col; c.row = this->row;
        for (int i = 0; i < this->row; i++) {
            for (int j = 0; j < this->col; j++) {
                c.m[i][j] = this->m[i][j]-b.m[i][j];
            }
        }
        return c;
    }
    Matrix operator + (Matrix b) {
        Matrix c;
        c.col = this->col; c.row = this->row;
        for (int i = 0; i < this->row; i++) {
            for (int j = 0; j < this->col; j++) {
                c.m[i][j] = this->m[i][j] + b.m[i][j];
            }
        }
        return c;
    }
    Matrix operator * (Matrix b) {
        Matrix c;
        memset(c.m, 0, sizeof(c.m));
        c.row = this->row; c.col = b.col;
        for (int i = 0; i < this->row; i++) {
            for (int k = 0; k < this->col; k++) {
                for (int j = 0; j < b.col; j++) {
                    c.m[i][j] += this->m[i][k]*b.m[k][j];
                }
            }
        }
        return c;
    }
    Matrix operator & (double a) {
        Matrix c = *this;
        for (int i = 0; i < this->row; i++) {
            for (int j = 0; j < this->col; j++) {
                c.m[i][j] = a*this->m[i][j];
            }
        }
        return c;
    }
    Matrix operator !() {
        Matrix c = *this;
        c.row = this->col, c.col = this->row;
        for (int i = 0; i < this->row; i++) {
            for (int j = 0; j < this->col; j++) {
                c.m[j][i] = this->m[i][j];
            }
        }
        return c;
    }
    int dcmp() {
        for (int i = 0; i < this->row; i++) {
            for (int j = 0; j < this->col; j++) {
                if (fabs(this->m[i][j]) > eps) return 1;
            }
        }
        return 0;
    }
    void show_Matrix() {
        for (int i = 0; i < this->row; i++) {
            printf("%.6f", this->m[i][0]);
            for (int j = 1; j < this->col; j++) {
                printf(" %.8f", this->m[i][j]);
            }
            printf("\n");
        }
    }
    void set_matrix(int n, int m) {
        this->row = n, this->col = m;
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                scanf("%lf", &this->m[i][j]);
            }
        }
    }
    void init(int n, int m) {
        this->row = n, this->col = m;
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                this->m[i][j] = 1.0;
            }
        }
    }
    int Col(){return this->col;}
    int Row(){return this->row;}
    double get(int i, int j) {return this->m[i-1][j-1];}
};

Matrix Conjugate_gradient(Matrix A, Matrix B) {
    Matrix x0 = B, p0, r1, p1;
    double alf, bet;
    x0.init(A.Col(), 1);
    Matrix r0 = B-(A*x0);
    int k = 0;
    r1 = r0;
    while (r1.dcmp()) {
        k++;
        if (k == 1) p1 = r0;
        else {
            bet = (!r1*r1).get(1, 1)/(!r0*r0).get(1, 1);
            p1 = r1 + (p0&bet);
        }
        alf = (!r1*r1).get(1, 1)/(!p1*A*p1).get(1, 1);
        x0 = x0 + (p1&alf);
        r0 = r1;
        r1 = r1-((A&alf)*p1);
    }
    return x0;
}

int main(int argc, char const *argv[])
{
   // freopen("in.txt", "r", stdin);
    Matrix A, B, C;
    int n, m;
    printf("输入系数矩阵的维数和矩阵A\n");
    scanf("%d %d", &n, &m);
    A.set_matrix(n, m);
    printf("输入矩阵B\n");
    B.set_matrix(n, 1);
    C = Conjugate_gradient(A, B);
    printf("利用共轭梯度法得到的解为\n");
    C.show_Matrix();
    return 0;
}