描述

本来是帮助同寝室的测绘专业的同学去完成“摄影测量”的“后方交会”作业，其中用到了很多线性代数的矩阵运算，于是我就写了这个cpp，心想也能帮他们减轻运算负担。

Matrix类

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
class Matrix{
    private :
        double Val[10][10];
        int Column,Row;        
    public :
        double Val[10][10];
        int Column,Row;
        void show(){
            puts("****");
            cout<<"Row = "<<Row<<" Column = "<<Column<<endl;
            puts("-----------------");
            for(int i = 1 ; i <= Row  ; i++){
                for(int j = 1 ; j <= Column ; j++){
                    printf("%.3lf\t%c",Val[i][j],j==Column?'\n':' ');
                }
            }
            puts("****");
        }
        void init(){
            cout<<"gets row and column first"<<endl;
            cin>>Row>>Column;
            cout<<"gets values now"<<endl;
            for(int i = 1 ; i <= Row  ; i++){
                for(int j = 1 ; j <= Column ; j++){
                   cin>>Val[i][j];
                }
            }
        }
        Matrix Transpose(){
            Matrix ret;
            ret.Column=Row;
            ret.Row=Column;
            for(int i=1 ; i<=Row ; i++){
                for(int j=1 ; j<=Column ; j++){
                    ret.Val[j][i]=Val[i][j];
                }
            }
            return ret;
        }
        Matrix Multiply(Matrix a){
            Matrix ret;
            ret.Row=Row;
            ret.Column=a.Column;
            for(int i=1 ; i<= ret.Row ; i++){
                for(int j=1 ; j<= ret.Column ; j++){
                    double tmp=0;
                    for(int k=1 ; k<= Column ; k++){
                        tmp+=Val[i][k]*a.Val[k][j];
                    }
                    ret.Val[i][j]=tmp;
                }
            }
            return ret;
        }
        //取绝对值必须是方阵，即行数等于列数
        double Absolute(){
            double ret=0;
            if(Row==2){
                ret=Val[1][1]*Val[2][2]-Val[1][2]*Val[2][1];
                return ret;
            }
            else{
                Matrix tmp[Row+1][Row+1];
                for(int i=1 ; i<= Row ; i++){
                    for(int j=1 ; j<= Column ; j++){
                        int pre;
                        if((i+j)%2==0)pre=1;
                        else pre=-1;
                        tmp[i][j].Row=Row-1;
                        tmp[i][j].Column=Column-1;
                        for(int ii=1,xi=1 ; ii <= Row ; ii++){
                            if(ii==i)continue;
                            for(int jj=1,xj=1 ; jj <= Column ; jj++){
                                 if(jj==j)continue;
                                 tmp[i][j].Val[xi][xj]=Val[ii][jj];
                                 xj++;
                            }
                            xi++;
                        }
//                        cout<<" i = "<<i<<" j = "<<j<<endl;
//                        tmp[i][j].show();
                        ret+=pre*Val[i][j]*tmp[i][j].Absolute();
                    }
                }
                return ret/Row;
            }
        }
        Matrix Adjoint_matrix (){
            Matrix ret;
            Matrix tmp[Row+1][Row+1];
            ret.Row=Row;
            ret.Column=Column;
            for(int i=1 ; i<= Row ; i++){
                for(int j=1 ; j<= Column ; j++){
                    double values;
                    int pre;
                    if((i+j)%2==0)pre=1;
                    else pre=-1;
                    tmp[i][j].Row=Row-1;
                    tmp[i][j].Column=Column-1;
                    for(int ii=1,xi=1 ; ii <= Row ; ii++){
                        if(ii==i)continue;
                        for(int jj=1,xj=1 ; jj <= Column ; jj++){
                             if(jj==j)continue;
                             tmp[i][j].Val[xi][xj]=Val[ii][jj];
                             xj++;
                        }
                        xi++;
                    }
//                  cout<<" i = "<<i<<" j = "<<j<<endl;
//                  tmp[i][j].show();
                    values=pre*tmp[i][j].Absolute();
                    ret.Val[i][j]=values;
                }
            }
            //正好反了，不知道为什么，以后再看
            return ret.Transpose();
        }
        Matrix Inverse (Matrix a){
            Matrix ret=a.Adjoint_matrix();
            double ab=a.Absolute();
            for(int i=1 ; i<= Row ; i++){
                for(int j=1 ; j<= Row ; j++){
                    ret.Val[i][j]/=ab;
                }
            }
            return ret;
        }
};

其中init初始化矩阵时，先输入行数，再输入列数。然后输入矩阵。

  • void show() 显示矩阵
  • void init() 初始化矩阵
  • Matrix Transpose() 矩阵转置
  • Matrix Multiply(Matrix a) 矩阵相乘
  • double Absolute() 矩阵(行列式)的绝对值
  • Matrix Adjoint_matrix () 求伴随矩阵
  • Matrix Inverse (Matrix a) 求逆矩阵 用法： b=a.Inverse(a)

出错

有些人用VS编译会发现不能成功，是因为

解决方法

用GCC的编译器，codeblocks，dev都能编译通过的。

后记

我的代码只适用于9阶一下的矩阵运算，由于代码的算法和完成方式比较低级，不适应高阶。
如果要高阶的话，要用高斯——但丁法，我还不会。……果然弱渣

鸣谢

在编写过程中得到了孔民同学的帮助。
代码完成后得到了邓仕宏同学的校验。
十分感谢这两位同学的付出。