数学库（持续更新中 18-05-06）
命名约定
/** 命名类型统一 ** ** ** ** **/ #ifndef CPLX_NAMED #define CPLX_NAMED #include <iostream> using namespace std; #pragma message("cplx_name头文件在调试") typedef int cint; typedef double creal; typedef long double creall; typedef void* cpoint; enum CBusiness{ Row = 1, Col = 2 }; enum CFunction{ Z, EXP_Z, CONJ_Z, INVERS_Z }; enum COperateType { Operate_Type_Unknow = -1, // Operate_Type_Sy = 1, // Symbol Operate_Type_Num, // Number Operate_Type_Fun, // Function Operate_Type_Op, // Opearte Operate_Type_Bu, // Bound }; enum CVar_Type { Var_Type_Unknow = 0, Var_Type_Num, Var_Type_Fun, Var_Type_Handle, Var_Type_Array }; #endif // CPLX_NAMED
基本库
基本库的大部分是开学初写的，现在看看觉得当时的代码好幼稚，可能这就是成长吧（被自己的厚脸皮逗笑）
因为存在一些项目已经用到了里面的一些内容，所以在此就不对其进行更新了，术语向后兼容
新添加了一个vecII 的容器类，虽然没有vector支持多线程直接安全使用那么屌的功能，但也还算不错啦（vector的 at() 比 operator[] 快是因为 at() 使用的是读锁，不用读互斥排队啥的，像我这种渣渣是写不出来这样的程序的）
之前的vec最好不要用了吧，说不定有bug（直白嫌弃开学时的我），但因为某些项目还受他的支持，就暂时不删吧。
1 /* 复数操作的核心模块（2018-03-26） 2 * 3 * 定义内容包括： 4 * 基本复数、复数容器对象 5 * 复函数对象 6 * 2\3维函数定义域对象 7 * 矩阵操作基本容器 vecII 18-05-06 8 * 9 * 已使用及可兼容的接口包括： 10 * 11 * 12 * 13 * 14 * 15 * dusty-cjh@qq.com 16 */ 17 #ifndef CPLX_CORE 18 #define CPLX_CORE 19 #include "cplx_named.h" 20 #include <initializer_list> 21 #include <stdexcept> 22 #include <math.h> 23 #include <string> 24 #ifdef QT_CORE_LIB 25 #include <QString> 26 #include <QDebug> 27 #endif 28 29 30 using std::runtime_error; 31 using std::initializer_list; 32 using std::string; 33 34 35 /*宏命名*/ 36 #define CORE_INT cint 37 #define CORE_REAL creal 38 #define CORE_REALL creall 39 #define CORE_POINT cpoint 40 41 /*宏量*/ 42 #define Pi 3.14159265358979323846 43 #define E 2.71828182845904523536 44 #define Cir (2*Pi) 45 #define ROUND 360 46 #define I(type) (cplx<type>(0,1)) 47 #define IS_INF(x) isinf((x)) 48 #define IS_NAN(x) isnan((x)) 49 #define IS_NORMAL(x) isnormal((x)) 50 51 52 /*宏函数*/ 53 #define C_MOVE(x) (std::move((x))) 54 #define RADIAN(x, y) std::atan2(y, x); 55 #define ANGLE(rad) (rad)/Cir*ROUND 56 57 58 59 // 前置 60 template<typename _Tp> struct _polar; 61 template<typename _Tp> struct _cplx; 62 template<typename _Tp> class cplx; 63 template<typename _Tp> class region_mult; 64 template<typename _Tp> class region_const; 65 template<typename _Geo> struct _geometry; 66 67 // 外部前置 68 template<typename _T> class _mtx; 69 70 // 常用类型 71 typedef cplx<double> Cplx; 72 73 74 75 /*定义域参数方程 point type */ 76 template<typename _T> struct defintion_refer{ 77 typedef _cplx<_T>(*fun_type)(_T); 78 }; 79 #define DEFINE_FUN_CAST(type) static_cast<defintion_refer<type>::fun_type> 80 /*映射函数 point type */ 81 template<typename _T> struct mapping_refer{ 82 typedef _T(*fun_type)(_T); 83 }; 84 #define MAPPING_FUN_CAST(type) static_cast<mapping_refer<type>::fun_type> 85 86 87 88 89 // base function 90 /* @模 */ 91 template<typename _T> inline _T cmold(const _cplx<_T>& c1) throw(){ 92 return sqrt(std::pow(c1.real, 2) + std::pow(c1.imag, 2)); 93 } 94 /* @角度 */ 95 template<typename _T> inline _T cangle(const _cplx<_T>& c1) throw(){ 96 return ANGLE(atan2(c1.imag, c1.real)); 97 } 98 /* @弧度 */ 99 template<typename _T> inline _T cradian(const _cplx<_T>& c1) throw(){ 100 return atan2(c1.imag, c1.real); 101 } 102 /* @点乘 */ 103 template<typename _T> inline _T cdot(const _cplx<_T>& c1, const _cplx<_T>& c2) throw(){ 104 return c1.real*c2.real + c1.imag*c2.imag; 105 } 106 /* @叉乘 */ 107 template<typename _T> inline _T ccross(const _cplx<_T>& c1, const _cplx<_T>& c2) throw(){ 108 // cross(c1,c2) equal x1*y2-x2*y1 109 return c1.real*c2.imag - c1.imag*c2.real; 110 } 111 /* @极坐标 */ 112 template<typename _T> inline _polar<_T> cto_polar(const _cplx<_T>& c1) throw(){ 113 return{ 114 cmold(c1), 115 atan2(c1.imag, c1.real) 116 }; 117 } 118 /* @复平面坐标 */ 119 template<typename _T> inline cplx<_T> cto_cplx(const _polar<_T>& p1) throw(){ 120 return{ 121 p1.r*cos(p1.a), 122 p1.r*sin(p1.a) 123 }; 124 } 125 /* @实数化尝试 */ 126 template<typename _T> inline _T cto_real(const _cplx<_T>& c1){ 127 return c1.imag ? nan(nullptr) : c1.real; 128 } 129 template<typename _T> inline _T cto_real(const _polar<_T>& p1){ 130 return p1.a%Cir ? p1.r : nan(nullptr); 131 } 132 133 134 //------------------------------------------- math mathod ----------------------------------------- 135 136 template<typename _T> inline cplx<_T> exp(cplx<_T> z)throw() 137 { 138 return{ 139 exp(z.real)*cos(z.imag), 140 exp(z.real)*sin(z.imag) 141 }; 142 } 143 template<typename _T> inline cplx<_T> sinh(cplx<_T> z)throw() 144 { 145 return (exp(z) - exp(-z)) / 2; 146 } 147 template<typename _T> inline cplx<_T> cosh(cplx<_T> z)throw() 148 { 149 return (exp(z) + exp(-z)) / 2; 150 } 151 template<typename _T> inline cplx<_T> sin(cplx<_T> z)throw() 152 { 153 return sinh(z / (cplx<_T>(0, 1)))*(cplx<_T>(0, 1)); 154 } 155 template<typename _T> inline cplx<_T> cos(cplx<_T> z)throw() 156 { 157 auto a = exp(z*I(_T)), 158 b = exp(-z*I(_T)), 159 c = a + b, 160 d = c / 2, 161 f = exp(z*I(_T)) + exp(-z*I(_T)), 162 g = (exp(z*I(_T)) + exp(-z*I(_T))) / 2; 163 return ((exp(z*I(_T)) + exp(-z*I(_T))) / 2); 164 } 165 template<typename _T> inline cplx<_T> tan(cplx<_T> z){ 166 return sin(z) / cos(z); 167 } 168 template<typename _T> inline cplx<_T> log(cplx<_T> z)throw(){ 169 cplx<_T> Tmp; 170 Tmp = I(_T)*z.imag + log(cto_polar(z).r); 171 return Tmp; 172 } 173 template<typename _T> inline cplx<_T> pow(cplx<_T> z, _T n){ 174 _polar<_T> por = cto_polar(z); 175 por.r = pow(por.r, n); 176 por.a *= n; 177 178 return cto_cplx(por); 179 } 180 template<typename _T> inline cplx<_T> sqrt(cplx<_T> z){ 181 return pow(z, 0.5); 182 } 183 184 185 186 // regulrary 187 template<typename _T> std::string cto_string(_T c) 188 { 189 return std::to_string(c); 190 } 191 template<typename _T> std::string cto_string(const cplx<_T>& c) 192 { 193 return c.to_string(); 194 } 195 template<typename _T> std::string cto_string(const region_const<_T>& reg) 196 { 197 string str; 198 for (int r = 0; r < reg.rows; ++r){ 199 str += "\nRows " + cto_string(r + 1); 200 for (int c = 0; c < reg.cols; ++c){ 201 str += "\n"; 202 str += std::to_string(reg(r, c)); 203 } 204 } 205 206 return str; 207 } 208 template<typename _T> std::string cto_string(const region_const<cplx<_T>>& reg) 209 { 210 string str; 211 for (int r = 0; r < reg.rows; ++r){ 212 str += "\nRows " + cto_string(r + 1); 213 for (int c = 0; c < reg.cols; ++c){ 214 str += "\n"; 215 str += cto_string(reg(r, c)); 216 } 217 str += '\n'; 218 } 219 220 return str; 221 } 222 template<typename _T> std::string cto_string(const _geometry<_T>& geo){ 223 string str; 224 str += "geometry(" + cto_string(geo.x) + ", " + 225 cto_string(geo.y) + ", " + 226 cto_string(geo.width) + ", " + 227 cto_string(geo.height) + ") "; 228 return str; 229 } 230 #ifdef QT_CORE_LIB 231 template<typename _T> QString cto_qstring(region_const<_T>& reg) 232 { 233 QString str; 234 for (int r = 0; r < reg.rows; ++r){ 235 for (int c = 0; c < reg.cols; ++c){ 236 str += '\t'; 237 str += QString::number(reg(r, c)); 238 } 239 str += '\n'; 240 } 241 242 return str; 243 } 244 template<typename _T> QString cto_qstring(region_const<cplx<_T>>& reg) 245 { 246 QString str; 247 for (int r = 0; r < reg.rows; ++r){ 248 for (int c = 0; c < reg.cols; ++c){ 249 str += "\t"; 250 str += cto_qstring(reg(r, c)); 251 str += cto_qstring(reg(r, c)); 252 } 253 str += '\n'; 254 } 255 256 return str; 257 } 258 template<typename _T> QString cto_qstring(_T& c) 259 { 260 return QString::number(c); 261 } 262 template<typename _T> QString cto_qstring(cplx<_T>& c) 263 { 264 QString cplx_str = QString("complex(real:") + QString::number(c.real) 265 + QString(", imag:") + QString::number(c.imag) + QString(") "); 266 267 return cplx_str; 268 } 269 #endif 270 271 272 // base complex only support the container of complex. 273 template<typename _Tp = CORE_REAL> 274 struct _cplx 275 { 276 _Tp real; 277 _Tp imag; 278 }; 279 // describtion of complex in polar 280 template<typename _Tp = CORE_REAL> 281 struct _polar 282 { 283 _Tp r; 284 _Tp a; 285 }; 286 // support like conj,to\set_polar,radius,angle... 287 template<typename _Tp> 288 class cplx : public _cplx<_Tp> 289 { 290 public: 291 typedef _Tp refer; 292 typedef const refer& const_refer; 293 typedef refer&& move_refer; 294 typedef _cplx<refer> base_refer; 295 typedef const _cplx<refer> const_base_refer; 296 typedef _cplx<refer>&& move_base_refer; 297 298 cplx(refer r = 0, refer i = 0) 299 :_cplx({ r, i }) 300 {} 301 cplx(_cplx& c) 302 :_cplx(c) 303 {} 304 cplx(const cplx& c) = default; 305 cplx& operator=(const cplx& c) = default; 306 cplx& operator=(_cplx& c){ 307 real = c.real; 308 imag = c.imag; 309 } 310 operator _cplx<refer>(){ 311 return{ 312 real, 313 imag 314 }; 315 } 316 317 318 void set_polar(refer r, refer a){ 319 real = r*cos(a); 320 imag = r*sin(a); 321 } 322 void set_polar(_polar<refer> p){ 323 real = p.r*cos(p.a); 324 imag = p.r*sin(p.a); 325 } 326 _polar<refer> polar()const{ 327 return{ 328 cmold(*this), 329 cradian(*this) 330 }; 331 } 332 refer r()const{ 333 return cmold(*this); 334 } 335 refer a()const{ 336 return cradian(*this); 337 } 338 cplx<refer> conj(){ 339 return cplx(real, -imag); 340 } 341 342 string to_string(void)const{ 343 string cplx_str = string("complex(real:") + std::to_string(real) 344 + string(", imag:") + std::to_string(imag) + string(") "); 345 346 return cplx_str; 347 } 348 #ifdef QT_CORE_LIB 349 QString to_qstring(void){ 350 return QString("complex(real:") + QString::number(real) + QString(", imag:") 351 + QString::number(imag) + QString(") "); 352 } 353 #endif 354 355 static cplx i(void){ 356 return cplx(0, 1); 357 } 358 359 public: // 基本运算 360 inline cplx operator*(const cplx<_Tp>& c1)const throw(){ 361 cplx Tmp = *this; 362 Tmp *= c1; 363 364 return Tmp; 365 } 366 inline cplx operator*(cplx<_Tp>&& c1)const throw(){ 367 cplx Tmp = *this; 368 Tmp *= c1; 369 370 return Tmp; 371 } 372 inline cplx operator*(_Tp x)const throw(){ 373 cplx Tmp = *this; 374 Tmp *= x; 375 376 return Tmp; 377 } 378 inline cplx operator/(const cplx& c1)const throw(){ 379 cplx Tmp = *this; 380 Tmp /= c1; 381 382 return Tmp; 383 } 384 inline cplx operator/(cplx<_Tp>&& c1)const throw(){ 385 cplx Tmp = *this; 386 Tmp /= c1; 387 388 return Tmp; 389 } 390 inline cplx operator/(_Tp x)const throw(){ 391 cplx Tmp = *this; 392 Tmp /= x; 393 394 return Tmp; 395 } 396 inline cplx operator+(const cplx& c1)const throw(){ 397 cplx Tmp = *this; 398 Tmp += c1; 399 400 return Tmp; 401 } 402 inline cplx operator+(cplx&& c1)const throw(){ 403 cplx Tmp = *this; 404 Tmp += c1; 405 406 return Tmp; 407 } 408 inline cplx operator+(_Tp x)const throw(){ 409 cplx Tmp = *this; 410 Tmp += x; 411 412 return Tmp; 413 } 414 inline cplx operator-(const cplx& c1)const throw(){ 415 cplx Tmp = *this; 416 Tmp -= c1; 417 418 return Tmp; 419 } 420 inline cplx operator-(_Tp x)const throw(){ 421 cplx Tmp = *this; 422 Tmp -= x; 423 424 return Tmp; 425 } 426 inline cplx operator-(cplx&& c1)const throw(){ 427 cplx Tmp = *this; 428 Tmp -= c1; 429 430 return Tmp; 431 } 432 433 inline cplx operator-(void)const throw(){ 434 return{ 435 -real, 436 -imag 437 }; 438 } 439 inline cplx operator+(void)const throw(){ 440 return *this; 441 } 442 inline cplx& operator+=(const cplx& c1) throw(){ 443 real += c1.real; 444 imag += c1.imag; 445 446 return *this; 447 } 448 inline cplx& operator+=(_Tp x) throw(){ 449 real += x; 450 451 return *this; 452 } 453 inline cplx& operator-=(const cplx& c1) throw(){ 454 real -= c1.real; 455 imag -= c1.imag; 456 457 return *this; 458 } 459 inline cplx& operator-=(_Tp x) throw(){ 460 real -= x; 461 462 return *this; 463 } 464 inline cplx& operator*=(const cplx& c1) throw(){ 465 _Tp tmp = c1.real*real - c1.imag*imag; 466 imag = c1.real*imag + c1.imag*real; 467 real = tmp; 468 469 return *this; 470 } 471 inline cplx& operator*=(_Tp x) throw(){ 472 real *= x; 473 imag *= x; 474 475 return *this; 476 } 477 inline cplx& operator/=(const cplx& c1) throw(){ 478 _Tp mu = pow(c1.real, 2) + pow(c1.imag, 2); 479 _Tp tmp; 480 tmp = (c1.real*real + c1.imag*imag) / mu; 481 imag = (c1.real*imag - c1.imag*real) / mu; 482 real = tmp; 483 484 return *this; 485 } 486 inline cplx& operator/=(_Tp x) throw(){ 487 real /= x; 488 imag /= x; 489 490 return *this; 491 } 492 inline bool operator==(const cplx& c1)const{ 493 return (real == c1.real) && (imag == c1.imag); 494 } 495 inline bool operator!=(const cplx& c1)const{ 496 return !operator==(c1); 497 } 498 499 500 }; 501 502 503 504 505 //--------------------------------------- defintion --------------------------------------------- 506 507 /* @定义域在坐标系中的矩形 */ 508 template<typename _Geo> 509 struct _geometry 510 { 511 _Geo x; 512 _Geo y; 513 _Geo width; 514 _Geo height; 515 }; 516 517 518 /*! 网格定义域 */ 519 template<typename _T> region_const<cplx<_T>> grid_region(cplx<_T> tl, cplx<_T> br, _T prec = 0.01){ 520 // 网格区域 grid_region(左上顶点坐标, 右下顶点坐标, 精度 = 0.01) 521 int r = static_cast<int>(abs(tl.imag - br.imag) / prec) + 1, 522 c = static_cast<int>(abs(tl.real - br.real) / prec) + 1; 523 524 region_const<cplx<_T>> region(r, c); 525 cplx<_T>* p = region.data; 526 for (int i = 0; i < r; ++i) 527 for (int j = 0; j < c; ++j) 528 *p++ = cplx<_T>(tl.real + j*prec, tl.imag - i*prec); 529 530 return region; 531 } 532 /*! 椭圆定义域 */ 533 template<typename _T> region_const<cplx<_T>> ellipse_region(_T a, _T b, _T angle_prec = 0.01, _T len_prec = 0.1) 534 {// 椭圆网格区域 ellipse_region (长轴， 短轴, 精度 = 0.01 )； 535 if (!angle_prec || !len_prec) 536 return region_const<cplx<_T>>(); 537 538 region_const<cplx<_T>> reg(static_cast<_T>(Cir / angle_prec), static_cast<_T>((a > b ? b : a) / len_prec)); 539 _T a_step = a / reg.cols, 540 b_step = b / reg.cols; 541 for (int r = 0; r<reg.rows; ++r){ 542 for (int c = 0; c<reg.cols; ++c){ 543 reg(r, c) = cplx<_T>(a_step*c*cos(angle_prec*r), b_step*c*sin(angle_prec*r)); 544 } 545 } 546 547 return reg; 548 } 549 550 /*! 计算定义域矩形geometry*/ 551 template<typename _T> _geometry<_T> geometry(region_const<cplx<_T>>& reg){ 552 if (!reg.amount()) 553 return{ 554 0, 0, 555 0, 0 556 }; 557 558 auto co = reg.data[0]; 559 _T x_min = co.real, 560 x_max = co.real, 561 y_min = co.imag, 562 y_max = co.imag; 563 for (int i = 0; i<reg.rows*reg.cols; ++i){ 564 co = reg.data[i]; 565 if (co.real > x_max) 566 std::swap(co.real, x_max); 567 else if (co.real < x_min) 568 std::swap(co.real, x_min); 569 if (co.imag > y_max) 570 std::swap(co.imag, y_max); 571 else if (co.imag < y_min) 572 std::swap(co.imag, y_min); 573 } 574 575 _geometry<_T> g = { 576 x_min, 577 y_max, 578 abs(x_max - x_min), 579 abs(y_max - y_min) 580 }; 581 return g; 582 } 583 584 585 586 587 588 /* 基本容器类 589 * 590 */ 591 template<typename _Cnt> 592 class _vec 593 { 594 public: 595 typedef _Cnt refer; 596 typedef _Cnt& using_refer; 597 typedef const _Cnt& const_refer; 598 typedef _Cnt&& move_refer; 599 explicit _vec(int length = 0, refer* d = nullptr) 600 :len(length), 601 data(d) 602 { 603 if (d) 604 d = nullptr; 605 else 606 data = length ? new refer[static_cast<int>(length)] : nullptr; 607 raise = 0; 608 } 609 _vec(_vec&& v) 610 :len(std::move(v.len)), 611 data(std::move(v.data)), 612 raise(std::move(v.raise)) 613 {} 614 _vec(const _vec&) = delete; 615 _vec(initializer_list<refer> il){ 616 data = nullptr; 617 this->operator =(il); 618 } 619 ~_vec(){ 620 delete[]data; 621 } 622 623 _vec& operator=(_vec&& v){ 624 len = std::move(v.len); 625 data = std::move(v.data); 626 raise = std::move(v.raise); 627 628 return *this; 629 } 630 _vec& operator=(const _vec& v){ 631 this->len = v.len; 632 this->raise = v.raise; 633 634 if (data) 635 delete[]data; 636 data = new refer[len + raise]; 637 for (int i = 0; i < len; ++i) 638 data[i] = v.data[i]; 639 640 return *this; 641 } 642 _vec& operator=(initializer_list<refer> il){ 643 resize(il.size()); 644 refer* p = data; 645 for (auto beg = il.begin(); beg != il.end(); ++beg, ++p) 646 *p = *beg; 647 648 return *this; 649 } 650 651 652 653 const_refer at(int i)const{ 654 return data[i]; 655 } 656 refer* atP(int i){ 657 return &data[i]; 658 } 659 using_refer operator[](int i){ 660 return data[i]; 661 } 662 663 // move 664 void append(const _vec& v){ 665 if (v.len > raise){ 666 int length = (int)(len + v.len); 667 raise = static_cast<int>(length*0.5); 668 669 refer* d = new refer[length + raise]; 670 for (int i = 0; i < length; ++i){ 671 d[i] = data[i]; 672 } 673 for (int i = len; i < len + v.len; ++i){ 674 d[i] = v.data[i - len]; 675 } 676 677 delete[]data; 678 data = d; 679 len = length; 680 } 681 else{ 682 for (int i = len; i < len + v.len; ++i){ 683 data[i] = v.data[i - len]; 684 } 685 raise -= v.len; 686 } 687 } 688 void append(const refer& c){ 689 if (1 > raise){ 690 raise = static_cast<int>((1 + len)*0.5); 691 692 refer* d = new refer[len + 1 + raise]; 693 for (int i = 0; i < len; ++i){ 694 d[i] = data[i]; 695 } 696 d[len] = c; 697 698 delete[]data; 699 data = d; 700 ++len; 701 } 702 else{ 703 data[len] = c; 704 raise -= 1; 705 } 706 } 707 void remove(int front, int endafter){ 708 if (endafter < len){ 709 for (int i = endafter; i < len; ++i){ 710 data[front + i - endafter] = data[i]; 711 } 712 len -= endafter - front; 713 raise += endafter - front; 714 } 715 else{ 716 len = front; 717 raise += endafter - front; 718 } 719 } 720 void insert(int before, const _vec& v){ 721 int length = static_cast<int>(len + v.len); 722 refer* d = new refer[static_cast<int>(length*1.5)]; 723 for (int i = 0; i < before; ++i){ 724 d[i] = data[i]; 725 } 726 for (int i = 0; i < v.len; ++i){ 727 d[before + i] = v.data[i]; 728 } 729 for (int i = 0; i < len - before; ++i){ 730 d[before + v.len + i] = data[before + i]; 731 } 732 733 delete[]data; 734 data = d; 735 len = length; 736 raise = static_cast<int>(length*0.5); 737 } 738 739 740 // redo 741 void clear(void){ 742 len = raise = 0; 743 delete[]data; 744 data = nullptr; 745 } 746 void resize(int size){ 747 len = size; 748 raise = 0; 749 delete[]data; 750 data = new refer[size]; 751 } 752 _vec clone(void){ 753 _vec cln(this->len); 754 for (int i = 0; i < len; ++i){ 755 cln[i] = (*this)[i]; 756 } 757 } 758 759 public: 760 refer* data; 761 int len; 762 int raise; 763 764 }; 765 766 /* 为标准矩阵设计的vecII（推荐使用vecII） 767 * 2018-04-25 18:12 （电工课已迟到） 768 * 769 * warning！ 模板参数需要可安全析构多次（在初始化内存时，直接拷贝了已初始化的对象数据） 770 * 判别标准为：凡调用除copy外的任何赋值操作时，数据所有权都将被转移 771 * 772 * 与所有我写的库一样，vecII不检验任何参数错误，绝对相信程序员 773 * 774 * C++ 的运算符重载对面向对象编程有很大影响。 775 * 关键在于不要期望代码能达到用户操作的方便程度，那样编写只会轻而易举地落入混乱的圈套 776 * 正确的方式是把运算符重载看成是加速编程的一种手段，除非对象的操作十分明显（迭代器用++，vec用[]）否则绝不用运算符重载 777 */ 778 template<typename _Cnt> 779 class _vecII 780 { 781 /*宏定义了一些重复代码，看参数和函数名就知道意思了*/ 782 #define VECII_SET(obj, length, raise_size) { \ 783 (obj).dt_p = (refer*)calloc((length + raise_size), sizeof(refer)); \ 784 (obj).len = length; \ 785 (obj).raise = raise_size; \ 786 refer *it = (obj).dt_p, \ 787 *end = (obj).dt_p + (obj).len; \ 788 while (it != end) \ 789 *it++ = refer(); \ 790 } /*order the default constructer for safe*/ 791 #define VECII_SET_VAL(obj, data_p, length, raise_size) { \ 792 (obj).dt_p = data_p; \ 793 (obj).len = length; \ 794 (obj).raise = raise_size; \ 795 } 796 #define VECII_CLONE(obj, src) { \ 797 (obj).dt_p = (refer*)calloc((src).len, sizeof(refer)); \ 798 (obj).len = (src).len; \ 799 (obj).raise = 0; \ 800 for (int i = 0; i < (src).len; ++i) \ 801 (obj)[i] = (src).at(i); \ 802 } 803 #define VECII_MOVE(obj, v) { \ 804 (obj).dt_p = (v).dt_p; \ 805 (obj).len = (v).len; \ 806 (obj).raise = (v).raise; \ 807 (v).dt_p = nullptr; \ 808 (v).len = 0; \ 809 (v).raise = 0; \ 810 } 811 #define VECII_COPY(dst, src, ele_count) { \ 812 for (int i = 0; i < (ele_count); ++i) \ 813 (dst)[i] = (src)[i]; \ 814 } 815 #define VECII_MOVE_MEMORY(dst, src, ele_count) { \ 816 for (int i = 0; i < (ele_count); ++i) \ 817 (dst)[i] = std::move((src)[i]); \ 818 } 819 #define VECII_NEW_SIZE(cur_len) (static_cast<int>(1.5*((cur_len)+1))) 820 #define VECII_NEW_MEMORY(ele_count) (_memory_init(ele_count)) 821 /*宏结束*/ 822 823 public: 824 typedef _Cnt refer; 825 typedef const refer const_refer; 826 typedef refer* iter; 827 typedef const_refer* const_iter; 828 template<typename _Other> friend class _mtx; 829 830 explicit _vecII(int len = 0){ 831 if (0 == len){ 832 this->len = 0; 833 this->raise = 0; 834 this->dt_p = nullptr; 835 }else 836 VECII_SET(*this, len, 0); 837 } 838 _vecII(const _vecII& v){ 839 VECII_SET_VAL(*this, VECII_NEW_MEMORY(v.len), v.len, 0); 840 VECII_COPY(dt_p, v.dt_p, v.len); 841 } 842 _vecII(_vecII&& v){ 843 VECII_MOVE(*this, v); 844 } 845 virtual ~_vecII(){ 846 clear(); 847 } 848 _vecII& operator=(const _vecII& v){ 849 clear(); 850 VECII_SET_VAL(*this, VECII_NEW_MEMORY(v.len), v.len, 0); 851 VECII_COPY(dt_p, v.dt_p, v.len); 852 853 return *this; 854 } 855 _vecII& operator=(_vecII&& v){ 856 clear(); 857 VECII_MOVE(*this, v); 858 859 return *this; 860 } 861 862 863 public: 864 // memory 865 void clear(){ 866 free(dt_p); 867 dt_p = nullptr; 868 len = raise = 0; 869 } 870 void scale_to_fit(){ 871 if (!raise) 872 return; 873 refer *p = VECII_NEW_MEMORY(this->len); 874 VECII_MOVE_MEMORY(p, dt_p, this->len); 875 free(dt_p); 876 VECII_SET_VAL(*this, p, len, 0); 877 } 878 void resize(int new_len){ // 重新设置大小 879 if (new_len == this->len + this->raise) 880 return; 881 if (new_len < this->len + this->raise){ 882 this->raise += len - new_len; 883 this->len = new_len; 884 return; 885 } 886 int data_len = new_len; 887 888 new_len = VECII_NEW_SIZE(data_len); // 计算新大小 889 refer* p = VECII_NEW_MEMORY(new_len); 890 VECII_MOVE_MEMORY(p, dt_p, data_len > len ? len : new_len) // 转移该块内存所有权 891 free(this->dt_p); 892 893 VECII_SET_VAL(*this, p, data_len, new_len - data_len); 894 } 895 void reset(int new_len){ // 格式化 896 if (new_len == this->len + this->raise) 897 return; 898 clear(); 899 VECII_SET(*this, new_len, 0); 900 } 901 902 // quote 903 inline refer& operator[](int pos)throw(){ 904 return memory(pos); 905 } 906 const_refer& at(int pos)const{ return dt_p[pos]; } 907 inline refer& memory(int pos)throw(){ // 返回索引项的引用 908 return this->dt_p[pos]; 909 } 910 iter begin(){ // 返回存储指针 911 return this->dt_p; 912 } 913 iter end(){ 914 return begin() + len; 915 } 916 const_iter cbegin()const{ // 返回存储指针 917 return this->dt_p; 918 } 919 const_iter cend()const{ 920 return begin() + len; 921 } 922 int size()const{ return this->len; } // 已有元素数量 923 924 // element 925 _vecII copy(int beg, int end){ // 此函数为唯一可复制值的函数，拷贝某一区域的值，返回新构建的容器 926 _vecII Tmp; 927 VECII_SET(Tmp, (end - beg), 0); 928 VECII_COPY(Tmp.dt_p, &this->memory(beg), end - beg); 929 930 return Tmp; 931 } 932 _vecII clone(void)const{ // 返回当前对象的克隆体 933 _vecII Tmp; 934 VECII_CLONE(Tmp, *this); 935 936 return Tmp; 937 } 938 _vecII move(void){ 939 return C_MOVE(*this); 940 } 941 942 // element operate 943 // warning! copy itself equal to suicide 944 int insert(int before, const _vecII& v){ 945 if (this->raise < v.len){ // no enough space 946 int new_size = static_cast<int>(1.5*(len +v.len)); 947 refer *p = VECII_NEW_MEMORY(new_size); 948 949 // copy 950 VECII_MOVE_MEMORY(p, this->dt_p, before); 951 VECII_MOVE_MEMORY(p + before + v.len, dt_p+before, this->len - before); 952 953 // free memory 954 free(this->dt_p); 955 this->dt_p = p; 956 this->len += v.len; 957 this->raise = new_size - this->len; 958 959 } 960 else{ // 961 int L = this->len - before; 962 int x = v.len; 963 refer *p = &memory(before + x); 964 965 while (0 <= L / x){// 多拷一次，省一堆事，对性能影响不大 966 VECII_MOVE_MEMORY(p, p - x, x);/*这样的话存在析构未初始化区域的问题，所以最好在使用内存前就用calloc全部归零，不过依然不安全*/ 967 L -= x; 968 p -= x; 969 } 970 } 971 972 VECII_COPY(dt_p + before, v.dt_p, v.len); // 插入数据 973 return before; 974 } 975 int insert(int before, _vecII&& v){ 976 if (this->raise < v.len){ // no enough space 977 int new_size = static_cast<int>(1.5*(len + v.len)); 978 refer *p = VECII_NEW_MEMORY(new_size); 979 980 // copy 981 VECII_MOVE_MEMORY(p, this->dt_p, before); 982 VECII_MOVE_MEMORY(p + before + v.len, dt_p + before, this->len - before); 983 984 // free memory 985 free(this->dt_p); 986 this->dt_p = p; 987 this->len += v.len; 988 this->raise = new_size - this->len; 989 990 } 991 else{ // 992 int L = this->len - before; 993 int x = v.len; 994 refer *p = &memory(before + x); 995 996 while (0 <= L / x){// 多拷一次，省一堆事，对性能影响不大 997 VECII_MOVE_MEMORY(p, p - x, x);/*这样的话存在析构未初始化区域的问题，所以最好在使用内存前就用calloc全部归零，不过依然不安全*/ 998 L -= x; 999 p -= x; 1000 } 1001 } 1002 1003 VECII_MOVE_MEMORY(dt_p + before, v.dt_p, v.len); // 插入数据 1004 free(v.dt_p); 1005 VECII_SET_VAL(v, nullptr, 0, 0); 1006 return before; 1007 } 1008 int insert(int before, refer& ele){ 1009 _vecII Tmp(1); 1010 Tmp.dt_p[0] = ele; 1011 1012 return insert(before, Tmp); 1013 } 1014 int append(const _vecII& v){ // 插入对象 1015 if (v.len > this->raise){ 1016 // expand memory 1017 int new_size = static_cast<int>(1.5*(v.len + this->len)); 1018 refer *p = VECII_NEW_MEMORY(new_size); 1019 1020 // copy 1021 VECII_MOVE_MEMORY(p, dt_p, this->len); 1022 1023 // clean 1024 free(dt_p); 1025 VECII_SET_VAL(*this, p, len + v.len, new_size - len - v.len); 1026 1027 // append data 1028 VECII_COPY(dt_p + len - v.len, v.dt_p, v.len); 1029 } 1030 else{ 1031 VECII_COPY(dt_p + len, v.dt_p, v.len); 1032 } 1033 1034 return len - v.len; 1035 } 1036 int append(_vecII&& v){ 1037 if (v.len > this->raise){ 1038 // expand memory 1039 int new_size = static_cast<int>(1.5*(v.len + this->len)); 1040 refer *p = VECII_NEW_MEMORY(new_size); 1041 1042 // copy 1043 VECII_MOVE_MEMORY(p, dt_p, this->len); 1044 1045 // clean 1046 free(dt_p); 1047 VECII_SET_VAL(*this, p, len + v.len, new_size - len - v.len); 1048 1049 // append data 1050 VECII_COPY(dt_p + len - v.len, v.dt_p, v.len); 1051 } 1052 else{ 1053 VECII_MOVE_MEMORY(dt_p + len, v.dt_p, v.len); 1054 free(v.dt_p); 1055 VECII_SET_VAL(v, nullptr, 0, 0); 1056 } 1057 1058 return len - v.len; 1059 } 1060 int append(refer& ele){ // 插入元素 1061 if (raise){ 1062 --raise; 1063 dt_p[len++] = ele; 1064 } 1065 else{ 1066 int new_size = VECII_NEW_SIZE(len); 1067 refer* p = VECII_NEW_MEMORY(new_size); 1068 VECII_MOVE_MEMORY(p, dt_p, this->len); 1069 p[this->len] = ele; 1070 1071 // clean 1072 free(dt_p); 1073 VECII_SET_VAL(*this, p, len + 1, new_size - len - 1); 1074 } 1075 1076 return len - 1; 1077 } 1078 int append(refer ele){ 1079 if (raise){ 1080 --raise; 1081 dt_p[len++] = ele; 1082 } 1083 else{ 1084 int new_size = VECII_NEW_SIZE(len); 1085 refer* p = VECII_NEW_MEMORY(new_size); 1086 VECII_MOVE_MEMORY(p, dt_p, this->len); 1087 p[this->len] = ele; 1088 1089 // clean 1090 free(dt_p); 1091 VECII_SET_VAL(*this, p, len + 1, new_size - len - 1); 1092 } 1093 1094 return len - 1; 1095 } 1096 int remove(int beg, int end){ // 移除某一区域的所有元素 1097 VECII_MOVE_MEMORY(dt_p + beg, dt_p + end, len - end); 1098 raise += end - beg; 1099 len -= end - beg; 1100 1101 return beg; 1102 } 1103 inline int remove(int pos){ // 移除该点元素 1104 return remove(pos, pos + 1); 1105 } 1106 void exchange(int pos1, int pos2){ // 交换两元素位置 1107 refer Tmp = dt_p[pos1]; 1108 dt_p[pos1] = dt_p[pos2]; 1109 dt_p[pos2] = Tmp; 1110 } 1111 void exchange(int dst_pos, int beg, int length){ // 交换两位置一定长度的元素 1112 refer *Tmp = VECII_NEW_MEMORY(length); 1113 1114 VECII_MOVE_MEMORY(Tmp, &dt_p[beg], length); 1115 VECII_MOVE_MEMORY(&dt_p[beg], &dt_p[dst_pos], length); 1116 VECII_MOVE_MEMORY(&dt_p[dst_pos], Tmp, length); 1117 1118 free(Tmp); 1119 } 1120 1121 private: 1122 int len; 1123 int raise; 1124 refer* dt_p; 1125 1126 private: /** memory safe request and destroy **/ 1127 static refer* _memory_init(int size){ 1128 refer Tmp = refer(); // 先初始化模板，使其具有良好的初始结构，要求对象的析构函数可安全调用多次 1129 refer *p = (refer*)malloc(sizeof(refer)*size); 1130 for (refer *it = p, *end = p + size; 1131 it != end; ++it) 1132 memcpy(it, &Tmp, sizeof(refer));// 较为安全的方法，直接拷贝初始状态的对象内存 1133 1134 return p; 1135 } 1136 static refer* _memory_move(refer* src, int ele_count){ 1137 refer *p = (refer*)malloc(sizeof(refer)*ele_count); 1138 for (int i = 0; i < ele_count; ++i) 1139 p[i] = std::move(src[i]); 1140 1141 return p; 1142 } 1143 static refer* _memory_copy(refer* src, int ele_count){ 1144 refer *p = (refer*)malloc(sizeof(refer)*ele_count); 1145 for (int i = 0; i < ele_count; ++i) 1146 p[i] = src[i]; 1147 1148 return p; 1149 } 1150 static void _memory_destroy(refer* p){ 1151 free(p); 1152 } 1153 1154 public:/** 后续添加项 **/ 1155 // 1156 // 1157 // 请在此处添加自定义函数 1158 // 1159 1160 #undef VECII_SET // 初始化各项 1161 #undef VECII_SET_VAL // 直接设置新值 1162 #undef VECII_MOVE // 移动某对象所有权 1163 #undef VECII_COPY // 等同memcpy，但是对对象的操作，使用拷贝构造函数 1164 #undef VECII_NEW_SIZE // 容器自然增长的新大小 1165 #undef VECII_NEW_MEMORY // 申请内存 1166 }; 1167 1168 /* 一个可变大小的网格定义域 1169 * 1170 */ 1171 template<typename _Df> 1172 class region_mult 1173 { 1174 public: 1175 friend class region_const<_Df>; 1176 typedef _Df refer; 1177 typedef _Df& using_refer; 1178 typedef const _Df& const_refer; 1179 typedef _Df&& move_refer; 1180 typedef _vec<_vec<_Df>> vec_type; 1181 typedef _vec<_Df>& row_refer; 1182 1183 1184 // iter 1185 class row_mult 1186 { 1187 public: 1188 explicit row_mult(int r, region_mult& reg) 1189 :row(r), 1190 reg(reg) 1191 {} 1192 row_mult(const row_mult& rows) 1193 :row(rows.row), 1194 reg(rows.reg) 1195 {} 1196 row_mult(const region_mult& reg){ 1197 setValue(reg); 1198 } 1199 ~row_mult(){ 1200 1201 } 1202 1203 using_refer operator[](int i){ 1204 return reg(row, i); 1205 } 1206 row_mult& operator=(const row_mult& rows){ 1207 this->row = rows.row; 1208 this->reg = rows.reg; 1209 } 1210 row_mult& operator=(const region_mult& reg){ 1211 setValue(reg); 1212 } 1213 row_mult& operator*=(refer n){ 1214 for (int i = 0; i<reg.cols; ++i){ 1215 (*this)[i] *= n; 1216 } 1217 } 1218 row_mult& operator/=(refer n){ 1219 for (int i = 0; i<reg.cols; ++i){ 1220 (*this)[i] /= n; 1221 } 1222 } 1223 row_mult& operator+=(refer n){ 1224 for (int i = 0; i<reg.cols; ++i){ 1225 (*this)[i] += n; 1226 } 1227 } 1228 row_mult& operator-=(refer n){ 1229 for (int i = 0; i<reg.cols; ++i){ 1230 (*this)[i] -= n; 1231 } 1232 } 1233 void setValue(const row_mult& rows){ 1234 for (int i = 0; i<reg.cols; ++i) 1235 (*this)[i] = rows[i]; 1236 } 1237 void setValue(const region_mult& reg){ 1238 if (reg.rows != 1 || reg.cols != this->reg.cols) 1239 throw(runtime_error("region_mult is not confitrable for row_mult")); 1240 for (int i = 0; i<reg.cols; ++i) 1241 (*this)[i] = reg(0, i); 1242 } 1243 1244 1245 public: 1246 int row; 1247 region_mult& reg; 1248 }; 1249 class col_mult 1250 { 1251 public: 1252 explicit col_mult(int r, region_mult& reg) 1253 :col(r), 1254 reg(reg) 1255 {} 1256 col_mult(const col_mult& cols) 1257 :col(cols.col), 1258 reg(cols.reg) 1259 {} 1260 1261 public: 1262 int col; 1263 region_mult& reg; 1264 }; 1265 1266 1267 1268 1269 1270 1271 explicit region_mult(int rows = 0, int cols = 0) 1272 { 1273 resize(rows, cols); 1274 } 1275 1276 // read 1277 using_refer operator()(int r, int c){ 1278 return v.data[r].data[c]; 1279 } 1280 const_refer at(int r, int c){ 1281 return v.data[r].data[c]; 1282 } 1283 1284 // 1285 void resize(int r, int c){ 1286 v.resize(r); 1287 for (int i = 0; i < v.len; ++i) 1288 v.data[i].resize(c); 1289 rows = r; 1290 cols = c; 1291 } 1292 region_mult rowClone(int r){ 1293 region_mult<_Df> reg(1, cols); 1294 reg.v = v[r].clone(); 1295 1296 return reg; 1297 } 1298 region_mult colClone(int c){ 1299 region_mult<_Df> reg(rows, 1); 1300 for (int i = 0; i<rows; ++i){ 1301 reg[i] = this->operator ()(i, c); 1302 } 1303 1304 return reg; 1305 } 1306 row_mult row(int r){ 1307 return row_mult(r, *this); 1308 } 1309 col_mult col(int c){ 1310 return col_mult(c, *this); 1311 } 1312 1313 public: 1314 int rows; 1315 int cols; 1316 _vec<_vec<_Df>> v; 1317 1318 }; 1319 1320 template<typename _Df> 1321 class region_const 1322 { 1323 public: 1324 friend class region_mult<_Df>; 1325 typedef _Df refer; 1326 typedef refer* pointer; 1327 typedef refer& using_refer; 1328 typedef const refer&const_refer; 1329 typedef refer&& move_refer; 1330 1331 1332 // row iterator 1333 class row_iter 1334 { 1335 public: 1336 explicit row_iter(int r, region_const& vec) 1337 :row(r), 1338 vec(vec), 1339 it(&vec.data[r*vec.cols]) 1340 {} 1341 row_iter(const row_iter& iter) 1342 :row(iter.row), 1343 it(iter.it), 1344 vec(iter.vec) 1345 {} 1346 row_iter(row_iter&& iter) 1347 :row(C_MOVE(iter.row)), 1348 it(C_MOVE(iter.it)), 1349 vec(C_MOVE(iter.vec)) 1350 {} 1351 row_iter& operator=(const row_iter& iter){ 1352 row = iter.row; 1353 it = iter.it; 1354 vec = iter.vec; 1355 } 1356 row_iter& operator=(row_iter&& iter){ 1357 row = std::move(iter.row); 1358 it = std::move(iter.it); 1359 vec = std::move(iter.vec); 1360 } 1361 1362 row_iter& operator+=(int i){ 1363 it += i; 1364 1365 return *this; 1366 } 1367 row_iter& operator-=(int i){ 1368 it -= i; 1369 1370 return *this; 1371 } 1372 row_iter operator+(int i){ 1373 row_iter Tmp = *this; 1374 Tmp += i; 1375 1376 return Tmp; 1377 } 1378 row_iter operator-(int i){ 1379 row_iter Tmp = *this; 1380 Tmp -= i; 1381 1382 return Tmp; 1383 } 1384 row_iter& operator++(void){ 1385 *this += 1; 1386 1387 return *this; 1388 } 1389 row_iter operator++(int){ 1390 row_iter Tmp = *this; 1391 *this += 1; 1392 1393 return Tmp; 1394 } 1395 row_iter& operator--(void){ 1396 *this -= 1; 1397 1398 return *this; 1399 } 1400 row_iter operator--(int){ 1401 row_iter Tmp = *this; 1402 *this -= 1; 1403 1404 return Tmp; 1405 } 1406 1407 using_refer operator*(void){ 1408 return *it; 1409 } 1410 using_refer operator[](int i){ 1411 return vec(row, i); 1412 } 1413 bool operator==(const row_iter& iter){ 1414 return (it == iter.it); 1415 } 1416 bool operator!=(const row_iter& iter){ 1417 return (it != iter.it); 1418 } 1419 1420 1421 private: 1422 int row; 1423 pointer it; 1424 region_const& vec; 1425 1426 }; 1427 // col iterator 1428 class col_iter 1429 { 1430 public: 1431 explicit col_iter(int c, region_const& vec) 1432 :col(c), 1433 vec(vec), 1434 it(&vec.data[c]) 1435 {} 1436 col_iter(const col_iter& iter) 1437 :col(iter.col), 1438 it(iter.it), 1439 vec(iter.vec) 1440 {} 1441 col_iter(col_iter&& iter) 1442 :col(C_MOVE(iter.col)), 1443 it(C_MOVE(iter.it)), 1444 vec(C_MOVE(iter.vec)) 1445 {} 1446 col_iter& operator=(const col_iter& iter){ 1447 col = iter.col; 1448 it = iter.it; 1449 vec = iter.vec; 1450 } 1451 col_iter& operator=(col_iter&& iter){ 1452 col = std::move(iter.col); 1453 it = std::move(iter.it); 1454 vec = std::move(iter.vec); 1455 } 1456 1457 col_iter& operator+=(int i){ 1458 it += vec.cols*i; 1459 1460 return *this; 1461 } 1462 col_iter& operator-=(int i){ 1463 it -= vec.cols*i; 1464 1465 return *this; 1466 } 1467 col_iter operator+(int i){ 1468 col_iter Tmp = *this; 1469 Tmp += i; 1470 1471 return Tmp; 1472 } 1473 col_iter operator-(int i){ 1474 col_iter Tmp = *this; 1475 Tmp -= i; 1476 1477 return Tmp; 1478 } 1479 col_iter& operator++(void){ 1480 *this += 1; 1481 1482 return *this; 1483 } 1484 col_iter operator++(int){ 1485 col_iter Tmp = *this; 1486 *this += 1; 1487 1488 return Tmp; 1489 } 1490 col_iter& operator--(void){ 1491 *this -= 1; 1492 1493 return *this; 1494 } 1495 col_iter operator--(int){ 1496 col_iter Tmp = *this; 1497 *this -= 1; 1498 1499 return Tmp; 1500 } 1501 1502 using_refer operator*(void){ 1503 return *it; 1504 } 1505 using_refer operator[](int i){ 1506 return vec(i, col); 1507 } 1508 bool operator==(const col_iter& iter){ 1509 return (it == iter.it); 1510 } 1511 bool operator!=(const col_iter& iter){ 1512 return (it != iter.it); 1513 } 1514 1515 1516 private: 1517 int col; 1518 pointer it; 1519 region_const& vec; 1520 1521 }; 1522 1523 1524 explicit region_const(int r = 0, int c = 0, pointer d = nullptr) 1525 :rows(r), 1526 cols(c), 1527 data(d) 1528 { 1529 // row first 1530 if (nullptr == data) 1531 data = new refer[r*c]; 1532 } 1533 region_const(region_const& reg) 1534 :rows(reg.rows), 1535 cols(reg.cols), 1536 data(reg.data) 1537 { 1538 if (nullptr == data) 1539 data = new refer[rows*cols]; 1540 else{ 1541 reg.data = nullptr; 1542 reg.rows = reg.cols = 0; 1543 } 1544 } 1545 region_const(const initializer_list<_Df>& il){ 1546 rows = il.size(); 1547 cols = 1; 1548 data = new refer[rows]; 1549 pointer p = data; 1550 for (auto beg = il.begin(); beg != il.end(); ++beg){ 1551 *p++ = *beg; 1552 } 1553 1554 } 1555 region_const(const initializer_list<initializer_list<_Df>>& il){ 1556 rows = il.size(); 1557 cols = il.begin()->size(); 1558 data = new refer[rows*cols]; 1559 pointer it = data; 1560 for (auto beg = il.begin(); beg != il.end(); ++beg){ 1561 if (beg->size() != cols) 1562 throw runtime_error("region initlizer should be squire"); 1563 for (auto b = beg->begin(); b != beg->end(); ++b){ 1564 *it = *b; 1565 ++it; 1566 } 1567 } 1568 1569 } 1570 ~region_const(){ 1571 delete[]data; 1572 } 1573 region_const& operator=(region_const& reg){ 1574 delete[]data; 1575 rows = reg.rows; 1576 cols = reg.cols; 1577 data = reg.data; 1578 1579 if (nullptr == data) 1580 data = new refer[rows*cols]; 1581 else{ 1582 reg.data = nullptr; 1583 reg.rows = reg.cols = 0; 1584 } 1585 1586 return *this; 1587 } 1588 region_const& operator=(const initializer_list<_Df>& il){ 1589 rows = il.size(); 1590 cols = 1; 1591 delete[]data; 1592 data = new refer[rows]; 1593 pointer p = data; 1594 for (auto beg = il.begin(); beg != il.end(); ++beg){ 1595 *p++ = *beg; 1596 } 1597 1598 } 1599 region_const& operator=(const initializer_list<initializer_list<_Df>>& il){ 1600 rows = il.size(); 1601 cols = il.begin()->size(); 1602 delete[]data; 1603 data = new refer[rows*cols]; 1604 pointer it = data; 1605 for (auto beg = il.begin(); beg != il.end(); ++beg){ 1606 if (beg->size() != cols) 1607 throw runtime_error("region initlizer should be squire"); 1608 for (auto b = beg->begin(); b != beg->end(); ++b){ 1609 *it++ = *b; 1610 } 1611 } 1612 1613 } 1614 1615 using_refer operator()(int r, int c){ 1616 return data[r*cols + c]; 1617 } 1618 const_refer at(int r, int c){ 1619 return data[r*rows + c]; 1620 } 1621 row_iter row(int r){ 1622 return row_iter(r, *this); 1623 } 1624 col_iter col(int c){ 1625 return col_iter(c, *this); 1626 } 1627 1628 region_const clone(void)const{ 1629 pointer p = new refer[rows*cols]; 1630 for (int i = 0; i < rows*cols; ++i) 1631 p[i] = data[i]; 1632 1633 return region_const(rows, cols, p); 1634 } 1635 inline int amount()const{ 1636 return rows*cols; 1637 } 1638 inline bool operator==(const region_const& reg){ 1639 return (reg.rows == rows && reg.cols == cols); 1640 } 1641 inline bool operator!=(const region_const& reg){ 1642 return !(*this == reg); 1643 } 1644 1645 public: 1646 refer *data; 1647 int rows; 1648 int cols; 1649 1650 }; 1651 1652 1653 1654 1655 //------------------------------------------- mapping function ----------------------------------------- 1656 1657 1658 1659 /*! 参数方程初始化定义域 */ 1660 template<typename _T> 1661 void defintionMapped(region_const<_T>& seed, region_const<cplx<_T>>& field, CORE_POINT fun = nullptr) 1662 {// seed:参数路径 field:值域对象 fun:映射函数，默认为_info的全局实例对象的指针 1663 int seg = seed.rows*seed.cols, 1664 fig = field.rows*field.cols; 1665 1666 // 异常处理 1667 try{ 1668 if (seg != fig) 1669 throw runtime_error("seed and field must be the same size."); 1670 if (nullptr == fun) 1671 throw runtime_error("function point is nullptr."); 1672 } 1673 catch (int){ 1674 // do noting 1675 } 1676 1677 for (int i = 0; i < seg; ++i){ 1678 field.data[i] = DEFINE_FUN_CAST(_T)(fun)(seed.data[i]); 1679 } 1680 1681 return; 1682 } 1683 /*! 函数映射 */ 1684 template<typename _T> 1685 void mappingMapped(region_const<_T>& seed, region_const<_T>& field, CORE_POINT fun = nullptr) 1686 {// seed:原象路径 field:值域对象 fun:映射函数，默认为_info的全局实例对象的指针 1687 try{ 1688 if (seed.rows*seed.cols != field.rows*field.cols) 1689 throw runtime_error("seed and field must be the same size."); 1690 if (nullptr == fun) 1691 throw runtime_error("function point is nullptr."); 1692 } 1693 catch (int){ 1694 // do nothing 1695 } 1696 1697 for (int i = 0, len = seed.rows*seed.cols; i < len; ++i){ 1698 field.data[i] = MAPPING_FUN_CAST(_T)(fun)(seed.data[i]); 1699 } 1700 1701 return; 1702 } 1703 1704 1705 1706 1707 #endif // CPLX_CORE
基本矩阵存储及数值运算矩阵
实现了矩阵存储模板，一些方便的函数都已经定义好了，还算完善
实现了数值矩阵基本操作，以double为运算类型，vecII为容器，凭借pthread实现了简单的运算多线程，比单线程的效率高2倍（依然很差，matlab感觉秒算出来了，可能我算法不够好）
1 #pragma once 2 3 #ifndef C_MTX 4 #define C_MTX 5 #define C_KERNAL 4 6 #include "cplx_core.h" 7 #include <pthread.h> 8 #include <exception> 9 #include <memory> 10 #include <vector> 11 #include <Windows.h> 12 #pragma comment(lib, "x86/pthreadVC2.lib") 13 14 #define MIN(x1, x2) ( (x1)>(x2) ? (x2) : (x1) ) 15 #define MAX(x1, x2) ( (x1)<(x2) ? (x2) : (x1) ) 16 17 template<typename _T> 18 class _mtx; 19 class _martix; 20 21 typedef _mtx<creal> Mtx; // 默认数值矩阵为creal类型 22 23 24 25 class _mtx_err : public std::logic_error 26 { 27 #define MTX_ERR(reason) _mtx_err(reason) 28 #define MTX_ERR_DSC(fun_name, detail) _mtx_err((#fun_name" "#detail)) /*错误描述（所在函数 + 详细）*/ 29 #define MTX_ERR__DIM_ERR(fun_name) MTX_ERR_DSC(fun_name, dim not cope) 30 template<typename _Other> friend class _mtx; 31 friend class _martix; 32 explicit _mtx_err(const std::string &s) 33 :std::logic_error(s) 34 {} 35 }; 36 37 /* 标准矩阵存储模板 38 39 1.列优先模式 40 2.以vecII为列的存储对象 41 42 实现了： 43 行列的引用 44 新建子矩阵 45 添加、交换、删除矩阵元素 46 47 */ 48 template<typename _T> 49 class _mtx 50 { 51 #define MTX_SET(obj, r, c) { \ 52 (obj).data.reset(c); \ 53 for (int i = 0; i < c; ++i) \ 54 (obj).data[i].reset(r); \ 55 (obj).rows = r; \ 56 (obj).cols = c; \ 57 } 58 #define MTX_SET_VAL(obj, r, c, dt) { \ 59 (obj).rows = r; \ 60 (obj).cols = c; \ 61 (obj).data = dt;\ 62 } 63 #define MTX_MOVE(dst, src) { \ 64 (dst).data = std::move((src).data); \ 65 (dst).rows = (src).rows; \ 66 (dst).cols = (src).cols; \ 67 (src).rows = 0; \ 68 (src).cols = 0; \ 69 } 70 #define MTX_CLONE(dst, src) { \ 71 (dst).data = (src).data.clone(); \ 72 (dst).rows = (src).rows; \ 73 (dst).cols = (src).cols; \ 74 } 75 public: 76 typedef _T _m_type; 77 typedef _vecII<_T> _col; 78 typedef _vecII<_col> _m; 79 // 80 /* 81 以指针为样式设计，普通指针所拥有的功能他都有 82 所以使用该迭代器只需把他想象成指针就行 83 /*/ 84 typedef _m_type* _col_iter; 85 class _row_iter 86 { 87 public: 88 _row_iter(int row, int col, _mtx* parent) 89 :r(row), c(col), parent(parent) 90 {} 91 _row_iter(_row_iter& it) = default; 92 _row_iter& operator=(const _row_iter& it){ 93 this->c = it.c; 94 this->r = it.r; 95 this->parent = it.parent; 96 } 97 98 private: 99 const int r; // 指定行 100 int c; 101 _mtx* parent; 102 103 public: // boring function 104 inline virtual _m_type& operator*(){ 105 return parent->data[c][r]; 106 } 107 inline virtual _row_iter& operator+=(int i){ 108 c += i; 109 110 return *this; 111 } 112 inline virtual _row_iter& operator-=(int i){ 113 *this += -1 * i; 114 115 return *this; 116 } 117 inline virtual _row_iter& operator++(){ 118 *this += 1; 119 return *this; 120 } 121 inline virtual _row_iter& operator--(){ 122 *this -= 1; 123 return *this; 124 } 125 inline virtual _row_iter operator++(int){ 126 _row_iter Tmp = *this; 127 *this += 1; 128 129 return Tmp; 130 } 131 inline virtual _row_iter operator--(int){ 132 _row_iter Tmp = *this; 133 *this -= 1; 134 135 return Tmp; 136 } 137 inline virtual _row_iter operator+(int i){ 138 _row_iter Tmp = *this; 139 Tmp += i; 140 141 return Tmp; 142 } 143 inline virtual _row_iter operator-(int i){ 144 _row_iter Tmp = *this; 145 Tmp -= i; 146 147 return Tmp; 148 } 149 inline virtual bool operator!=(const _row_iter& it)const{ 150 return (c != it.c); 151 } 152 inline virtual bool operator==(const _row_iter& it)const{ 153 return (c == it.c); 154 } 155 }; 156 157 explicit _mtx(int r = 0, int c = 0){ 158 if (!r || !c){ 159 rows = cols = 0; 160 } 161 MTX_SET(*this, r, c); 162 } 163 _mtx(const _mtx& m){ 164 *this = m; 165 } 166 _mtx(_mtx&& m){ 167 *this = m; 168 } 169 _mtx& operator=(const _mtx& m){ 170 MTX_CLONE(*this, m); 171 172 return *this; 173 } 174 _mtx& operator=(_mtx&& m){ 175 MTX_MOVE(*this, m); 176 177 return *this; 178 } 179 virtual ~_mtx(){ 180 rows = cols = 0; 181 } 182 183 public: // 矩阵基本操作 184 _mtx transposition(void)const{/** 转置 **/ 185 _mtx Tmp(cols, rows); 186 for (int i = 0; i < cols; ++i) 187 for (int j = 0; j < rows; ++j) 188 Tmp.getElement(i, j) = this->at(j, i); 189 190 return Tmp; 191 } 192 193 public: // read by iter (多线程由此实现) 194 _mtx sub_mtx(int r, int c, int w, int h){ 195 _mtx Tmp(w, h); 196 for (int i = 0; i < w; ++i) 197 for (int j = 0; j < h; ++j) 198 Tmp.getElement(j, i) = getElement(r + j, c + i); 199 200 return Tmp; 201 } 202 inline _mtx sub_col(int beg, int width){ 203 return sub_mtx(0, beg, width, rows); 204 } 205 inline _mtx sub_row(int beg, int width){ 206 return sub_mtx(beg, 0, cols, width); 207 } 208 inline _m_type& getElement(int r, int c)throw(){ 209 return data.dt_p[c][r]; 210 } 211 inline const _m_type& at(int r, int c)const{ 212 return data.dt_p[c][r]; 213 } 214 inline _col_iter col_of(int i)throw(){ 215 return _col_iter(data[i].begin()); 216 } 217 inline _col_iter col_end_of(int i)throw(){ 218 return _col_iter(data[i].begin() + rows); 219 } 220 inline _row_iter row_of(int i)throw(){ 221 _row_iter Tmp = { 222 i, 223 0, 224 this 225 }; 226 227 return Tmp; 228 } 229 inline _row_iter row_end_of(int i)throw(){ 230 _row_iter Tmp = { 231 i, 232 cols, 233 this 234 }; 235 236 return Tmp; 237 } 238 239 public: // write by it's function 240 /** 交换 **/ 241 void exchange_col(int p1, int p2)throw(){ 242 _col Tmp = C_MOVE(data[p1]); 243 data[p1] = C_MOVE(data[p2]); 244 data[p2] = C_MOVE(Tmp); 245 } 246 void exchange_row(int p1, int p2){ 247 _m_type Tmp; 248 _row_iter r1 = row_of(p1), 249 r2 = row_of(p2), 250 r1_end = row_end_of(p1); 251 while (r1 != r1_end){ 252 Tmp = C_MOVE(*r1); 253 *r1 = C_MOVE(*r2); 254 *r2 = C_MOVE(Tmp); 255 ++r1; 256 ++r2; 257 } 258 } 259 /** 添加 **/ 260 void append_cols(const _mtx& m){ 261 if (m.rows != this->rows) 262 throw(MTX_ERR__DIM_ERR(append_cols)); 263 data.append(m.data); 264 this->cols = data.size(); 265 } 266 void append_cols(_mtx&& m){ 267 if (m.rows != this->rows) 268 throw(MTX_ERR__DIM_ERR(append_cols)); 269 data.append(C_MOVE(m.data)); 270 this->cols = data.size(); 271 } 272 void append_rows(const _mtx& m){ 273 if (m.cols != this->cols) 274 throw(MTX_ERR__DIM_ERR(append_rows)); 275 for (int i = 0; i < cols; ++i) 276 data[i].append(m.data.at(i)); 277 this->rows = data[0].size(); 278 } 279 void append_rows(_mtx&& m){ 280 if (m.cols != this->cols) 281 throw(MTX_ERR__DIM_ERR(append_rows)); 282 for (int i = 0; i < cols; ++i) 283 data[i].append(C_MOVE(m.data[i])); 284 this->rows = data[0].size(); 285 } 286 void insert_cols(int before, const _mtx& m){ 287 if (m.rows != rows) 288 throw(MTX_ERR__DIM_ERR(insert_cols)); 289 data.insert(before, m.data); 290 this->cols = data.size(); 291 } 292 void insert_cols(int before, _mtx&& m){ 293 if (m.rows != rows) 294 throw(MTX_ERR__DIM_ERR(insert_cols)); 295 data.insert(before, C_MOVE(m.data)); 296 this->cols = data.size(); 297 } 298 void insert_rows(int before, const _mtx& m){ 299 if (m.cols != cols) 300 throw(MTX_ERR__DIM_ERR(insert_rows)); 301 for (int i = 0; i < cols; ++i) 302 data[i].insert(before, m.data[i]); 303 this->rows = data[0].size(); 304 } 305 void insert_rows(int before, _mtx&& m){ 306 if (m.cols != cols) 307 throw(MTX_ERR__DIM_ERR(insert_rows)); 308 for (int i = 0; i < cols; ++i) 309 data[i].insert(before, C_MOVE(m.data[i])); 310 this->rows = data[0].size(); 311 } 312 /** 移除 **/ 313 void remove_cols(int beg, int end){ 314 data.remove(beg, end); 315 this->cols = data.size(); 316 } 317 void remove_rows(int beg, int end){ 318 for (int i = 0; i < cols; ++i) 319 data[i].remove(beg, end); 320 this->rows = data[0].size(); 321 } 322 inline void remove_col(int pos){ 323 remove_cols(pos, pos + 1); 324 } 325 inline void remove_row(int pos){ 326 remove_rows(pos, pos + 1); 327 } 328 /** 替换 **/ 329 void replace_sub_mtx(int r, int c, const _mtx& m){ 330 for (int i = 0; i < m.cols; ++i) 331 for (int j = 0; j < m.rows; ++j) 332 getElement(r + j, c + i) = m.at(j, i); 333 } 334 void replace_sub_mtx(int r, int c, _mtx&& m){ 335 for (int i = 0; i < m.cols; ++i) 336 for (int j = 0; j < m.rows; ++j) 337 getElement(r + j, c + i) = C_MOVE(m.getElement(j, i)); 338 } 339 inline void replace_cols(int pos, const _mtx& m){ 340 replace_sub_mtx(0, pos, m); 341 } 342 inline void replace_cols(int pos, _mtx&& m){ 343 replace_sub_mtx(0, pos, C_MOVE(m)); 344 } 345 inline void replace_rows(int pos, const _mtx& m){ 346 replace_sub_mtx(pos, 0, m); 347 } 348 inline void replace_rows(int pos, _mtx&& m){ 349 replace_sub_mtx(pos, 0, C_MOVE(m)); 350 } 351 //// 此处添加代码 352 /** **/ 353 354 355 356 public: // 对象操作 357 void resize(int r, int c){ // 重置大小到fit 358 if (r == rows && c == cols) 359 return; 360 rows = r; 361 cols = c; 362 data.resize(c); 363 for (int i = 0; i < cols; ++i) 364 data[i].resize(r); 365 } 366 void reset(){ // 数据归零 367 for (int i = 0; i < cols; ++i){ 368 _col_iter it = col_of(i), 369 end = col_end_of(i); 370 while (it != end){ 371 *it = _m_type(); 372 ++it; 373 } 374 } 375 } 376 virtual _mtx clone()const{ // 暂为单线程 377 _mtx Tmp; 378 MTX_CLONE(Tmp, *this); 379 380 return Tmp; 381 } 382 virtual _mtx move(){ // 383 _mtx Tmp(0, 0); 384 MTX_MOVE(Tmp, *this); 385 386 return Tmp; 387 } 388 389 390 private: 391 void _extend_row(int extend_size){ 392 for (int i = 0; i < cols; ++i){ 393 data[i].resize(rows + extend_size); 394 } 395 rows += extend_size; 396 } 397 void _extend_col(int extend_size){ 398 _col Tmp(this->rows); 399 data.resize(data.size() + extend_size); 400 for (int i = 0; i < extend_size; ++i) 401 data[cols + i] = Tmp; 402 cols += extend_size; 403 } 404 405 public: 406 int rows; 407 int cols; 408 private: 409 _m data; 410 411 #undef MTX_SET // 412 #undef MTX_SET_VAL 413 #undef MTX_MOVE 414 #undef MTX_CLONE 415 }; 416 417 418 // 常用矩阵函数 419 Mtx zeros(int r, int c); 420 Mtx ones(int r, int c); 421 Mtx eye(int r, int c); 422 Mtx random(int r, int c); 423 Mtx diag(const Mtx::_col& col); 424 string mtx_display(Mtx& m); 425 long mtx_set_getKernalNumber(void); 426 427 class _martix :public _mtx<creal> 428 { 429 public: 430 typedef creal refer; 431 432 explicit _martix(int r = 0, int c = 0) 433 :_mtx(r, c) 434 {} 435 _martix(const _martix& mtx) 436 :_mtx(mtx) 437 {} 438 _martix(_martix&& mtx) 439 :_mtx(C_MOVE(mtx)) 440 {} 441 _martix(const _mtx& mtx) 442 :_mtx(mtx) 443 {} 444 _martix(_mtx&& mtx) 445 :_mtx(C_MOVE(mtx)) 446 {} 447 _martix& operator=(const _martix& mtx){ 448 _mtx::operator=(mtx); 449 450 return *this; 451 } 452 _martix& operator=(_martix&& mtx){ 453 _mtx::operator=(C_MOVE(mtx)); 454 455 return *this; 456 } 457 _martix& operator=(const _mtx& mtx){ 458 _mtx::operator=(mtx); 459 460 return *this; 461 } 462 _martix& operator=(_mtx&& mtx){ 463 _mtx::operator=(C_MOVE(mtx)); 464 465 return *this; 466 } 467 virtual ~_martix() 468 {} 469 470 public: // 矩阵操作 471 _martix operator+(_martix& mtx); 472 _martix operator-(_martix& mtx); 473 _martix operator*(_martix& mtx); 474 _martix operator*(refer val); 475 _martix operator^(int i); 476 477 private: 478 }; 479 480 481 #endif C_MTX
1 #include "c_mtx.h" 2 3 #pragma disable(line:65) 4 5 #pragma region Reg_Thread_Init 6 7 // 线程全局变量 8 _martix* x_m;//结果量 x 9 _martix* a_m;//运行量 A 10 _martix* b_m;//参数量 b Ax = b 11 _martix::refer val_m; 12 long thread_count; 13 pthread_t threads[4]; 14 pthread_rwlock_t rwlock; 15 16 int thread_global_init() 17 { 18 thread_count = mtx_set_getKernalNumber(); 19 pthread_rwlock_init(&rwlock, NULL); 20 21 return 0; 22 } 23 int thread_global_destroy() 24 { 25 thread_count = 0; 26 pthread_rwlock_destroy(&rwlock); 27 28 return 0; 29 } 30 typedef int(*fun)(); 31 #pragma data_seg(".CRT$XCU") // 线程全局初始化 32 fun start = { thread_global_init }; 33 #pragma data_seg() 34 #pragma data_seg(".CRT$XPU") // 线程全局销毁 35 fun destroy = { thread_global_destroy }; 36 #pragma data_seg() 37 38 #pragma endregion 以上初始化了线程变量 39 40 41 // 常用矩阵函数 42 Mtx zeros(int r, int c){ 43 return Mtx(r, c); 44 } 45 Mtx ones(int r, int c){ 46 Mtx m = zeros(r, c); 47 for (int i = 0; i < m.cols; ++i){ 48 Mtx::_col_iter it = m.col_of(i), 49 end = m.col_end_of(i); 50 while (it != end) 51 *it++ = 1; 52 } 53 54 return m; 55 } 56 Mtx eye(int r, int c){ 57 Mtx m = zeros(r, c); 58 for (int i = 0, rand = MIN(r, c); i < rand; ++i) 59 m.getElement(i, i) = 1; 60 61 return m; 62 } 63 Mtx random(int r, int c) 64 { 65 srand(time(NULL)); 66 Mtx Tmp(r, c); 67 double rd; 68 for (int c = 0; c < Tmp.cols; ++c) 69 for (int r = 0; r < Tmp.rows; ++r){ 70 rd = (double)std::rand() / Pi; 71 Tmp.getElement(r, c) = rd - (int)rd; 72 } 73 74 return Tmp; 75 } 76 Mtx diag(const Mtx::_col& col) 77 { 78 Mtx m(col.size(), col.size()); 79 for (int i = 0, len = col.size(); i < len; ++i) 80 m.getElement(i, i) = col.at(i); 81 82 return m; 83 } 84 string mtx_display(Mtx& m) 85 { 86 string str; 87 for (int i = 0; i < m.rows; ++i){ 88 Mtx::_row_iter it = m.row_of(i), 89 end = m.row_end_of(i); 90 while (it != end){ 91 str += cto_string(*it) + '\t'; 92 ++it; 93 } 94 str += '\n'; 95 } 96 97 return str; 98 } 99 long mtx_set_getKernalNumber(void) 100 { 101 #if defined(WIN32) 102 SYSTEM_INFO info; 103 GetSystemInfo(&info); 104 return info.dwNumberOfProcessors; 105 #elif defined(LINUX) || defined(SOLARIS) || defined(AIX) 106 return get_nprocs(); 107 #else 108 #error 109 #endif 110 } 111 112 113 114 // 线程函数 115 116 #pragma region Reg_Thread_Define 117 118 #define THREAD_SET_VAL(x_p, a_p, b_p) { \ 119 x_m = (x_p); \ 120 a_m = (a_p); \ 121 b_m = (b_p); \ 122 } 123 #define THREAD_CLEAN THREAD_SET_VAL(nullptr, nullptr, nullptr) 124 #define THREAD_COL_ALLOC(mtx_size, thread_index, thread_count)/*平均分配*/( (mtx_size)/(thread_count) + ( (mtx_size)%(thread_count) > (thread_index) ) ) // n/t+ (less>thread? 1 : 0) 125 #define THREAD_COL_BEGIN(mtx_size, thread_index, thread_count)/*起始点*/( (mtx_size)/(thread_count)*(thread_index) + ( (mtx_size)%(thread_count) > (thread_index) ? (thread_index) : (mtx_size)%(thread_count) ) ) // n/t*i + (less>thread? i : less) 126 #define THREAD_CALCU_BEGIN__FOR__(col_begin, col_alloc)/* 列循环--起始 */for(int i = (col_begin), len = (col_alloc)+(col_begin); i < len; ++i){ 127 #define THREAD_CALCU_END__FOR__/* 列循环--结束 */ } /**** 以 i 为外层索引 ****/ 128 129 #pragma endregion 定义了一些方便的线程计算操作 130 131 void *pth_add(void *idx) 132 {/*在外层调用处加写锁*/ 133 long thread_index = (long)idx; 134 int col_begin = THREAD_COL_BEGIN(x_m->cols, thread_index, thread_count); 135 int col_alloc = THREAD_COL_ALLOC(x_m->cols, thread_index, thread_count); 136 THREAD_CALCU_BEGIN__FOR__(col_begin, col_alloc) 137 for (int j = 0; j < x_m->rows; ++j) 138 { 139 x_m->getElement(j, i) = a_m->getElement(j, i) + b_m->getElement(j, i); 140 } 141 THREAD_CALCU_END__FOR__ 142 143 return NULL; 144 } 145 void *pth_minus(void *idx) 146 {/*在外层调用处加写锁*/ 147 long thread_index = (long)idx; 148 int col_begin = THREAD_COL_BEGIN(x_m->cols, thread_index, thread_count); 149 int col_alloc = THREAD_COL_ALLOC(x_m->cols, thread_index, thread_count); 150 THREAD_CALCU_BEGIN__FOR__(col_begin, col_alloc) 151 for (int j = 0; j < x_m->rows; ++j) 152 { 153 x_m->getElement(j, i) = a_m->getElement(j, i) - b_m->getElement(j, i); 154 } 155 THREAD_CALCU_END__FOR__ 156 157 return NULL; 158 } 159 void *pth_multiply(void *idx) 160 { 161 long thread_index = (long)idx; 162 int col_begin = THREAD_COL_BEGIN(x_m->cols, thread_index, thread_count); 163 int col_alloc = THREAD_COL_ALLOC(x_m->cols, thread_index, thread_count); 164 int a_col_relative_begin = THREAD_COL_BEGIN(a_m->cols, thread_index, thread_count) - col_begin; // 循环中与x列位置相关的起始位置 165 166 THREAD_CALCU_BEGIN__FOR__(col_begin, col_alloc) // 计算乘法 167 for (int j = 0; j < x_m->rows; ++j) 168 { 169 val_m = 0; 170 for (int k = 0; k < b_m->rows; ++k) 171 { 172 val_m += a_m->getElement(j, k) * b_m->getElement(k, i); 173 } 174 x_m->getElement(j, i) = val_m; 175 } 176 THREAD_CALCU_END__FOR__ 177 178 return NULL; 179 } 180 void *pth_multiply_dot(void *idx) 181 { 182 long thread_index = (long)idx; 183 int col_begin = THREAD_COL_BEGIN(x_m->cols, thread_index, thread_count); 184 int col_alloc = THREAD_COL_ALLOC(x_m->cols, thread_index, thread_count); 185 186 _martix::_col_iter it, end; 187 THREAD_CALCU_BEGIN__FOR__(col_begin, col_alloc) 188 it = x_m->col_of(i), 189 end = x_m->col_end_of(i); 190 while (it != end){ 191 *it *= val_m; 192 ++it; 193 } 194 THREAD_CALCU_END__FOR__ 195 196 return NULL; 197 } 198 199 200 201 202 203 204 #pragma region Reg_Mtx_Define 205 206 // 用户接口 207 void m_add(_martix* x, _martix* a, _martix* b) 208 { 209 THREAD_SET_VAL(x, a, b); 210 pthread_rwlock_wrlock(&rwlock); 211 for (int i = 0; i < thread_count; ++i) 212 pthread_create(&threads[i], NULL, pth_add, (void*)i); 213 for (int i = 0; i < thread_count; ++i) 214 pthread_join(threads[i], NULL); 215 pthread_rwlock_unlock(&rwlock); 216 THREAD_CLEAN; 217 } 218 void m_minus(_martix* x, _martix* a, _martix* b) 219 { 220 THREAD_SET_VAL(x, a, b); 221 pthread_rwlock_wrlock(&rwlock); 222 for (int i = 0; i < thread_count; ++i) 223 pthread_create(&threads[i], NULL, pth_minus, (void*)i); 224 for (int i = 0; i < thread_count; ++i) 225 pthread_join(threads[i], NULL); 226 pthread_rwlock_unlock(&rwlock); 227 THREAD_CLEAN; 228 } 229 void m_multiply(_martix* x, _martix* a, _martix* b) 230 { 231 THREAD_SET_VAL(x, a, b); 232 pthread_rwlock_wrlock(&rwlock); 233 for (int i = 0; i < thread_count; ++i) 234 pthread_create(&threads[i], NULL, pth_multiply, (void*)i); 235 for (int i = 0; i < thread_count; ++i) 236 pthread_join(threads[i], NULL); 237 pthread_rwlock_unlock(&rwlock); 238 THREAD_CLEAN; 239 } 240 void m_multiply_dot(_martix* x, _martix::refer v) 241 { 242 THREAD_SET_VAL(x, NULL, NULL); 243 val_m = v; 244 pthread_rwlock_wrlock(&rwlock); 245 for (long i = 0; i < thread_count; ++i) 246 pthread_create(&threads[i], NULL, pth_multiply_dot, (void*)i); 247 for (long i = 0; i < thread_count; ++i) 248 pthread_join(threads[i], NULL); 249 pthread_rwlock_unlock(&rwlock); 250 } 251 252 253 _martix _martix::operator+(_martix& mtx){ 254 if (mtx.rows != rows || mtx.cols != cols) 255 throw(MTX_ERR__DIM_ERR(operator+=)); 256 257 _martix Tmp(rows, cols); 258 m_add(&Tmp, this, &mtx); 259 260 return Tmp; 261 } 262 _martix _martix::operator-(_martix& mtx){ 263 if (mtx.rows != rows || mtx.cols != cols) 264 throw(MTX_ERR__DIM_ERR(operator+=)); 265 266 _martix Tmp(rows, cols); 267 m_minus(&Tmp, this, &mtx); 268 269 return Tmp; 270 } 271 _martix _martix::operator*(_martix& mtx){ 272 if (cols != mtx.rows) 273 throw(MTX_ERR__DIM_ERR(operator*)); 274 275 _martix Tmp(rows, mtx.cols); 276 m_multiply(&Tmp, this, &mtx); 277 278 return Tmp; 279 } 280 _martix _martix::operator*(refer val){ 281 _martix Tmp = this->clone(); 282 m_multiply_dot(&Tmp, val); 283 284 return Tmp; 285 } 286 _martix _martix::operator^(int i){ 287 if (cols != rows) 288 throw(MTX_ERR__DIM_ERR(operator^)); 289 290 _martix Tmp = this->clone(); 291 for (int k = 0; k < i - 1; ++k) 292 Tmp = Tmp*(*this); 293 294 return Tmp; 295 } 296 297 298 #pragma endregion 这一段是数值矩阵运算的定义
接下来我可能会写一些多线程实现的数据结构，学学pthread和OpenMP
the end
----------------------2018-05-15-----------------------
说明一下效率，vs2013 release模式下500*500用时63mms，debug模式下vs忽略inline调用，因此性能差异巨大
每一次的函数调用包括了：
参数的拷贝（写入堆栈，从堆栈读取，函数返回时还要从堆栈恢复
环境的保存与恢复，包括寄存器，返回地址等等，而矩阵相乘每次都需要用_vecII的接口读取，性能差异10倍就很正常了。
验证这个结论很简单：在debug mode下，加入/Ob1编译选项（这个选项打开了inline功能），强迫编译器把inline函数作为inline处理，此时debugging information也不能选program database for edit and continue, 从而也无法进入相关的函数了。这时候运行，我们会发现速度和release mode一样了。参考大佬的文章。
说明一下pragma data_seg(".CRT$XCU")，因为里面用到了，用于在主函数前调用初始化线程，有点像全局的构造和析构函数。
运行一个程序，首先调用的部分必然是与相应操作系统和程序所需外部支持有关的代码，外部项配置好之后才会进入我们所写的main函数中；编译时会生成两个表，一个在程序跑之前调用，用于初始化，另一个在程序结束时调用，用于安排程序后事。万能的data_seg(".CRT$XCU")指令帮助我们把某些需要提前调用的函数插入表中。详细内容这里有一篇大佬的文章。
然后整体做了一些小改动，加了一个线程安全的链表，准备等多写一点的时候在一次性加进去
觉得好记得分享，想借用给我点个赞就行，哈哈哈，祝你快乐！
