MATRIX.H
上传用户:wyp_nj
上传日期:2022-05-03
资源大小:484k
文件大小:8k
源码类别:

GPS编程

开发平台:

Visual C++

  1. #ifndef MATRIX_H
  2. #define MATRIX_H
  3. // 定义适合数学运算的实数matrix类,数据将存放在buffer类中
  4. #include <iostream.h>
  5. #include <iomanip.h>
  6. #include <stdio.h>
  7. // buffer.h包含实数缓存类buffer的定义
  8. #include "buffer.h"
  9. // 实数矩阵类,每个矩阵有rownum行,colnum列,总共rownum乘colnum个实数
  10. // 下标一律从0开始起算,即存在第0行第0列
  11. class matrix {
  12.  public:
  13. buffer * buf;  // 缓存区指针,指向存放实数的地方,可以是任何媒体
  14. // buffer类是一抽象类,具体的子类可以将数据存放在内存或者
  15. // 临时文件中,用户还可以设计自己的缓存区子类
  16. size_t rownum, colnum;  // 矩阵的行数和列数
  17. unsigned istrans:1; // 转置标志,0为缺省,1为转置
  18. unsigned isneg:1; // 为负标志,0为非负,1为取负
  19. unsigned issym:1; // 对称标志,0为非对称,1为对称,不具有约束力
  20. // 在下面的所有构造函数中,缓存区指针b通常不设,缺省为0,这样程序将根据
  21. // 自动申请缺省的缓存区类型。如果用户定义了特殊的缓顾虑区类型,可以
  22. // 先申请后在构造函数中指明
  23. matrix(buffer * b=0); // 缺省构造函数,产生0行0列空矩阵
  24. matrix(size_t n, buffer * b=0); // 产生nX1列向量
  25. matrix(size_t r, size_t c, buffer * b=0); // 构造函数,产生r行c列矩阵
  26. matrix(matrix& m); //拷贝构造函数,产生的矩阵内容同m相同,
  27. // 新产生的矩阵并不申请新的缓存,而是指向与m同样的缓存,
  28. // 只是它们共同指向的缓存的引用数增1,这样可以节省存储资源
  29. matrix(const char * filename, buffer * b=0); /* 文件构造函数,
  30. 使用数据文件对矩阵初使化,数据文件为文本文件,格式:
  31. 以小写的matrix关键词打头,然后是空格,再跟矩阵的行数,再跟
  32. 空格,再跟列数,再跟空格,接着是内容,
  33. 内容按行分,首先是行号,从1开始计数,行号后面紧跟着冒号:,
  34. 然后是这行的数据,各个数据之间都用空格或者换行隔开 */
  35. matrix(void * data, size_t r, size_t c=1, buffer * b=0); /* 数据构造函数,
  36. 行数为r,列数为c,data必须指向一存放实数的内存,按先列后行的
  37. 顺序进行存放。构造函数进行数据拷贝,因此在构造函数完成后,
  38. data指向的数据区可以释放 */
  39. DOUBLE operator()(size_t r, size_t c); // 重载函数运算符()返回第r行第c列的实数,
  40. DOUBLE operator()(size_t r){return (*this)(r,0);}; // 取回第r行0列的实数
  41. virtual void set(size_t r, size_t c, DOUBLE v) {
  42. size_t l;
  43. v = isneg ? -v : v;
  44. l = istrans ? r+c*rownum : r*colnum+c;
  45. buf = buf->set(l, v); }; // 设置第r行第c列的值为v,
  46. // 在设置时如果缓存的引用数大于1,产生新的缓存
  47. // 因此每次设置返回的缓存区指针都要再赋给buf。
  48. void set(size_t r, DOUBLE v) {
  49. set(r,0,v);}; // 设置第r行0列的值为v,用作列向量
  50. virtual ~matrix(){ // 析构函数,如果缓存的引用数大于1,说明还有其它的
  51. // 矩阵在用此缓存,则不释放缓存,只是将引用数减1
  52. buf->refnum--;
  53. if(!buf->refnum) delete buf;};
  54. matrix& operator=(matrix& m); // 重载赋值运算符,使矩阵的内容等于另一个矩阵
  55. matrix& operator=(DOUBLE a); // 通过赋值运算符将矩阵所有元素设为a
  56. matrix operator*(DOUBLE a); // 重载乘法运算符,实现矩阵与数相乘,矩阵在前
  57. // 数在后
  58. matrix& operator*=(DOUBLE a);  // 重载*=运算符,功能同上,但*产生了一个新
  59. // 的矩阵,原矩阵不变,而本运算符修改了原矩阵
  60. friend matrix operator*(DOUBLE a, matrix& m);
  61. // 数乘矩阵,但数在前矩阵在后
  62. matrix operator+(matrix& m); // 矩阵相加,原二矩阵不变并产生一个新的和矩阵
  63. matrix& operator+=(matrix &m); // 矩阵相加,不产生新矩阵,有一矩阵的内容改变为和
  64. matrix operator+(DOUBLE a); // 矩阵加常数,指每一元素加一固定的常数,产生
  65. // 新矩阵,原矩阵不变
  66. friend matrix operator+(DOUBLE a, matrix& m); // 常数加矩阵,产生新的和矩阵
  67. matrix& operator+=(DOUBLE a); // 矩阵自身加常数,自身内容改变
  68. matrix operator*(matrix& m); // 矩阵相乘,产生出新的矩阵
  69. matrix& operator*=(matrix& m); // 矩阵相乘,结果修改原矩阵
  70. matrix operator-(); // 矩阵求负,产生新的负矩阵
  71. matrix & trans(){
  72. size_t r = rownum; rownum = colnum; colnum = r;
  73. istrans = !istrans; return (*this);}; // 矩阵自身转置
  74. matrix t(); // 矩阵转置,产生新的矩阵
  75. matrix & neg(){isneg = !isneg; return (*this);}; // 将自己求负
  76. matrix operator-(matrix& m); // 矩阵相减,产生新的矩阵
  77. matrix& operator-=(matrix& m); // 矩阵相减,结果修改原矩阵
  78. matrix operator-(DOUBLE a)  // 矩阵减常数,指每一元素减一固定的常数,产生
  79. // 新矩阵,原矩阵不变
  80. { return (*this)+(-a); };
  81. friend matrix operator-(DOUBLE a, matrix& m); // 常数减矩阵,产生新的矩阵
  82. matrix& operator-=(DOUBLE a) // 矩阵自身减常数,自身内容改变
  83. { return operator+=(-a);};
  84. int isnear(matrix& m, double e = defaulterr); // 检查两个矩阵是否近似相等,
  85. // 如近似相等则返回1,否则返回0,e为允许误差
  86. int isnearunit(double e = defaulterr); // 检查矩阵是否近似为单位矩阵
  87. // 如是则返回1,否则返回0,e为允许误差
  88. matrix row(size_t r); // 提取第r行行向量
  89. matrix col(size_t c); // 提取第c列列向量
  90. void checksym(); // 检查和尝试调整矩阵到对称矩阵
  91. void swapr(size_t r1, size_t r2, size_t k=0); // 交换两行,只交换k列以上
  92. void swapc(size_t c1, size_t c2, size_t k=0); // 交换两列,只交换k行以上
  93. void swap(size_t r1, size_t c1, size_t r2, size_t c2){ // 交换两个元素
  94. DOUBLE a; a=(*this)(r1,c1);set(r1,c1,(*this)(r2,c2));
  95. set(r2,c2,a);};
  96. DOUBLE maxabs(size_t &r, size_t &c, size_t k=0);
  97.  // 求主元值和位置,即k行k列之后的最大元素,
  98. // 元素所在的行列将放在r,c中,返回此元素值
  99. size_t zgsxy(matrix & m, int fn=0); // 主高斯消元运算
  100. matrix& operator/=(matrix m); // 用主高斯消元法求解线性方程的解
  101. // 矩阵本身为系数矩阵,m为常数向量,返回解向量
  102. matrix operator/(matrix m); // 左乘m的逆矩阵
  103. matrix& inv(); // 用全选主元高斯-约当法求逆矩阵
  104. matrix operator~(); // 求逆矩阵,但产生新矩阵
  105. friend matrix operator/(DOUBLE a, matrix& m); // 求逆矩阵再乘常数
  106. matrix& operator/=(DOUBLE a); // 所有元素乘a的倒数,自身改变
  107. matrix operator/(DOUBLE a); // 所有元素乘a的倒数,产生新的矩阵
  108. matrix cholesky(matrix& d); // 用乔里斯基分解法求对称正定阵的线性
  109. // 方程组ax=d返回方程组的解,本身为a,改变为分解阵u,d不变
  110. DOUBLE det(DOUBLE err=defaulterr); // 求行列式的值
  111. size_t rank(DOUBLE err=defaulterr); // 求矩阵的秩
  112. void house(buffer & b, buffer & c); // 用豪斯荷尔德变换将对称阵变为对称三对
  113. // 角阵,b返回主对角线元素,c返回次对角线元素
  114. void trieigen(buffer& b, buffer& c, size_t l=600, DOUBLE eps=defaulterr);
  115. // 计算三对角阵的全部特征值与特征向量
  116. matrix eigen(matrix & eigenvalue, DOUBLE eps=defaulterr, size_t l=600);
  117.  // 计算对称阵的全部特征向量和特征值
  118. // 返回按列排放的特征向量,而eigenvalue将返回一维矩阵,为所有特征值
  119. // 组成的列向量
  120. int ispositive(); // 判定矩阵是否对称非负定阵,如是返回1,否则返回0
  121. void hessenberg(); // 将一般实矩阵约化为赫申伯格矩阵
  122. void qreigen(matrix & a, matrix & b, size_t l=600, DOUBLE eps=defaulterr);
  123. // 求一般实矩阵的所有特征根
  124. // a和b均返回rownum行一列的列向量矩阵,返回所有特征根的实部和虚部
  125. friend ostream& operator<<(ostream& o, matrix& m);
  126. friend istream& operator>>(istream& in, matrix& m);
  127. virtual DOUBLE value(size_t r, size_t c){ // 返回第r行c列的数值,不做行列检查
  128. DOUBLE a;
  129. a = istrans ? (*buf)[r+c*rownum] : (*buf)[r*colnum + c];
  130. return isneg ? -a : a;
  131. };
  132.  protected:
  133. DOUBLE detandrank(size_t & r, DOUBLE err); // 求方阵的行列式和秩
  134. };
  135. inline matrix operator*(DOUBLE a, matrix& m) {
  136. return m*a;
  137. };
  138. inline matrix operator+(DOUBLE a, matrix& m) { // 常数加矩阵,产生新的和矩阵
  139. return m+a;
  140. };
  141. matrix operator-(DOUBLE a, matrix& m); // 常数减矩阵,产生新的矩阵
  142. matrix operator/(DOUBLE a, matrix& m); // 求逆矩阵再乘常数
  143. ostream& operator<<(ostream& o, matrix& m);
  144. istream& operator>>(istream& in, matrix& m);
  145. matrix unit(size_t n); // 产生n阶单位矩阵
  146. DOUBLE gassrand(int rr=0); // 返回一零均值单位方差的正态分布随机数
  147. // rr为种子
  148. class gassvector //返回零均值给定协方差阵的正态随机向量的类
  149. {
  150.  public:
  151. matrix a; // a是增益矩阵,由协方差矩阵算出
  152. gassvector(matrix & r); //r必须是正定对称阵,为正态随机向量的协方差
  153. matrix operator()(matrix & r); // 返回给定协方差矩阵的正态随机向量
  154. matrix operator()(); // 返回已设定的协方差矩阵的正态随机向量
  155. void setdata(matrix & r); // 根据协方差矩阵设置增益矩阵
  156. };
  157. DOUBLE lineopt(matrix& a,matrix& b, matrix& c, matrix & x);
  158. // 线性规划最优点寻找程序,a为mXn不等式约束条件左端系数矩阵,b为不等式约束
  159. // 条件的右端项,为m维向量,c为目标函数系数,n维向量,x返回极小点,为n维向量
  160. #endif // MATRIX_H