矩阵基本运算的实现(standard C++Version)
一年前 自己动手编写的代码，尽量使用了标准C++，其中用到了一些标准模板库STL，但没有用到很高级的部分，不懂STL的朋友也应该可以看懂，建议看一下《Generic Programming and STL》，有候捷译本《泛型编程与STL》。
头文件
实现文件
实现了很多基本的操作，需要的话还可以自己扩充。如果实现那些注释掉的代码，功能就比较强了。
另有标准C语言版本。
头文件
1
/***
2
*
3* author:XieXiaokui
5* purpose:Defines functions for matrix
6*
7
***/
8
#pragma once
9
10#include <iostream>
11#include <fstream>
12
13#include <string>
14#include <sstream>
15#include <algorithm>
16#include <functional>
17#include <numeric>
18#include <iterator>
19#include <cassert>
20
21#include "MatrixException.h"
22
23using namespace std;
24
25class Matrix
26{
27public:
28
29 // Constructors
30 explicit Matrix();
31 explicit Matrix(int size);
32 Matrix(int row,int col);
33 Matrix(const Matrix& m);
34
35 // Destructor
36 ~Matrix();
37
38 // Assignment operators
39 Matrix& operator= (const Matrix& m);
40
41 // Value extraction method
42 int GetRow() const;
43 int GetCol() const;
44
45 // Subscript operator
46 double operator()(int i,int j)const; //subscript operator to get individual elements
47 double& operator()(int i,int j); //subscript operator to set individual elements
48
49 // Unary operators
50 Matrix operator+() const; //unary negation operator
51 Matrix operator-() const; //unary negation operator
52
53 //Binary operators
54 Matrix operator+(const Matrix& m) const;
55 Matrix operator-(const Matrix& m) const;
56 Matrix operator*(const Matrix& m) const;
57// Matrix operator*(double d)const;
58 Matrix operator/(const Matrix& m) const;
59 Matrix operator/(double d) const;
60 Matrix operator^(int pow) const;
61
62 bool operator== (const Matrix& m) const; //logical equal-to operator
63 bool operator!= (const Matrix& m) const; //logical not-equal-to operator
64
65 friend Matrix operator* (double d,const Matrix& m);
66 friend Matrix operator/ (double d,const Matrix& m);
67
68
69
70 // Combined assignment - calculation operators
71 Matrix& operator +=(const Matrix& m) const;
72 Matrix& operator -=(const Matrix& m) const;
73 Matrix& operator *=(const Matrix& m) const;
74 Matrix& operator *=(double d) const;
75 Matrix& operator /=(const Matrix& m) const;
76 Matrix& operator /=(double d) const;
77 Matrix& operator ^=(int pow) const;
78
79 // Miscellaneous -methods
80 void SetZero() ; //zero matrix:零阵
81 void SetUnit() ; //unit matrix:单位阵
82 void SetSize(int size) ; //rsizing matrix
83 void SetSize(int row,int col) ; //resizing matrix
84
85 // Utility methods
86 Matrix Solve(const Matrix& m)const; //
87 Matrix Adjoin() const; //adjoin matrix:伴随矩阵
88 double Determinant() const; //determinant:行列式
89 double Norm() const; //norm:模
90 Matrix Inverse() const; //inverse:逆阵
91 Matrix Transpose() const; //transpose:转置
92 double Cofactor() const; //confactor
93 double Condition() const; //the condition number of a matrix
94 int Pivot(int row) const; // partial pivoting
95
96 //primary change
97 Matrix& Exchange(int i,int j);// 初等变换 对调两行：ri<-->rj
98 Matrix& Multiple(int index,double k); //初等变换 第index 行乘以k
99 Matrix& MultipleAdd(int index,int src,double k); //初等变换 第src行乘以k加到第index行
100
101
102 // Type of matrices
103 bool IsSquare() const; //determine if the matrix is square:方阵
104 bool IsSingular() const; //determine if the matrix is singular奇异阵
105 bool IsDiagonal() const; //determine if the matrix is diagonal对角阵
106 bool IsScalar() const; //determine if the matrix is scalar数量阵
107 bool IsUnit() const; //determine if the matrix is unit单位阵
108 bool IsZero() const; //determine if the matrix is zero零矩阵
109 bool IsSymmetric() const; //determine if the matrix is symmetric对称阵
110 bool IsSkewSymmetric() const; //determine if the matrix is skew-symmetric斜对称阵
111 bool IsUpperTriangular() const; //determine if the matrix is upper-triangular上三角阵
112 bool IsLowerTriangular() const; //determine if the matrix is lower-triangular下三角阵
113
114 // output stream function
115 friend ostream& operator<<(ostream& os,const Matrix& m);
116
117 // input stream function
118 friend istream& operator>>(istream& is,Matrix& m);
119
120 //conert to string
121 string ToString() const;
122
123protected:
124
125 //delete the matrix
126 void Create(int row,int col);
127 void Clear();
128
129private:
130
131 const static double epsilon;
132
133 double** m_data;
134 size_t m_row;
135 size_t m_col;
136
137};
/***
2
*3
* author:XieXiaokui5
* purpose:Defines functions for matrix6
*7
***/8
#pragma once9
10
#include <iostream>11
#include <fstream>12
13
#include <string>14
#include <sstream>15
#include <algorithm>16
#include <functional>17
#include <numeric>18
#include <iterator>19
#include <cassert>20
21
#include "MatrixException.h"22
23
using namespace std;24
25
class Matrix26
{
27
public:28
29
// Constructors30
explicit Matrix();31
explicit Matrix(int size);32
Matrix(int row,int col);33
Matrix(const Matrix& m);34
35
// Destructor36
~Matrix();37
38
// Assignment operators39
Matrix& operator= (const Matrix& m);40
41
// Value extraction method42
int GetRow() const;43
int GetCol() const;44
45
// Subscript operator46
double operator()(int i,int j)const; //subscript operator to get individual elements47
double& operator()(int i,int j); //subscript operator to set individual elements48
49
// Unary operators50
Matrix operator+() const; //unary negation operator51
Matrix operator-() const; //unary negation operator52
53
//Binary operators54
Matrix operator+(const Matrix& m) const;55
Matrix operator-(const Matrix& m) const;56
Matrix operator*(const Matrix& m) const;57
// Matrix operator*(double d)const;58
Matrix operator/(const Matrix& m) const;59
Matrix operator/(double d) const;60
Matrix operator^(int pow) const;61
62
bool operator== (const Matrix& m) const; //logical equal-to operator63
bool operator!= (const Matrix& m) const; //logical not-equal-to operator64
65
friend Matrix operator* (double d,const Matrix& m);66
friend Matrix operator/ (double d,const Matrix& m);67
68
69
70
// Combined assignment - calculation operators71
Matrix& operator +=(const Matrix& m) const;72
Matrix& operator -=(const Matrix& m) const;73
Matrix& operator *=(const Matrix& m) const;74
Matrix& operator *=(double d) const;75
Matrix& operator /=(const Matrix& m) const;76
Matrix& operator /=(double d) const;77
Matrix& operator ^=(int pow) const;78
79
// Miscellaneous -methods80
void SetZero() ; //zero matrix:零阵81
void SetUnit() ; //unit matrix:单位阵82
void SetSize(int size) ; //rsizing matrix83
void SetSize(int row,int col) ; //resizing matrix84
85
// Utility methods86
Matrix Solve(const Matrix& m)const; //87
Matrix Adjoin() const; //adjoin matrix:伴随矩阵88
double Determinant() const; //determinant:行列式89
double Norm() const; //norm:模90
Matrix Inverse() const; //inverse:逆阵91
Matrix Transpose() const; //transpose:转置92
double Cofactor() const; //confactor93
double Condition() const; //the condition number of a matrix94
int Pivot(int row) const; // partial pivoting95
96
//primary change97
Matrix& Exchange(int i,int j);// 初等变换 对调两行：ri<-->rj98
Matrix& Multiple(int index,double k); //初等变换 第index 行乘以k99
Matrix& MultipleAdd(int index,int src,double k); //初等变换 第src行乘以k加到第index行100
101
102
// Type of matrices103
bool IsSquare() const; //determine if the matrix is square:方阵104
bool IsSingular() const; //determine if the matrix is singular奇异阵105
bool IsDiagonal() const; //determine if the matrix is diagonal对角阵106
bool IsScalar() const; //determine if the matrix is scalar数量阵107
bool IsUnit() const; //determine if the matrix is unit单位阵108
bool IsZero() const; //determine if the matrix is zero零矩阵109
bool IsSymmetric() const; //determine if the matrix is symmetric对称阵110
bool IsSkewSymmetric() const; //determine if the matrix is skew-symmetric斜对称阵111
bool IsUpperTriangular() const; //determine if the matrix is upper-triangular上三角阵112
bool IsLowerTriangular() const; //determine if the matrix is lower-triangular下三角阵113
114
// output stream function115
friend ostream& operator<<(ostream& os,const Matrix& m);116
117
// input stream function118
friend istream& operator>>(istream& is,Matrix& m);119
120
//conert to string121
string ToString() const;122
123
protected:124
125
//delete the matrix126
void Create(int row,int col);127
void Clear();128
129
private:130
131
const static double epsilon;132
133
double** m_data;134
size_t m_row;135
size_t m_col;136
137
}; 1/***
2*
3*
4* author:XieXiaokui
5* purpose:Defines functions for matrix
6*
7***/
8#include "stdafx.h"
9
10#include <iostream>
11#include <fstream>
12
13#include <string>
14#include <sstream>
15#include <algorithm>
16#include <functional>
17#include <numeric>
18#include <iterator>
19#include <cmath>
20#include <cassert>
21
22#include "matrix.h"
23
24using namespace std;
25
26const double Matrix::epsilon = 1e-7;
27
28// constructor
29
30Matrix::Matrix():m_row(0),m_col(0),m_data(0)
31{
32}
33
34
35// constructor
36Matrix::Matrix(int size):m_row(size),m_col(size)
37{
38 if(size<=0)
39
{
40 m_row=0;
41 m_col=0;
42 m_data=0;
43
44 ostringstream oss;
45 oss<<"In Matrix::Matrix(int size) size "<<size<<" <=0.Please check it.";
46 throw oss.str();
47
48
}
49
50 m_data=new double*[size];
51
52 for(int i=0;i<size;i++)
53 {
54 m_data[i]=new double[size];
55 }
56
57}
58
59// constructor
60Matrix::Matrix(int row,int col):m_row(row),m_col(col)
61{
62 if(row<=0)
63 {
64 m_row=0;
65 m_col=0;
66 m_data=0;
67
68 ostringstream oss;
69 oss<<"In Matrix::Matrix(int row,int col),row "<<row<<" <=0.Please check it.";
70 throw oss.str();
71 }
72 if(col<=0)
73 {
74 m_row=0;
75 m_col=0;
76 m_data=0;
77
78 ostringstream oss;
79 oss<<" In Matrix::Matrix(int row,int col),col "<<col<<" <=0.Pleasecheck it.";
80 throw oss.str();
81 }
82
83 m_data=new double*[row];
84 for(int i=0;i<row;i++)
85 {
86 m_data[i]=new double[col];
87 }
88}
89
90//copy constructor
91Matrix::Matrix(const Matrix& m):m_row(m.m_row),m_col(m.m_col)
92{
93 m_data = new double*[m_row];
94 for(int i=0;i<m_row;i++)
95 {
96 m_data[i] = new double[m_col];
97 }
98
99 for(int i=0;i<m_row;i++)
100 {
101 copy(m.m_data[i],m.m_data[i]+m_col,m_data[i]);
102 }
103
104}
105
106Matrix::~Matrix()
107{
108 if(m_data != 0)
109 {
110 for(int i=0;i<m_row;i++)
111 {
112 delete[] m_data[i];
113 }
114 delete[] m_data;
115 }
116}
117
118void Matrix::SetZero()
119{
120 for(int i=0;i<m_row;i++)
121 fill(m_data[i],m_data[i]+m_col,0);
122}
123
124void Matrix::SetUnit()
125{
126 for(int i=0;i<m_row;i++)
127 for(int j=0;j<m_col;j++)
128 m_data[i][j] = (i==j)?1:0;
129}
130
131void Matrix::SetSize(int size)
132{
133 Clear();
134 Create(size,size);
135}
136
137
138void Matrix::SetSize(int row,int col)
139{
140 Clear();
141 Create(row,col);
142}
143
144
145void Matrix::Clear()
146{
147
148 if(m_data != 0)
149 {
150 for(int i=0;i<m_row;i++)
151 {
152 delete[] m_data[i];
153 }
154 delete[] m_data;
155 }
156
157 m_row=0;
158 m_col=0;
159 m_data=0;
160}
161
162/*
163Matrix& Matrix::Create(int size)
164{
165 if(size<=0)
166 {
167 ostringstream oss;
168 oss<<"In Matrix::Create(int size),size "<<size<<" <=0.Please check it.";
169
170 throw oss.str();
171 }
172
173 if(m_data != 0)
174 {
175 for(int i=0;i<m_row;i++)
176 {
177 delete[] m_data[i];
178 }
179 delete[] m_data;
180 }
181
182 m_row=size;
183 m_col=size;
184
185 m_data=new double*[size];
186 for(int i=0;i<size;i++)
187 {
188 m_data[i]=new double[size];
189 }
190
191 for(int i=0;i<size;i++)
192 {
193 for(int j=0;j<size;j++)
194 {
195 m_data[i][j] = ((i==j)?1:0);
196 }
197 }
198
199 return *this;
200}
201*/
202
203void Matrix::Create(int row,int col)
204{
205 if(row<=0)
206 {
207 ostringstream oss;
208 oss<<"In Matrix::Create(int row,int col),row "<<row<<" <=0.Please check it.";
209 throw oss.str();
210 }
211 if(col<=0)
212 {
213 ostringstream oss;
214 oss<<"In Matrix::Create(int row,int col),col "<<col<<" <=0.Please check it.";
215 throw oss.str();
216 }
217
218 if(m_data != 0)
219 {
220 for(int i=0;i<m_row;i++)
221 {
222 delete[] m_data[i];
223 }
224 delete[] m_data;
225 }
226
227 m_row=row;
228 m_col=col;
229 m_data=new double*[row];
230 for(int i=0;i<row;i++)
231 {
232 m_data[i]=new double[col];
233 }
234}
235
236
237int Matrix::GetRow() const
238{
239 return m_row;
240}
241
242int Matrix::GetCol() const
243{
244 return m_col;
245}
246
247//transpose 转置
248Matrix Matrix::Transpose() const
249{
250 Matrix ret(m_col,m_row);
251 for(int i=0;i<m_row;i++)
252 {
253 for(int j=0;j<m_col;j++)
254 {
255 ret.m_data[j][i] = m_data[i][j];
256 }
257 }
258 return ret;
259}
260
261int Matrix::Pivot(int row) const
262{
263 int index=row;
264
265 for(int i=row+1;i<m_row;i++)
266 {
267 if(m_data[i][row] > m_data[index][row])
268 index=i;
269 }
270
271 return index;
272}
273
274
275Matrix& Matrix::Exchange(int i,int j) // 初等变换：对调两行ri<-->rj
276{
277 if(i<0 || i>=m_row)
278 {
279 ostringstream oss;
280 oss<<"In void Matrix::Exchange(int i,int j) ,i "<<i<<" out of bounds.Please check it.";
281 throw oss.str();
282 }
283
284 if(j<0 || j>=m_row)
285 {
286 ostringstream oss;
287 oss<<"In void Matrix::Exchange(int i,int j) ,j "<<j<<" out of bounds.Please check it.";
288 throw oss.str();
289 }
290
291 for(int k=0;k<m_col;k++)
292 {
293 swap(m_data[i][k],m_data[j][k]);
294 }
295
296 return *this;
297}
298
299
300Matrix& Matrix::Multiple(int index,double mul) //初等变换 第index 行乘以mul
301{
302 if(index <0 || index >= m_row)
303 {
304 ostringstream oss;
305 oss<<"In void Matrix::Multiple(int index,double k) index "<<index<<" out of bounds.Please check it.";
306 throw oss.str();
307 }
308
309 transform(m_data[index],m_data[index]+m_col,m_data[index],bind2nd(multiplies<double>(),mul));
310 return *this;
311}
312
313
314Matrix& Matrix::MultipleAdd(int index,int src,double mul) //初等变换 第src行乘以mul加到第index行
315{
316 if(index <0 || index >= m_row)
317 {
318 ostringstream oss;
319 oss<<"In void MultipleAdd(int index,int src,double k) index "<<index<<" out of bounds.Please check it.";
320 throw oss.str();
321 }
322
323 if(src < 0 || src >= m_row)
324 {
325 ostringstream oss;
326 oss<<"In void MultipleAdd(int index,int src,double k) src "<<src<<" out of bounds.Please check it.";
327 throw oss.str();
328 }
329
330 for(int j=0;j<m_col;j++)
331 {
332 m_data[index][j] += m_data[src][j] * mul;
333 }
334
335 return *this;
336
337}
338
339
340//inverse 逆阵：使用矩阵的初等变换，列主元素消去法
341Matrix Matrix::Inverse() const
342{
343 if(m_row != m_col) //非方阵
344 {
345 ostringstream oss;
346 oss<<"In Matrix Matrix::Invert() const. m_row "<<m_row<<" != m_col "<<m_col<<" Please check it";
347 throw oss.str();
348 }
349
350 Matrix tmp(*this);
351 Matrix ret(m_row); //单位阵
352 ret.SetUnit();
353
354 int maxIndex;
355 double dMul;
356
357 for(int i=0;i<m_row;i++)
358 {
359
360 maxIndex = tmp.Pivot(i);
361
362 if(tmp.m_data[maxIndex][i]==0)
363 {
364 ostringstream oss;
365 oss<<"In Matrix Matrix::Invert() const 行列式的值等于0.Please check it";
366 throw oss.str();
367 }
368
369 if(maxIndex != i) //下三角阵中此列的最大值不在当前行，交换
370 {
371 tmp.Exchange(i,maxIndex);
372 ret.Exchange(i,maxIndex);
373
374 }
375
376 ret.Multiple(i,1/tmp.m_data[i][i]);
377 tmp.Multiple(i,1/tmp.m_data[i][i]);
378
379
380 for(int j=i+1;j<m_row;j++)
381 {
382 dMul = -tmp.m_data[j][i];
383 tmp.MultipleAdd(j,i,dMul);
384 ret.MultipleAdd(j,i,dMul);
385
386 }
387
388 }//end for
389
390 for(int i=m_row-1;i>0;i--)
391 {
392 for(int j=i-1;j>=0;j--)
393 {
394 dMul = -tmp.m_data[j][i];
395 tmp.MultipleAdd(j,i,dMul);
396 ret.MultipleAdd(j,i,dMul);
397 }
398 }//end for
399
400 return ret;
401}
402
403
404
405// assignment operator 賦值运算符
406Matrix& Matrix::operator= (const Matrix& m)
407{
408 if(m_data != 0)
409 {
410 for(int i=0;i<m_row;i++)
411 {
412 delete[] m_data[i];
413 }
414 delete[] m_data;
415 }
416
417 m_row = m.m_row;
418 m_col = m.m_col;
419
420 m_data = new double*[m_row];
421 for(int i=0;i< m_row;i++)
422 {
423 m_data[i] = new double[m_col];
424 }
425
426 for(int i=0;i<m_row;i++)
427 {
428 copy(m.m_data[i],m.m_data[i]+m_col,m_data[i]);
429 }
430 return *this;
431}
432
433
434//binary addition 矩阵加
435Matrix Matrix::operator+ (const Matrix& m) const
436{
437 if(m_row != m.m_row)
438 {
439 ostringstream oss;
440 oss<<"In Matrix::operator+(const Matrix& m) this->m_row "<<m_row<<" != m.m_row "<<m.m_row<<" Please check it.";
441 throw oss.str();
442 }
443 if(m_col != m.m_col)
444 {
445 ostringstream oss;
446 oss<<" Matrix::operator+ (const Matrix& m) this->m_col "<<m_col<<" != m.m_col "<<m.m_col<<" Please check it.";
447 throw oss.str();
448 }
449
450 Matrix ret(m_row,m_col);
451 for(int i=0;i<m_row;i++)
452 {
453 transform(m_data[i],m_data[i]+m_col,m.m_data[i],ret.m_data[i],plus<double>());
454 }
455 return ret;
456}
457
458//unary addition 求正
459Matrix Matrix::operator+() const
460{
461 return *this;
462}
463
464//binary subtraction 矩阵减
465Matrix Matrix::operator- (const Matrix& m) const
466{
467 if(m_row != m.m_row)
468 {
469 ostringstream oss;
470 oss<<"In Matrix::operator-(const Matrix& m) this->m_row "<<m_row<<" != m.m_row "<<m.m_row<<" Please check it.";
471 throw oss.str();
472 }
473 if(m_col != m.m_col)
474 {
475 ostringstream oss;
476 oss<<"In Matrix::operator- (const Matrix& m) this->m_col "<<m_col<<" != m.m_col "<<m.m_col<<" Please check it.";
477 throw oss.str();
478 }
479
480 Matrix ret(m_row,m_col);
481 for(int i=0;i<m_row;i++)
482 {
483 transform(m_data[i],m_data[i]+m_col,m.m_data[i],ret.m_data[i],minus<double>());
484 }
485 return ret;
486}
487
488//unary substraction 求负
489Matrix Matrix::operator-() const
490{
491 Matrix ret(*this);
492
493 for(int i=0;i<m_row;i++)
494 for(int j=0;j<m_col;j++)
495 ret.m_data[i][j] *= -1;
496
497 return ret;
498
499}
500
501
502//binary multiple 矩阵乘
503Matrix Matrix::operator*(const Matrix& m) const
504{
505 if(m_col != m.m_row)
506 {
507 ostringstream oss;
508 oss<<"In Matrix::operator*(const Matrix& m) this->m_col "<<m_col<<" != m.m_row "<<m.m_row<<" Please check it.";
509 throw oss.str();
510 }
511
512 Matrix ret(m_row,m.m_col);
513 Matrix tmp=m.Transpose();
514
515 for(int i=0;i<m_row;i++)
516 {
517 for(int j=0;j<m.m_col;j++)
518 {
519 ret.m_data[i][j]=inner_product(m_data[i],m_data[i]+m_col,tmp.m_data[j],0.0);
520 }
521 }
522
523 return ret;;
524}
525
526//friend scalar multiple 数乘
527Matrix operator* (double d,const Matrix& m)
528{
529 Matrix ret(m);
530
531 for(int i=0;i<ret.m_row;i++)
532 for(int j=0;j<ret.m_col;j++)
533 ret.m_data[i][j] *= d;
534
535 return ret;
536
537}
538
539//binary matrix division equivalent to multiple inverse矩阵除，等价于乘以逆阵
540Matrix Matrix::operator/(const Matrix& m) const
541{
542 return *this * m.Inverse();
543}
544
545//binary scalar division equivalent to multiple reciprocal数除，等价于乘以此数的倒数
546Matrix Matrix::operator/(double d) const
547{
548 return 1/d * (*this);
549}
550
551//friend division
552Matrix operator/(double d,const Matrix& m)
553{
554 return d * m.Inverse();
555}
556
557
558// subscript operator to get individual elements 下标运算符
559double Matrix::operator()(int i,int j) const
560{
561 if(i<0 || i>= m_row)
562 {
563 ostringstream oss;
564 oss<<"In double Matrix::operator()(int i,int j) const.i "<<i<<" out of bounds.Please check it";
565 throw oss.str();
566 }
567 if(j<0 || j>= m_col)
568 {
569 ostringstream oss;
570 oss<<"In double Matrix::operator()(int i,int j) const.j "<<j<<" out of bounds.Please check it";
571 throw oss.str();
572 }
573
574 return m_data[i][j];
575}
576
577// subscript operator to set individual elements 下标运算符
578double& Matrix::operator ()(int i,int j)
579{
580 if(i<0 || i>= m_row)
581 {
582 ostringstream oss;
583 oss<<"In double Matrix::operator()(int i,int j) const.i "<<i<<" out of bounds.Please check it";
584 throw oss.str();
585 }
586 if(j<0 || j>= m_col)
587 {
588 ostringstream oss;
589 oss<<"In double Matrix::operator()(int i,int j) const.j "<<j<<" out of bounds.Please check it";
590 throw oss.str();
591 }
592
593 return m_data[i][j];
594}
595
596//to string 化为标准串
597string Matrix::ToString() const
598{
599 ostringstream oss;
600 for(int i=0;i<m_row;i++)
601 {
602 copy(m_data[i],m_data[i]+m_col,ostream_iterator<double>(oss," "));
603 oss<<"\n";
604 }
605
606 return oss.str();
607}
608
609// outputt stream function 输出
610
611ostream& operator<<(ostream& os,const Matrix& m)
612{
613 for(int i=0;i<m.m_row;i++)
614 {
615 copy(m.m_data[i],m.m_data[i]+m.m_col,ostream_iterator<double>(os," "));
616 os<<'\n';
617 }
618 return os;
619}
620
621// input stream function 输入
622istream& operator>>(istream& is,Matrix& m)
623{
624 for(int i=0;i<m.m_row;i++)
625 {
626 for(int j=0;j<m.m_col;j++)
627 {
628 is>>m.m_data[i][j];
629 }
630 }
631 return is;
632}
633
634
635// reallocation method
636
637// public method for resizing matrix
638
639// logical equal-to operator
640
641// logical no-equal-to operator
642
643// combined addition and assignment operator
644
645// combined subtraction and assignment operator
646
647// combined scalar multiplication and assignment operator
648
649
650// combined matrix multiplication and assignment operator
651
652// combined scalar division and assignment operator
653
654// combined power and assignment operator
655
656// unary negation operator
657
658// binary addition operator
659
660// binary subtraction operator
661
662// binary scalar multiplication operator
663
664// binary scalar multiplication operator
665
666// binary matrix multiplication operator
667
668// binary scalar division operator
669
670
671// binary scalar division operator
672
673// binary matrix division operator
674
675// binary power operator
676
677// unary transpose operator
678
679// unary Inverse operator
680
681// Inverse function
682
683// solve simultaneous equation
684
685// set zero to all elements of this matrix
686
687// set this matrix to unity
688
689// private partial pivoting method
690
691// calculate the determinant of a matrix
692
693// calculate the norm of a matrix
694
695// calculate the condition number of a matrix
696
697// calculate the cofactor of a matrix for a given element
698
699// calculate adjoin of a matrix
700
701// Determine if the matrix is singular
702
703// Determine if the matrix is diagonal
704
705// Determine if the matrix is scalar
706
707// Determine if the matrix is a unit matrix
708
709// Determine if this is a null matrix
710
711// Determine if the matrix is symmetric
712
713// Determine if the matrix is skew-symmetric
714// Determine if the matrix is upper triangular
715
716// Determine if the matrix is lower triangular
/***2
*3
* 4
* author:XieXiaokui5
* purpose:Defines functions for matrix6
*7
***/8
#include "stdafx.h"9
10
#include <iostream>11
#include <fstream>12
13
#include <string>14
#include <sstream>15
#include <algorithm>16
#include <functional>17
#include <numeric>18
#include <iterator>19
#include <cmath>20
#include <cassert>21
22
#include "matrix.h"23
24
using namespace std;25
26
const double Matrix::epsilon = 1e-7;27
28
// constructor29
30
Matrix::Matrix():m_row(0),m_col(0),m_data(0)31
{32
}33
34
35
// constructor36
Matrix::Matrix(int size):m_row(size),m_col(size)37
{38
if(size<=0)39
{
40
m_row=0;41
m_col=0;42
m_data=0;43
44
ostringstream oss;45
oss<<"In Matrix::Matrix(int size) size "<<size<<" <=0.Please check it.";46
throw oss.str();47
48
}49
50
m_data=new double*[size];51
52
for(int i=0;i<size;i++)53
{54
m_data[i]=new double[size];55
}56
57
}58
59
// constructor60
Matrix::Matrix(int row,int col):m_row(row),m_col(col)61
{62
if(row<=0)63
{64
m_row=0;65
m_col=0;66
m_data=0;67
68
ostringstream oss;69
oss<<"In Matrix::Matrix(int row,int col),row "<<row<<" <=0.Please check it.";70
throw oss.str();71
}72
if(col<=0)73
{74
m_row=0;75
m_col=0;76
m_data=0;77
78
ostringstream oss;79
oss<<" In Matrix::Matrix(int row,int col),col "<<col<<" <=0.Pleasecheck it.";80
throw oss.str();81
}82
83
m_data=new double*[row];84
for(int i=0;i<row;i++)85
{86
m_data[i]=new double[col];87
}88
}89
90
//copy constructor91
Matrix::Matrix(const Matrix& m):m_row(m.m_row),m_col(m.m_col)92
{93
m_data = new double*[m_row];94
for(int i=0;i<m_row;i++)95
{96
m_data[i] = new double[m_col];97
}98
99
for(int i=0;i<m_row;i++)100
{101
copy(m.m_data[i],m.m_data[i]+m_col,m_data[i]);102
}103
104
}105
106
Matrix::~Matrix()107
{108
if(m_data != 0)109
{110
for(int i=0;i<m_row;i++)111
{112
delete[] m_data[i];113
}114
delete[] m_data;115
}116
}117
118
void Matrix::SetZero()119
{120
for(int i=0;i<m_row;i++)121
fill(m_data[i],m_data[i]+m_col,0);122
}123
124
void Matrix::SetUnit()125
{126
for(int i=0;i<m_row;i++)127
for(int j=0;j<m_col;j++)128
m_data[i][j] = (i==j)?1:0;129
}130
131
void Matrix::SetSize(int size)132
{133
Clear();134
Create(size,size);135
}136
137
138
void Matrix::SetSize(int row,int col)139
{140
Clear();141
Create(row,col);142
}143
144
145
void Matrix::Clear()146
{147
148
if(m_data != 0)149
{ 150
for(int i=0;i<m_row;i++)151
{152
delete[] m_data[i];153
}154
delete[] m_data;155
}156
157
m_row=0;158
m_col=0;159
m_data=0;160
}161
162
/*163
Matrix& Matrix::Create(int size)164
{165
if(size<=0)166
{167
ostringstream oss;168
oss<<"In Matrix::Create(int size),size "<<size<<" <=0.Please check it.";169
170
throw oss.str();171
}172
173
if(m_data != 0)174
{175
for(int i=0;i<m_row;i++)176
{177
delete[] m_data[i];178
}179
delete[] m_data;180
}181
182
m_row=size;183
m_col=size;184
185
m_data=new double*[size];186
for(int i=0;i<size;i++)187
{188
m_data[i]=new double[size];189
}190
191
for(int i=0;i<size;i++)192
{193
for(int j=0;j<size;j++)194
{195
m_data[i][j] = ((i==j)?1:0);196
}197
}198
199
return *this;200
}201
*/202
203
void Matrix::Create(int row,int col)204
{205
if(row<=0)206
{207
ostringstream oss;208
oss<<"In Matrix::Create(int row,int col),row "<<row<<" <=0.Please check it.";209
throw oss.str();210
}211
if(col<=0)212
{213
ostringstream oss;214
oss<<"In Matrix::Create(int row,int col),col "<<col<<" <=0.Please check it.";215
throw oss.str();216
}217
218
if(m_data != 0)219
{220
for(int i=0;i<m_row;i++)221
{222
delete[] m_data[i];223
}224
delete[] m_data;225
}226
227
m_row=row;228
m_col=col;229
m_data=new double*[row];230
for(int i=0;i<row;i++)231
{232
m_data[i]=new double[col];233
}234
}235
236
237
int Matrix::GetRow() const238
{239
return m_row;240
}241
242
int Matrix::GetCol() const243
{244
return m_col;245
}246
247
//transpose 转置248
Matrix Matrix::Transpose() const249
{250
Matrix ret(m_col,m_row);251
for(int i=0;i<m_row;i++)252
{253
for(int j=0;j<m_col;j++)254
{255
ret.m_data[j][i] = m_data[i][j];256
}257
}258
return ret;259
}260
261
int Matrix::Pivot(int row) const262
{263
int index=row;264
265
for(int i=row+1;i<m_row;i++)266
{267
if(m_data[i][row] > m_data[index][row])268
index=i;269
}270
271
return index;272
}273
274
275
Matrix& Matrix::Exchange(int i,int j) // 初等变换：对调两行ri<-->rj276
{277
if(i<0 || i>=m_row)278
{279
ostringstream oss;280
oss<<"In void Matrix::Exchange(int i,int j) ,i "<<i<<" out of bounds.Please check it.";281
throw oss.str();282
}283
284
if(j<0 || j>=m_row)285
{286
ostringstream oss;287
oss<<"In void Matrix::Exchange(int i,int j) ,j "<<j<<" out of bounds.Please check it.";288
throw oss.str();289
}290
291
for(int k=0;k<m_col;k++)292
{293
swap(m_data[i][k],m_data[j][k]);294
}295
296
return *this;297
}298
299
300
Matrix& Matrix::Multiple(int index,double mul) //初等变换 第index 行乘以mul301
{302
if(index <0 || index >= m_row)303
{304
ostringstream oss;305
oss<<"In void Matrix::Multiple(int index,double k) index "<<index<<" out of bounds.Please check it.";306
throw oss.str();307
}308
309
transform(m_data[index],m_data[index]+m_col,m_data[index],bind2nd(multiplies<double>(),mul));310
return *this;311
}312
313
314
Matrix& Matrix::MultipleAdd(int index,int src,double mul) //初等变换 第src行乘以mul加到第index行315
{316
if(index <0 || index >= m_row)317
{318
ostringstream oss;319
oss<<"In void MultipleAdd(int index,int src,double k) index "<<index<<" out of bounds.Please check it.";320
throw oss.str();321
}322
323
if(src < 0 || src >= m_row)324
{325
ostringstream oss;326
oss<<"In void MultipleAdd(int index,int src,double k) src "<<src<<" out of bounds.Please check it.";327
throw oss.str();328
}329
330
for(int j=0;j<m_col;j++)331
{332
m_data[index][j] += m_data[src][j] * mul;333
}334
335
return *this;336
337
}338
339
340
//inverse 逆阵：使用矩阵的初等变换，列主元素消去法341
Matrix Matrix::Inverse() const342
{343
if(m_row != m_col) //非方阵344
{345
ostringstream oss;346
oss<<"In Matrix Matrix::Invert() const. m_row "<<m_row<<" != m_col "<<m_col<<" Please check it";347
throw oss.str();348
}349
350
Matrix tmp(*this);351
Matrix ret(m_row); //单位阵352
ret.SetUnit();353
354
int maxIndex;355
double dMul;356
357
for(int i=0;i<m_row;i++)358
{359
360
maxIndex = tmp.Pivot(i);361
362
if(tmp.m_data[maxIndex][i]==0)363
{364
ostringstream oss;365
oss<<"In Matrix Matrix::Invert() const 行列式的值等于0.Please check it";366
throw oss.str();367
}368
369
if(maxIndex != i) //下三角阵中此列的最大值不在当前行，交换370
{371
tmp.Exchange(i,maxIndex);372
ret.Exchange(i,maxIndex);373
374
}375
376
ret.Multiple(i,1/tmp.m_data[i][i]);377
tmp.Multiple(i,1/tmp.m_data[i][i]);378
379
380
for(int j=i+1;j<m_row;j++)381
{382
dMul = -tmp.m_data[j][i];383
tmp.MultipleAdd(j,i,dMul);384
ret.MultipleAdd(j,i,dMul);385
386
}387
388
}//end for389
390
for(int i=m_row-1;i>0;i--)391
{392
for(int j=i-1;j>=0;j--)393
{394
dMul = -tmp.m_data[j][i];395
tmp.MultipleAdd(j,i,dMul);396
ret.MultipleAdd(j,i,dMul);397
}398
}//end for399
400
return ret;401
}402
403
404
405
// assignment operator 賦值运算符406
Matrix& Matrix::operator= (const Matrix& m)407
{408
if(m_data != 0)409
{410
for(int i=0;i<m_row;i++)411
{412
delete[] m_data[i];413
}414
delete[] m_data;415
}416
417
m_row = m.m_row;418
m_col = m.m_col;419
420
m_data = new double*[m_row];421
for(int i=0;i< m_row;i++)422
{423
m_data[i] = new double[m_col];424
}425
426
for(int i=0;i<m_row;i++)427
{428
copy(m.m_data[i],m.m_data[i]+m_col,m_data[i]);429
}430
return *this;431
}432
433
434
//binary addition 矩阵加435
Matrix Matrix::operator+ (const Matrix& m) const436
{437
if(m_row != m.m_row)438
{439
ostringstream oss;440
oss<<"In Matrix::operator+(const Matrix& m) this->m_row "<<m_row<<" != m.m_row "<<m.m_row<<" Please check it.";441
throw oss.str();442
}443
if(m_col != m.m_col)444
{445
ostringstream oss;446
oss<<" Matrix::operator+ (const Matrix& m) this->m_col "<<m_col<<" != m.m_col "<<m.m_col<<" Please check it.";447
throw oss.str();448
}449
450
Matrix ret(m_row,m_col);451
for(int i=0;i<m_row;i++)452
{453
transform(m_data[i],m_data[i]+m_col,m.m_data[i],ret.m_data[i],plus<double>());454
}455
return ret;456
}457
458
//unary addition 求正459
Matrix Matrix::operator+() const460
{461
return *this;462
}463
464
//binary subtraction 矩阵减465
Matrix Matrix::operator- (const Matrix& m) const466
{467
if(m_row != m.m_row)468
{469
ostringstream oss;470
oss<<"In Matrix::operator-(const Matrix& m) this->m_row "<<m_row<<" != m.m_row "<<m.m_row<<" Please check it.";471
throw oss.str();472
}473
if(m_col != m.m_col)474
{475
ostringstream oss;476
oss<<"In Matrix::operator- (const Matrix& m) this->m_col "<<m_col<<" != m.m_col "<<m.m_col<<" Please check it.";477
throw oss.str();478
}479
480
Matrix ret(m_row,m_col);481
for(int i=0;i<m_row;i++)482
{483
transform(m_data[i],m_data[i]+m_col,m.m_data[i],ret.m_data[i],minus<double>());484
}485
return ret;486
}487
488
//unary substraction 求负489
Matrix Matrix::operator-() const490
{491
Matrix ret(*this);492
493
for(int i=0;i<m_row;i++)494
for(int j=0;j<m_col;j++)495
ret.m_data[i][j] *= -1;496
497
return ret;498
499
}500
501
502
//binary multiple 矩阵乘503
Matrix Matrix::operator*(const Matrix& m) const504
{505
if(m_col != m.m_row)506
{507
ostringstream oss;508
oss<<"In Matrix::operator*(const Matrix& m) this->m_col "<<m_col<<" != m.m_row "<<m.m_row<<" Please check it.";509
throw oss.str();510
}511
512
Matrix ret(m_row,m.m_col);513
Matrix tmp=m.Transpose();514
515
for(int i=0;i<m_row;i++)516
{517
for(int j=0;j<m.m_col;j++)518
{519
ret.m_data[i][j]=inner_product(m_data[i],m_data[i]+m_col,tmp.m_data[j],0.0);520
}521
}522
523
return ret;;524
}525
526
//friend scalar multiple 数乘527
Matrix operator* (double d,const Matrix& m)528
{529
Matrix ret(m);530
531
for(int i=0;i<ret.m_row;i++)532
for(int j=0;j<ret.m_col;j++)533
ret.m_data[i][j] *= d;534
535
return ret;536
537
}538
539
//binary matrix division equivalent to multiple inverse矩阵除，等价于乘以逆阵540
Matrix Matrix::operator/(const Matrix& m) const541
{542
return *this * m.Inverse();543
}544
545
//binary scalar division equivalent to multiple reciprocal数除，等价于乘以此数的倒数546
Matrix Matrix::operator/(double d) const547
{548
return 1/d * (*this);549
}550
551
//friend division552
Matrix operator/(double d,const Matrix& m)553
{554
return d * m.Inverse();555
}556
557
558
// subscript operator to get individual elements 下标运算符559
double Matrix::operator()(int i,int j) const560
{561
if(i<0 || i>= m_row)562
{563
ostringstream oss;564
oss<<"In double Matrix::operator()(int i,int j) const.i "<<i<<" out of bounds.Please check it";565
throw oss.str();566
}567
if(j<0 || j>= m_col)568
{569
ostringstream oss;570
oss<<"In double Matrix::operator()(int i,int j) const.j "<<j<<" out of bounds.Please check it";571
throw oss.str();572
}573
574
return m_data[i][j];575
}576
577
// subscript operator to set individual elements 下标运算符578
double& Matrix::operator ()(int i,int j)579
{580
if(i<0 || i>= m_row)581
{582
ostringstream oss;583
oss<<"In double Matrix::operator()(int i,int j) const.i "<<i<<" out of bounds.Please check it";584
throw oss.str();585
}586
if(j<0 || j>= m_col)587
{588
ostringstream oss;589
oss<<"In double Matrix::operator()(int i,int j) const.j "<<j<<" out of bounds.Please check it";590
throw oss.str();591
}592
593
return m_data[i][j];594
}595
596
//to string 化为标准串597
string Matrix::ToString() const598
{599
ostringstream oss;600
for(int i=0;i<m_row;i++)601
{602
copy(m_data[i],m_data[i]+m_col,ostream_iterator<double>(oss," "));603
oss<<"\n";604
}605
606
return oss.str();607
}608
609
// outputt stream function 输出610
611
ostream& operator<<(ostream& os,const Matrix& m)612
{613
for(int i=0;i<m.m_row;i++)614
{615
copy(m.m_data[i],m.m_data[i]+m.m_col,ostream_iterator<double>(os," "));616
os<<'\n';617
}618
return os;619
}620
621
// input stream function 输入622
istream& operator>>(istream& is,Matrix& m)623
{624
for(int i=0;i<m.m_row;i++)625
{626
for(int j=0;j<m.m_col;j++)627
{628
is>>m.m_data[i][j];629
}630
}631
return is;632
}633
634
635
// reallocation method636
637
// public method for resizing matrix638
639
// logical equal-to operator640
641
// logical no-equal-to operator642
643
// combined addition and assignment operator644
645
// combined subtraction and assignment operator646
647
// combined scalar multiplication and assignment operator648
649
650
// combined matrix multiplication and assignment operator651
652
// combined scalar division and assignment operator653
654
// combined power and assignment operator655
656
// unary negation operator657
658
// binary addition operator659
660
// binary subtraction operator661
662
// binary scalar multiplication operator663
664
// binary scalar multiplication operator665
666
// binary matrix multiplication operator667
668
// binary scalar division operator669
670
671
// binary scalar division operator672
673
// binary matrix division operator674
675
// binary power operator676
677
// unary transpose operator678
679
// unary Inverse operator680
681
// Inverse function682
683
// solve simultaneous equation684
685
// set zero to all elements of this matrix686
687
// set this matrix to unity688
689
// private partial pivoting method690
691
// calculate the determinant of a matrix692
693
// calculate the norm of a matrix694
695
// calculate the condition number of a matrix696
697
// calculate the cofactor of a matrix for a given element698
699
// calculate adjoin of a matrix700
701
// Determine if the matrix is singular702
703
// Determine if the matrix is diagonal704
705
// Determine if the matrix is scalar706
707
// Determine if the matrix is a unit matrix708
709
// Determine if this is a null matrix710
711
// Determine if the matrix is symmetric712
713
// Determine if the matrix is skew-symmetric714
// Determine if the matrix is upper triangular715
716
// Determine if the matrix is lower triangular另有标准C语言版本。
