数学库（持续更新中 18-05-06）
2024-05-01 05:07:33
命名约定
/** 命名类型统一
**
**
**
**
**/
#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, // SymbolOperate_Type_Num, // NumberOperate_Type_Fun, // FunctionOperate_Type_Op, // OpearteOperate_Type_Bu, // Bound
};enum CVar_Type
{Var_Type_Unknow = 0,Var_Type_Num,Var_Type_Fun,Var_Type_Handle,Var_Type_Array
};#endif // CPLX_NAMEDcplx_named.h
基本库
基本库的大部分是开学初写的，现在看看觉得当时的代码好幼稚，可能这就是成长吧（被自己的厚脸皮逗笑）
因为存在一些项目已经用到了里面的一些内容，所以在此就不对其进行更新了，术语向后兼容
新添加了一个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_COREcplx_core.h
基本矩阵存储及数值运算矩阵
实现了矩阵存储模板，一些方便的函数都已经定义好了，还算完善
实现了数值矩阵基本操作，以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_MTXc_mtx.h
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 这一段是数值矩阵运算的定义c_mtx.cpp
接下来我可能会写一些多线程实现的数据结构，学学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")指令帮助我们把某些需要提前调用的函数插入表中。详细内容这里有一篇大佬的文章。
然后整体做了一些小改动，加了一个线程安全的链表，准备等多写一点的时候在一次性加进去
觉得好记得分享，想借用给我点个赞就行，哈哈哈，祝你快乐！
转载于:https://www.cnblogs.com/dusty-cjh/p/dusty-math-lab.html
数学库（持续更新中 18-05-06）相关推荐
- python基础和第三方库 笔记（python基础完结包括高级用法，第三方库持续更新中...）
python基础 注:本笔记面向有一定基础的人 本笔记是本人快速复习python过程中记录的,不适合零基础的人学习python的主工具,可以作为辅工具,本笔记记录了入门阶段常用操作,如有错误的地方,希 ...
- 嵌入式相关开源项目、库、资料------持续更新中
学习初期最难找的就是找学习资料了,本贴精心汇总了一些嵌入式相关资源,包括但不限于编程语言.单片机.开源项目.物联网.操作系统.Linux.计算机等资源,并且在不断地更新中,致力于打造全网最全的嵌入式资 ...
- STM32 之十五 奇怪问题处理及驱动库 BUG 汇总（持续更新中）
在使用 STM32 的 MCU 开发过程中,难免遇到各种各样的奇葩问题.或许是开发环境的问题,或许是 MCU 使用的问题,也或许是驱动库的 BUG 等等.这些问题可能不局限于某一种具体型号的 MC ...
- 最好的Vue组件库之Vuetify的入坑指南（持续更新中）
目录 安装Vuetify 文档结构 快速入门 特性 样式和动画 首先先声明,个人不是什么很牛逼的大佬,只是想向那些想入坑Vuetify的前端新手或者嫌文档太长不知如何入手的人提供一些浅显的建议而已,能 ...
- 【Spring Boot官方文档原文理解翻译-持续更新中】
[Spring Boot官方文档原文理解翻译-持续更新中] 文章目录 [Spring Boot官方文档原文理解翻译-持续更新中] Chapter 4. Getting Started 4.1. Int ...
- 面试1：Java、微服务、架构常见面试题（持续更新中）
Java.微服务.架构常见面试题(持续更新中) 文章目录 Java.微服务.架构常见面试题(持续更新中) ==**Java**== 1.Java概述 (1)JVM.JRE和JDK (2)Java特点 ...
- 我学习 Java 的历程和体会（写给新手看，欢迎老司机批评和建议，持续更新中）
我学习 Java 的历程和体会(写给新手看,欢迎老司机批评和建议,持续更新中) 最初写这篇文章的时候,是在今年的 9 月中旬.今天,我想再写写这将近两个多月以来的感受. 在今年的 10 月我来到北京求 ...
- linux运维常用命令一句话（持续更新中）
整理收集一些Linux运维管理.系统管理的常用命令,太多了记不住,只能记录下来方便日后查看.也可以和大家分享.如果你有好的一句话命令也贴出来吧.本文持续更新中. 1.linux启动过程 开启电源 -- ...
- db2dual_DB2常用SQL的写法(持续更新中...)
DB2常用SQL的写法(持续更新中...) -- Author: lavasoft -- Date : 2006-12-14 -- 创建一个自定义单值类型 create distinct typ ...
- 前端知识点总结——JS高级（持续更新中）
前端知识点总结--JS高级(持续更新中) 1.字符串 什么是: 连续存储多个字符的字符数组 相同: 1. 下标 2. .length 3. 遍历 4. 选取: slice(starti[, endi] ...
最新文章
- 限制 TEdit 只能接收 数字 输入...
- Android----获取activity上所有的控件
- FastJson的使用方法总结
- Linux移植随笔：终于解决Tslib的问题了
- 非平衡电桥电阻计算_直流双臂电桥使用方法，统统告诉你
- 算法导论的道与术、工程师思维奠定能走多远
- PowerDesigner 数据字典模板
- 数据库左连接、右连接
- 数学实验基于matlab软件,数学实验：基于MATLAB软件
- 当贝投影F5发布：3000ANSI流明亮度，6S疾速开机，16ms超低延迟
- ural 1818 Fair Fishermen
- DataMatrix 编码生成和译码原理即方法
- 2022茶艺师（中级）考试模拟100题及模拟考试
- 今后，若你的公众号还按老方式发广告，罚款高的可达百万。有公众号已经被惩处了。
- html标题字体重叠,如何在HTML和CSS的段落标题上方显示堆叠的字体超棒图标？
- 中科大软院数据挖掘(LJL)考试回忆+课堂笔记+复习资料
- 存储系统测试工具学习之vdbench
- “数据结构基础”系列网络课程主页
- 使用jquery判断是否为数字
- hpuoj 1715: 感恩节KK专场——2015年的第一场雪 （暴力）
热门文章
- Python--Redis实战：第五章：使用Redis构建支持程序：第2节：计数器和统计数据...
- 小白理解k-svd算法
- RC522RFID读卡寻卡失败原因
- 18 | 安全标准和框架：怎样依“葫芦”画出好“瓢”？
- html 支持ssi,nginx下开启ssi支持
- python中uppercase是什么意思_Python string.ascii_uppercase方法代码示例
- boost::stacktrace::stacktrace相关的测试程序
- python calu_Python中*args，**kwargs的使用
- 2021年12月27日多载波NOMA论文阅读
- 用这个工具，让人抓狂的领导驾驶舱报表五步搞定