C++ mpi module for stochmagnet_main Package
functions_array.h
1 #ifndef FUNCTIONS_ARRAY_H
2 #define FUNCTIONS_ARRAY_H
3 
4 //types of the code
5 #include "types.h"
6 #include "functions_type.h"
7 
8 //standart library
9 #include <stdlib.h>
10 
11 //for type name
12 #include <cxxabi.h>
13 
14 //include the Template condition
15 #include<type_traits>
16 
17 //regex header for string manipulation
18 #include <regex>
19 
20 //template type verification header
21 #include<type_traits>
22 #include<concepts>
23 
24 //math header isnan
25 #include <math.h>
26 
27 //numeric function utilities header
28 #include "functions_numeric.h"
29 #include "functions_string.h"
30 
31 #ifdef DEBUG
32 //define the assert method
33 #include<cassert>
34 #endif
35 
36 namespace functions_array {
37 
38 
39  //SET method
42  template<typename T,size_t D> requires (functions_type::isArithmeticType<T>)
43  inline void reset(std::array<T,D>& a) {
44  memset(a.data(),0,sizeof(T)*D);
45  }
48  template<typename T> requires (functions_type::isArithmeticType<T>)
49  inline void reset(std::valarray<T>& a) {
50  if (a.size()>0) memset(&a[0],0,sizeof(T)*a.size());
51  }
52 
53  //arithmetic vector & array methods
54  //==================================
55 
56  //mixed product
57  //---------------
58 
61  template<typename T,size_t D,size_t L> requires ((L==D) && functions_type::isArithmeticType<T>)
62  inline T levelMixedProduct(const T* A,const T* B,const T* C) {
63  return (*C)*(A[1-D]*B[-1]-A[-1]*B[1-D]);
64  }
67  template<typename T,size_t D,size_t L> requires ((L>1) && (L<D) && functions_type::isArithmeticType<T>)
68  inline T levelMixedProduct(const T* A,const T* B,const T* C) {
69  return (*C)*(A[1]*B[-1]-A[-1]*B[1])+levelMixedProduct<T,D,L+1>(A+1,B+1,C+1);
70  }
73  template<typename T,size_t D,size_t L> requires ((L==1) && functions_type::isArithmeticType<T>)
74  inline T levelMixedProduct(const T* A,const T* B,const T* C) {
75  return (*C)*(A[1]*B[D-1]-A[D-1]*B[1])+levelMixedProduct<T,D,2>(A+1,B+1,C+1);
76  }
77 
78 
79 
80 
83  template<typename T,size_t D> requires (functions_type::isArithmeticType<T>)
84  inline T mixedProduct(const T* A,const T* B,const T* C) {
85  return levelMixedProduct<T,D,1>(A,B,C);
86  }
89  template<typename T,size_t D> requires (functions_type::isArithmeticType<T>)
90  inline T mixedProduct(const std::array<T,D>& A,const std::array<T,D>& B,const std::array<T,D>& C) {
91  return levelMixedProduct<T,D,1>(A.data(),B.data(),C.data());
92  }
93 
94  //vector product
95  //---------------
96 
99  template<typename T,size_t D,size_t L> requires ((L==D) && functions_type::isArithmeticType<T>)
100  inline void levelVectorProduct(const T* A,const T* B,T* R) {
101  (*R)+=(A[1-D]*B[-1]-A[-1]*B[1-D]);
102  }
105  template<typename T,size_t D,size_t L> requires ((L>1) && (L<D) && functions_type::isArithmeticType<T>)
106  inline void levelVectorProduct(const T* A,const T* B,T* R) {
107  (*R)+=(A[1]*B[-1]-A[-1]*B[1]);
108  levelVectorProduct<T,D,L+1>(A+1,B+1,R+1);
109  }
112  template<typename T,size_t D,size_t L> requires ((L==1) && functions_type::isArithmeticType<T>)
113  inline void levelVectorProduct(const T* A,const T* B,T* R) {
114  (*R)+=(A[1]*B[D-1]-A[D-1]*B[1]);
115  levelVectorProduct<T,D,2>(A+1,B+1,R+1);
116  }
117 
118 
119 
122  template<typename T,size_t D> requires (functions_type::isArithmeticType<T>)
123  inline void vectorProduct(const T* A,const T* B, T* R) {
124  ASSERT_IN((R!=A) && (R!=B));
125  levelVectorProduct<T,D,1>(A,B,R);
126  }
129  template<typename T,size_t D> requires (functions_type::isArithmeticType<T>)
130  inline void vectorProduct(const T* A,const T* B,const tBoolean& alpha,T* R) {
131  if (alpha==0) memset(R,0,sizeof(T)*D);
132  vectorProduct<T,D,1>(A,B,R);
133  }
134 
137  template<typename T,size_t D> requires (functions_type::isArithmeticType<T>)
138  inline void vectorProduct(const std::array<T,D>& A,const std::array<T,D>& B,const tBoolean& alpha, std::array<T,D>& R) {
139  vectorProduct<T,D,1>(A.data(),B.data(),alpha,R.data());
140  }
141 
142 
143  //scalar product
144  //---------------
148  template<typename T,typename Q,size_t L,size_t D> requires ( (L==D) && functions_type::isArithmeticType<T>&& functions_type::isArithmeticType<Q>)
149  inline T levelScalarProduct(const T * M,const Q *H) {
150  return (*M)*(*H);
151  }
152 
156  template<typename T,typename Q,size_t L,size_t D> requires ( (L<D) && functions_type::isArithmeticType<T> && functions_type::isArithmeticType<Q>)
157  inline T levelScalarProduct(const T * M,const Q *H) {
158  return (*M)*(*H)+levelScalarProduct<T,Q,L+1,D>(M+1,H+1);
159  }
160 
161 
162 
166  template<typename T,typename Q,size_t D> requires (functions_type::isArithmeticType<T> && functions_type::isArithmeticType<Q>)
167  inline T scalarProduct(const T* M,const Q* H) {
168  return levelScalarProduct<T,Q,1,D>(M,H);
169  }
173  template<typename T,size_t D> requires functions_type::isArithmeticType<T>
174  inline T scalarProduct(const T* M,const T* H) {
175  return levelScalarProduct<T,T,1,D>(M,H);
176  }
177 
181  template<typename T,typename Q,size_t D> requires (functions_type::isArithmeticType<T> && functions_type::isArithmeticType<Q>)
182  inline T scalarProduct(const std::array<T,D>& M,const std::array<Q,D>& H) {
183  return levelScalarProduct<T,Q,1,D>(M.data(),H.data());
184  }
188  template<typename T,size_t D> requires functions_type::isArithmeticType<T>
189  inline T scalarProduct(const std::array<T,D>& M,const std::array<T,D>& H) {
190  return levelScalarProduct<T,T,1,D>(M.data(),H.data());
191  }
192 
193  //norm
194  //=====
195 
199  template<typename T,size_t D> requires functions_type::isArithmeticType<T>
200  inline T norm2(const std::array<T,D>& M) {
201  return levelScalarProduct<T,T,1,D>(M.data(),M.data());
202  }
206  template<typename T,size_t D> requires functions_type::isArithmeticType<T>
207  inline tReal norm(const std::array<T,D>& M) {
208  return sqrt(norm2(M));
209  }
210 
214  template<typename T, size_t D> requires functions_type::isArithmeticType<T>
215  inline tReal norm2(const T* M) {
216  return levelScalarProduct<T,T,1,D>(M,M);
217  }
221  template<typename T,size_t D> requires functions_type::isArithmeticType<T>
222  inline tReal norm(const T* M) {
223  return sqrt(norm2(M));
224  }
225 
226  //L2 Distance
227  //-------------
231  template<typename T,size_t L,size_t D> requires ( (L==D) && functions_type::isArithmeticType<T>)
232  inline T levelL2Distance2(const T * M,const T *H) {
233  return ((*M)-(*H))*((*M)-(*H));
234  }
235 
239  template<typename T,size_t L,size_t D> requires ( (L<D) && functions_type::isArithmeticType<T>)
240  inline T levelL2Distance2(const T * M,const T *H) {
241  return ((*M)-(*H))*((*M)-(*H))+levelL2Distance2<T,L+1,D>(M+1,H+1);
242  }
248  template<typename T,size_t D>
249  inline T L2Distance2(const std::array<T,D>& M,const std::array<T,D>& H) {
250  return levelL2Distance2<T,1,D>(M.data(),H.data());
251  }
257  template<typename T,size_t D>
258  inline T L2Distance2(const T* M,const T* H) {
259  return levelL2Distance2<T,1,D>(M,H);
260  }
261 
267  template<typename T,size_t D>
268  inline T L2Distance(const std::array<T,D>& M,const std::array<T,D>& H) {
269  return sqrt(L2Distance2(M,H));
270  }
276  template<typename T,size_t D>
277  inline T L2Distance(const T* M,const T* H) {
278  return sqrt(L2Distance2<T,D>(M,H));
279  }
280 
281  // product of elements
282  // ------------------
286  template<typename S,typename T,size_t L, size_t D> requires ((L==D) && functions_type::isArithmeticType<T> && functions_type::isArithmeticType<S>)
287  inline S levelProduct(const T* a) {
288  return *a;
289  }
290 
294  template<typename S,typename T,size_t L, size_t D> requires ((L<D) && functions_type::isArithmeticType<T> && functions_type::isArithmeticType<S>)
295  inline S levelProduct(const T* a) {
296  return ((S)(*a))*levelProduct<S,T,L+1,D>(a+1);
297  }
301  template<typename S,typename T,size_t D> requires (functions_type::isArithmeticType<T> && functions_type::isArithmeticType<S>)
302  inline S product(const std::array<T,D>& a) {
303  return levelProduct<S,T,1,D>(a.data());
304  }
308  template<typename T,size_t D> requires (functions_type::isArithmeticType<T>)
309  inline T product(const std::array<T,D>& a) {
310  return levelProduct<T,T,1,D>(a.data());
311  }
315  template<typename T> requires (functions_type::isArithmeticType<T>)
316  inline T product(const std::vector<T>& as) {
317  T p=1;
318  for(const auto& a:as) {
319  p*=a;
320  }
321  return p;
322  }
323 
324 
328  template<typename S,typename T,size_t D> requires (functions_type::isArithmeticType<T> && functions_type::isArithmeticType<S>)
329  inline S product(const T*a) {
330  return levelProduct<S,T,1,D>(a);
331  }
332 
333 
334  //reference to min value of an array
335  //-----------------------------------
340  template<typename T,size_t L, size_t D> requires ((L==D) && functions_type::isArithmeticType<T>)
341  inline const T& levelMin(const T* a) {
342  return *a;
343  }
348  template<typename T,size_t L, size_t D> requires ((L<D) && functions_type::isArithmeticType<T>)
349  inline const T& levelMin(const T* a) {
350  return functions_numeric::min(levelMin<T,L+1,D>(a+1),*a);
351  }
352 
357  template<typename T,size_t D> requires functions_type::isArithmeticType<T>
358  inline const T& min(const std::array<T,D>& a) {
359  return levelMin<T,1,D>(a.data());
360  }
365  template<typename T,size_t D> requires functions_type::isArithmeticType<T>
366  inline const T& min(const T* a) {
367  return levelMin<T,1,D>(a);
368  }
369 
370  //reference to max value of an array
371  //-----------------------------------
376  template<typename T,size_t L, size_t D> requires ((L==D) && functions_type::isArithmeticType<T>)
377  inline const T& levelMax(const T* a) {
378  return *a;
379  }
380 
385  template<typename T,size_t L, size_t D> requires ((L<D) && functions_type::isArithmeticType<T>)
386  inline const T& levelMax(const T* a) {
387  return functions_numeric::max(levelMax<T,L+1,D>(a+1),*a);
388  }
389 
394  template<typename T,size_t D> requires functions_type::isArithmeticType<T>
395  inline const T& max(const std::array<T,D>& a) {
396  return levelMax<T,1,D>(a.data());
397  }
402  template<typename T,size_t D> requires functions_type::isArithmeticType<T>
403  inline const T& max(const T* a) {
404  return levelMax<T,1,D>(a);
405  }
406 
407 
408  //a:=max(a,b)
409  //============
415  template<typename T,size_t L, size_t D> requires ((L==D) && functions_type::isArithmeticType<T>)
416  inline void levelMaxArray(const T* al,const T* bl,T* rl) {
417  (*rl)=std::max(*al,*bl);
418  }
424  template<typename T,size_t L, size_t D> requires ((L<D) && functions_type::isArithmeticType<T>)
425  inline void levelMaxArray(const T* al,const T* bl,T* rl) {
426  (*rl)=std::max(*al,*bl);
427  levelMaxArray<T,L+1,D>(al+1,bl+1,rl+1);
428  }
434  template<typename T,size_t D>
435  inline std::array<T,D>& max(std::array<T,D>& a,const std::array<T,D>& b) {
436  levelMaxArray<T,1,D>(a.data(),b.data(),a.data());
437  return a;
438  }
444  template<typename T>
445  inline std::valarray<T>& max(std::valarray<T>& a,const std::valarray<T>& b) {
446  const auto *ib=&b[0];
447  for(auto& ak:a) {
448  ak=(ak>(*ib))?ak:(*ib);
449  ib++;
450  }
451  return a;
452  }
453 
454  //a:=min(a,b)
455  //===========
461  template<typename T,size_t L, size_t D> requires ((L==D) && functions_type::isArithmeticType<T>)
462  inline void levelMinArray(const T* al,const T* bl,T* rl) {
463  (*rl)=std::min(*al,*bl);
464  }
465 
471  template<typename T,size_t L, size_t D> requires ((L<D) && functions_type::isArithmeticType<T>)
472  inline void levelMinArray(const T* al,const T* bl,T* rl) {
473  (*rl)=std::min(*al,*bl);
474  levelMinArray<T,L+1,D>(al+1,bl+1,rl+1);
475  }
476 
482  template<typename T,size_t D>
483  inline std::array<T,D>& min(std::array<T,D>& a,const std::array<T,D>& b) {
484  levelMinArray<T,1,D>(a.data(),b.data(),a.data());
485  return a;
486  }
492  template<typename T>
493  inline std::valarray<T>& min(std::valarray<T>& a,const std::valarray<T>& b) {
494  const auto *ib=&b[0];
495  for(auto& ak:a) {
496  ak=(ak<(*ib))?ak:(*ib);
497  ib++;
498  }
499  return a;
500  }
501  //a:=b
502  //=====
503 
508  template<typename T,typename Q,size_t L, size_t D> requires (L==D)
509  inline void levelCopy(T* al,const Q* bl) {
510  (*al)=(T)(*bl);
511  }
512 
517  template<typename T,typename Q,size_t L, size_t D> requires (L<D)
518  inline void levelCopy(T* al,const Q* bl) {
519  (*al)=(T) (*bl);
520  levelCopy<T,Q,L+1,D>(al+1,bl+1);
521  }
522 
528  template<typename T,typename Q,size_t D>
529  inline std::array<T,D>& copy(std::array<T,D>& a,const std::array<Q,D>& b) {
530  levelCopy<T,Q,1,D>(a.data(),b.data());
531  return a;
532  }
533 
534  //a:=a+b
535  //======
541  template<typename T,typename Q,size_t L, size_t D> requires ((L==D) && functions_type::isArithmeticType<T> && functions_type::isArithmeticType<Q>)
542  inline void levelAdd(T* al,const Q* bl) {
543  (*al)+=(T) (*bl);
544  }
545 
551  template<typename T,typename Q,size_t L, size_t D> requires ((L<D) && functions_type::isArithmeticType<T> && functions_type::isArithmeticType<Q>)
552  inline void levelAdd(T* al,const Q* bl) {
553  (*al)+=(T) (*bl);
554  levelAdd<T,Q,L+1,D>(al+1,bl+1);
555  }
556 
562  template<typename T,typename Q,size_t D> requires ( functions_type::isArithmeticType<T> && functions_type::isArithmeticType<Q>)
563  inline std::array<T,D>& add(std::array<T,D>& a,const std::array<Q,D>& b) {
564  levelAdd<T,Q,1,D>(a.data(),b.data());
565  return a;
566  }
567 
568  //a:=a-b
569  //======
575  template<typename T,typename Q,size_t L, size_t D> requires ((L==D) && functions_type::isArithmeticType<T> && functions_type::isArithmeticType<Q>)
576  inline void levelSub(T* al,const Q* bl) {
577  (*al)-=(T) (*bl);
578  }
579 
585  template<typename T,typename Q,size_t L, size_t D> requires ((L<D) && functions_type::isArithmeticType<T> && functions_type::isArithmeticType<Q>)
586  inline void levelSub(T* al,const Q* bl) {
587  (*al)-=(T) (*bl);
588  levelSub<T,Q,L+1,D>(al+1,bl+1);
589  }
595  template<typename T,typename Q,size_t D> requires ( functions_type::isArithmeticType<Q> && functions_type::isArithmeticType<T>)
596  inline std::array<T,D>& sub(std::array<T,D>& a,const std::array<Q,D>& b) {
597  levelSub<T,Q,1,D>(a.data(),b.data());
598  return a;
599  }
605  template<typename T,typename Q> requires ( functions_type::isArithmeticType<Q> && functions_type::isArithmeticType<T>)
606  inline std::valarray<T>& sub(std::valarray<T>& a,const std::valarray<Q>& b) {
607  const Q* ib=&b[0];
608  for(auto& ak:a) {
609  ak-=(T) (*ib);
610  ib++;
611  }
612  return a;
613  }
614 
615  //y:=a.x+y
616  //=============
623  template<typename T,typename Q,typename R,size_t L, size_t D> requires ((L==D) && functions_type::isArithmeticType<T>)
624  inline void levelAXPY(const R& a,const Q* xl,T* yl) {
625  (*yl)+=(T) (((Q)a)*(*xl));
626  }
627 
634  template<typename T,typename Q,typename R,size_t L, size_t D> requires ((L<D) && functions_type::isArithmeticType<T>)
635  inline void levelAXPY(const R& a,const Q* xl,T* yl) {
636  (*yl)+=(T) (((Q)a)*(*xl));
637  levelAXPY<T,Q,R,L+1,D>(a,xl+1,yl+1);
638  }
639 
646  template<typename T,typename Q,typename R,size_t D>
647  inline std::array<T,D>& axpy(const R& alpha,const std::array<Q,D>& x,std::array<T,D>& y) {
648  levelAXPY<T,Q,R,1,D>(alpha,x.data(),y.data());
649  return y;
650  }
651  //y[i]=a[i].x[i]+y[i]
652  //===================
659  template<typename T,typename Q,typename R,size_t L, size_t D> requires (L==D)
660  inline void levelAXPY(const R* al,const Q* xl,T* yl) {
661  (*yl)+=(T) (((R)(*al))*(*xl));
662  }
663 
670  template<typename T,typename Q,typename R,size_t L, size_t D> requires (L<D)
671  inline void levelAXPY(const R* al,const Q* xl,T* yl) {
672  (*yl)+=(T) (((R)(*al))*(*xl));
673  levelAXPY<T,Q,R,L+1,D>(al+1,xl+1,yl+1);
674  }
675 
682  template<typename T,typename Q,typename R,size_t D> requires ( functions_type::isArithmeticType<Q> && functions_type::isArithmeticType<T>)
683  inline std::array<T,D>& axpy(const std::array<R,D>& alpha,const std::array<Q,D>& x,std::array<T,D>& y) {
684  levelAXPY<T,Q,R,1,D>(alpha.data(),x.data(),y.data());
685  return y;
686  }
687 
688 
689  //TO STRING
690  //=========
691 
692 
695  template<typename Q,typename T,size_t D>
696  inline tString toString(const std::array<T,D>& a) {
697  std::stringstream ss;
698  ss<<"{";
699  for(const auto & ai : a) {
700  ss<<std::setprecision(12)<<((Q)ai)<<",";
701  }
702  if (a.size()>0)ss.seekp(-1, std::ios_base::end);
703  ss<<"}";
704  return ss.str();
705  }
706 
707 
710  template<typename T,size_t D>
711  inline tString toString(const std::array<T,D>& a) {
712  std::stringstream ss;
713  ss<<"{";
714  for(const auto & ai : a) {
715  ss<<std::setprecision(12)<<ai<<",";
716  }
717  if (a.size()>0)ss.seekp(-1, std::ios_base::end);
718  ss<<"}";
719  return ss.str();
720  }
721 
724  template<typename Q,typename T>
725  inline tString toString(const size_t& asize,const T* a,const tIndex& nPrintedValues) {
726  std::stringstream ss;
727  //end iterator on a
728  const T* ea=a;ea+=asize;
729  //begin iterator on a
730  const T* ia=a;
731  ss<<"{";
732  if (nPrintedValues>=asize) {
733  while (ia!=ea) {
734  ss<<std::setprecision(12)<<((Q)(*ia))<<",";
735  //next iterator
736  ia++;
737  }
738  }else {
739  tIndex index;
740  for (tInteger k=0;k<=nPrintedValues;k++) {
741  index=((asize-1)*k)/nPrintedValues;
742  ss<<index<<":"<<std::setprecision(12)<<((Q)a[index])<<",";
743  }
744  }
745  if (asize>0)ss.seekp(-1, std::ios_base::end);
746  ss<<"}";
747  return ss.str();
748  }
751  template<typename Q,typename T>
752  inline tString toString(const size_t& asize,const T* a) {
753  return toString<Q,T>(asize,a,asize);
754  }
757  template<typename T>
758  inline tString toString(const size_t& asize,const T* a) {
759  return toString<T,T>(asize,a,asize);
760  }
761 
762  template<typename S,typename T>
763  inline tString toString(const std::valarray<T>& a,const tIndex& nPrintedValues) {
764  return toString<S,T>(a.size(),&a[0],nPrintedValues);
765  }
766  template<typename T>
767  inline tString toString(const std::valarray<T>& a,const tIndex& nPrintedValues) {
768  return toString<T,T>(a.size(),&a[0],nPrintedValues);
769  }
770 
773  template<typename T>
774  inline tString toString(const std::valarray<T>& a) {
775  return toString<T,T>(a.size(),&a[0],a.size());
776  }
779  template<typename T>
780  inline tString toString(const std::vector<T>& a) {
781  std::stringstream ss;
782  ss<<"{";
783  for(const auto & ai : a) {
784  ss<<ai<<",";
785  }
786  if (a.size()>0)ss.seekp(-1, std::ios_base::end);
787  ss<<"}";
788  return ss.str();
789  }
792  template<typename Kd,typename Vd,typename K,typename V>
793  inline tString toString(const std::map<K,V>& a) {
794  std::stringstream ss;
795  ss<<"{";
796  for(const auto & ai : a) {
797  ss<<((Kd)ai.first)<<"->"<<((Vd)ai.second)<<",";
798  }
799  if (a.size()>0) ss.seekp(-1, std::ios_base::end);
800  ss<<"}";
801  return ss.str();
802  }
804  template<typename K,typename V>
805  inline tString toString(const std::map<K,V>& a) {
806  std::stringstream ss;
807  ss<<"{";
808  for(const auto & ai : a) {
809  ss<<ai.first<<"->"<<ai.second<<",";
810  }
811  if (a.size()>0) ss.seekp(-1, std::ios_base::end);
812  ss<<"}";
813  return ss.str();
814  }
817  template<typename T> requires (functions_type::isArithmeticType<T>)
818  inline tBoolean parse(const tString& s,std::vector<T>& vs) {
819  //clear vs
820  vs.clear();
821  //array delim
822  std::array<tChar,3> beginArrayDelim={'{','[','('};
823  std::array<tChar,3> endArrayDelim={'}',']',')'};
824  //trim s
825  tString sexpr=s;
826  functions_string::trim(sexpr);
827 
828  //find the array delim
829  tIndex c=0;
830  for(const auto& b:beginArrayDelim) {
831  if (sexpr[0]==b) break;
832  c++;
833  }
834  if (c>=beginArrayDelim.size()) return false;
835  //get the end of the array delim
836  tIndex i=sexpr.find_last_of(endArrayDelim[c]);
837  if (i>=sexpr.length()) return false;
838 
839  //get the string between the array delim
840  T v;
841  tString token;
842  sexpr=sexpr.substr(1,i-1);
843  while (sexpr.length()>0) {//while tokens are in sexpr
844  i=sexpr.find_first_of(',');
845  if (i<sexpr.length()) {
846  //get the current token
847  token=sexpr.substr(0,i);
848  //set sexpr to next tokens
849  sexpr=sexpr.substr(i+1);
850  } else {
851  //get th current token
852  token=sexpr;
853  //set sexpr to null
854  sexpr="";
855  }
856  //convert the token to T type
857  functions_numeric::parse(token,v);
858  //register the token
859  vs.push_back(v);
860  }
861  return true;
862 
863  }
864 }
865 
866 #endif
867 
868 
Generated by doxygen 1.9.1