ConicBundle
matop.hxx
Go to the documentation of this file.
1 
2 
3 #ifndef CH_MATRIX_CLASSES__MATOP_HXX
4 #define CH_MATRIX_CLASSES__MATOP_HXX
5 
16 #include <limits.h>
17 #include <float.h>
18 #include <assert.h>
19 #if defined(WITH_BLAS) || defined(BLAS_GENMULT)
20 extern "C"{
21 #include "cblas.h"
22 }
23 #endif
24 #include <random>
25 #include "gb_rand.hxx"
26 #include "CBconfig.hxx"
27 
30 namespace CH_Matrix_Classes {
31 
32 
33  //---------------------------------------------------------------------------
34 
38 
40  typedef int Integer;
41 
43  const Integer max_Integer= INT_MAX;
44 
46  const Integer min_Integer= INT_MIN;
47 
48 
50  typedef double Real;
51 
53  const Real max_Real= DBL_MAX;
54 
56  const Real min_Real= -DBL_MAX;
57 
59  const Real eps_Real= DBL_EPSILON;
60 
62 
63  //---------------------------------------------------------------------------
64 
69 
70 
76 template<class Val>
77 inline void mat_xea( Integer len, Val *x, const Val a)
78 {
79  const Val* const xend=x+len;
80  for(;x!=xend;)
81  (*x++)=a;
82 }
83 
90 template<class Val>
91 inline void mat_xea( Integer len, Val *x, const Integer incx,
92  const Val a)
93 {
94  const Val* const xend=x+len*incx;
95  for(;x!=xend;x+=incx)
96  (*x)=a;
97 }
98 
99 #ifdef WITH_BLAS
100 
105  inline void mat_xey( Integer len, double *x, const double* y)
106  {
107  cblas_dcopy(len,y,1,x,1);
108  }
109 #endif
110 
116 template<class Val>
117 inline void mat_xey( Integer len, Val *x, const Val* y)
118 {
119  const Val* const xend=x+len;
120  for(;x!=xend;)
121  (*x++)=(*y++);
122 }
123 
124 #ifdef WITH_BLAS
125 
132  inline void mat_xey( Integer len, double *x, const Integer incx,
133  const double* y, const Integer incy)
134  {
135  cblas_dcopy(len,y,incy,x,incx);
136  }
137 #endif
138 
146 template<class Val>
147 inline void mat_xey( Integer len, Val *x, const Integer incx,
148  const Val *y, const Integer incy)
149 {
150  const Val* const xend=x+len*incx;
151  for(;x!=xend;x+=incx,y+=incy)
152  (*x)=(*y);
153 }
154 
155 #ifdef WITH_BLAS
156 
161 inline void mat_xmey( Integer len, double *x,
162  const double *y)
163  {
164  cblas_daxpy(len,-1.,y,1,x,1);
165  }
166 #endif
167 
173 template<class Val>
174 inline void mat_xmey( Integer len, Val *x, const Val *y)
175 {
176  const Val* const xend=x+len;
177  for(;x!=xend;)
178  (*x++)-=(*y++);
179 }
180 
181 #ifdef WITH_BLAS
182 
189 inline void mat_xmey( Integer len, double *x, const Integer incx,
190  const double *y, const Integer incy)
191  {
192  cblas_daxpy(len,-1.,y,incy,x,incx);
193  }
194 #endif
195 
203 template<class Val>
204 inline void mat_xmey( Integer len, Val *x, const Integer incx,
205  const Val *y, const Integer incy)
206 {
207  const Val* const xend=x+len*incx;
208  for(;x!=xend;x+=incx,y+=incy)
209  (*x)-=(*y);
210 }
211 
212 #ifdef WITH_BLAS
213 
217 inline void mat_xemx( Integer len, double *x)
218  {
219  cblas_dscal(len,-1.,x,1);
220  }
221 #endif
222 
227 template<class Val>
228 inline void mat_xemx( Integer len, Val *x)
229 {
230  const Val* const xend=x+len;
231  for(;x!=xend;x++)
232  (*x)=-(*x);
233 }
234 
235 #ifdef WITH_BLAS
236 
241 inline void mat_xemx( Integer len, double *x, const Integer incx)
242  {
243  cblas_dscal(len,-1.,x,incx);
244  }
245 #endif
246 
252 template<class Val>
253 inline void mat_xemx( Integer len, Val *x, const Integer incx)
254 {
255  const Val* const xend=x+len*incx;
256  for(;x!=xend;x+=incx)
257  (*x)=-(*x);
258 }
259 
260 #ifdef WITH_BLAS
261 
266 inline void mat_xemy( Integer len, double *x, const double *y)
267  {
268  cblas_dcopy(len,y,1,x,1);
269  cblas_dscal(len,-1.,x,1);
270  }
271 #endif
272 
278 template<class Val>
279 inline void mat_xemy( Integer len, Val *x, const Val *y)
280 {
281  const Val* const xend=x+len;
282  for(;x!=xend;)
283  (*x++)=-(*y++);
284 }
285 
286 #ifdef WITH_BLAS
287 
294 inline void mat_xemy( Integer len, double *x, const Integer incx,
295  const double *y, const Integer incy)
296  {
297  cblas_dcopy(len,y,incy,x,incx);
298  cblas_dscal(len,-1.,x,incx);
299  }
300 #endif
301 
309 template<class Val>
310 inline void mat_xemy( Integer len, Val *x, const Integer incx,
311  const Val *y, const Integer incy)
312 {
313  const Val* const xend=x+len*incx;
314  for(;x!=xend;x+=incx,y+=incy)
315  (*x)=-(*y);
316 }
317 
318 #ifdef WITH_BLAS
319 
325 inline void mat_xeya( Integer len, double *x,
326  const double *y, double a)
327  {
328  cblas_dcopy(len,y,1,x,1);
329  cblas_dscal(len,a,x,1);
330  }
331 #endif
332 
339 template<class Val>
340 inline void mat_xeya( Integer len, Val *x,
341  const Val *y, const Val a)
342 {
343  if (a==Val(0)){
344  mat_xea(len,x,Val(0));
345  return;
346  }
347  if (a==Val(1)){
348  mat_xey(len,x,y);
349  return;
350  }
351  if (a==Val(-1)){
352  mat_xemy(len,x,y);
353  return;
354  }
355  const Val* const xend=x+len;
356  for(;x!=xend;)
357  (*x++)=a*(*y++);
358 }
359 
360 #ifdef WITH_BLAS
361 
369 inline void mat_xeya( Integer len, double *x, const Integer incx,
370  const double *y, const Integer incy,double a)
371  {
372  cblas_dcopy(len,y,incy,x,incx);
373  cblas_dscal(len,a,x,incx);
374  }
375 #endif
376 
385 template<class Val>
386 inline void mat_xeya( Integer len, Val *x, const Integer incx,
387  const Val *y, const Integer incy, const Val a)
388 {
389  if (a==Val(0)){
390  mat_xea(len,x,incx,Val(0));
391  return;
392  }
393  if (a==Val(1)){
394  mat_xey(len,x,incx,y,incy);
395  return;
396  }
397  if (a==Val(-1)){
398  mat_xemy(len,x,incx,y,incy);
399  return;
400  }
401  const Val* const xend=x+len*incx;
402  for(;x!=xend;x+=incx,y+=incy)
403  (*x)=a*(*y);
404 }
405 
406 #ifdef WITH_BLAS
407 
412 inline void mat_xpey( Integer len, double *x,
413  const double *y)
414  {
415  cblas_daxpy(len,1.,y,1,x,1);
416  }
417 #endif
418 
424 template<class Val>
425 inline void mat_xpey( Integer len, Val *x, const Val *y)
426 {
427  const Val* const xend=x+len;
428  for(;x!=xend;)
429  (*x++)+=(*y++);
430 }
431 
432 #ifdef WITH_BLAS
433 
440 inline void mat_xpey( Integer len, double *x, const Integer incx,
441  const double *y, const Integer incy)
442  {
443  cblas_daxpy(len,1.,y,incy,x,incx);
444  }
445 #endif
446 
454 template<class Val>
455 inline void mat_xpey( Integer len, Val *x, const Integer incx,
456  const Val *y, const Integer incy)
457 {
458  const Val* const xend=x+len*incx;
459  for(;x!=xend;x+=incx,y+=incy)
460  (*x)+=(*y);
461 }
462 
463 
469 template<class Val>
470 inline void mat_xhadey( Integer len, Val *x, const Val *y)
471 {
472  const Val* const xend=x+len;
473  for(;x!=xend;)
474  (*x++)*=(*y++);
475 }
476 
482 template<class Val>
483 inline void mat_xinvhadey( Integer len, Val *x, const Val *y)
484 {
485  const Val* const xend=x+len;
486  for(;x!=xend;)
487  (*x++)/=(*y++);
488 }
489 
497 template<class Val>
498 inline void mat_xhadey( Integer len, Val *x, const Integer incx,
499  const Val *y, const Integer incy)
500 {
501  const Val* const xend=x+len*incx;
502  for(;x!=xend;x+=incx,y+=incy)
503  (*x)*=(*y);
504 }
505 
513 template<class Val>
514 inline void mat_xinvhadey( Integer len, Val *x, const Integer incx,
515  const Val *y, const Integer incy)
516 {
517  const Val* const xend=x+len*incx;
518  for(;x!=xend;x+=incx,y+=incy)
519  (*x)/=(*y);
520 }
521 
522 #ifdef WITH_BLAS
523 
529 inline void mat_xpeya( Integer len, Real *x,
530  const Real *y, const Real a)
531 {
532  cblas_daxpy(len,a,y,1,x,1);
533  }
534 #endif
535 
542 template<class Val>
543 inline void mat_xpeya( Integer len, Val *x,
544  const Val *y, const Val a)
545 {
546  if (a==Val(0))
547  return;
548  if (a==Val(1)){
549  mat_xpey(len,x,y);
550  return;
551  }
552  if (a==Val(-1)){
553  mat_xmey(len,x,y);
554  return;
555  }
556  const Val* const xend=x+len;
557  for(;x!=xend;)
558  (*x++)+=a*(*y++);
559 }
560 
561 #ifdef WITH_BLAS
562 
570 inline void mat_xpeya( Integer len, Real *x, const Integer incx,
571  const Real *y, const Integer incy, const Real a)
572 {
573  cblas_daxpy(len,a,y,incy,x,incx);
574  }
575 #endif
576 
585 template<class Val>
586 inline void mat_xpeya( Integer len, Val *x, const Integer incx,
587  const Val *y, const Integer incy, const Val a)
588 {
589  if (a==Val(0))
590  return;
591  if (a==Val(1)){
592  mat_xpey(len,x,incx,y,incy);
593  return;
594  }
595  if (a==Val(-1)){
596  mat_xmey(len,x,incx,y,incy);
597  return;
598  }
599  const Val* const xend=x+len*incx;
600  for(;x!=xend;x+=incx,y+=incy)
601  (*x)+=a*(*y);
602 }
603 
604 #ifdef WITH_BLAS
605 
612 inline void mat_xbpeya( Integer len, Real *x,
613  const Real *y, const Real a, const Real b)
614 {
615  cblas_dscal(len,b,x,1);
616  cblas_daxpy(len,a,y,1,x,1);
617 }
618 #endif
619 
627 template<class Val>
628 inline void mat_xbpeya( Integer len, Val *x,
629  const Val *y, const Val a,const Val b)
630 {
631  const Val* const xend=x+len;
632  if (b!=Val(1)){
633  for(;x!=xend;x++)
634  (*x)=b*(*x)+a*(*y++);
635  }
636  else {
637  for(;x!=xend;)
638  (*x++)+=a*(*y++);
639  }
640 }
641 
642 #ifdef WITH_BLAS
643 
652 inline void mat_xbpeya( Integer len, Real *x, const Integer incx,
653  const Real *y, const Integer incy, const Real a, const Real b)
654 {
655  cblas_dscal(len,b,x,incx);
656  cblas_daxpy(len,a,y,incy,x,incx);
657 }
658 #endif
659 
669 template<class Val>
670 inline void mat_xbpeya( Integer len, Val *x, const Integer incx,
671  const Val *y, const Integer incy, const Val a,const Val b)
672 {
673  const Val* const xend=x+len*incx;
674  if (b!=Val(1)){
675  for(;x!=xend;x+=incx,y+=incy)
676  (*x)=b*(*x)+a*(*y);
677  }
678  else {
679  for(;x!=xend;x+=incx,y+=incy)
680  (*x)+=a*(*y);
681  }
682 }
683 
689 template<class Val>
690 inline void mat_xpea( Integer len, Val *x, const Val a)
691 {
692  if (a==Val(0))
693  return;
694  const Val* const xend=x+len;
695  for(;x!=xend;)
696  (*x++)+=a;
697 }
698 
705 template<class Val>
706 inline void mat_xpea( Integer len, Val *x, const Integer incx,
707  const Val a)
708 {
709  if (a==Val(0))
710  return;
711  const Val* const xend=x+len*incx;
712  for(;x!=xend;x+=incx)
713  (*x)+=a;
714 }
715 
716 #ifdef WITH_BLAS
717 
722 inline void mat_xmultea( Integer len, double *x,double a)
723  {
724  cblas_dscal(len,a,x,1);
725  }
726 #endif
727 
733 template<class Val>
734 inline void mat_xmultea( Integer len, Val *x, const Val a)
735 {
736  if (a==Val(1))
737  return;
738  const Val* const xend=x+len;
739  for(;x!=xend;)
740  (*x++)*=a;
741 }
742 
743 #ifdef WITH_BLAS
744 
750 inline void mat_xmultea( Integer len, double *x, const Integer incx,double a)
751  {
752  cblas_dscal(len,a,x,incx);
753  }
754 #endif
755 
762 template<class Val>
763 inline void mat_xmultea( Integer len, Val *x, const Integer incx,
764  const Val a)
765 {
766  if (a==Val(1))
767  return;
768  const Val* const xend=x+len*incx;
769  for(;x!=xend;x+=incx)
770  (*x)*=a;
771 }
772 
773 #ifdef WITH_BLAS
774 
779 inline void mat_xdivea( Integer len, double *x,double a)
780  {
781  cblas_dscal(len,1./a,x,1);
782  }
783 #endif
784 
790 template<class Val>
791 inline void mat_xdivea( Integer len, Val *x, const Val a)
792 {
793  if (a==Val(1))
794  return;
795  const Val* const xend=x+len;
796  for(;x!=xend;)
797  (*x++)/=a;
798 }
799 
800 #ifdef WITH_BLAS
801 
807 inline void mat_xdivea( Integer len, double *x, const Integer incx,double a)
808  {
809  cblas_dscal(len,1./a,x,incx);
810  }
811 #endif
812 
819 template<class Val>
820 inline void mat_xdivea( Integer len, Val *x, const Integer incx,
821  const Val a)
822 {
823  if (a==Val(1))
824  return;
825  const Val* const xend=x+len*incx;
826  for(;x!=xend;x+=incx)
827  (*x)/=a;
828 }
829 
835 template<class Val>
836 inline void mat_xmodea( Integer len, Val *x, const Val a)
837 {
838  const Val* const xend=x+len;
839  for(;x!=xend;)
840  (*x++)%=a;
841 }
842 
849 template<class Val>
850 inline void mat_xmodea( Integer len, Val *x, const Integer incx,
851  const Val a)
852 {
853  const Val* const xend=x+len*incx;
854  for(;x!=xend;x+=incx)
855  (*x)%=a;
856 }
857 
858 #ifdef WITH_BLAS
859 
865 inline void mat_xeypz( Integer len, Real *x,
866  const Real *y,
867  const Real *z)
868 {
869  cblas_dcopy(len,y,1,x,1);
870  cblas_daxpy(len,1.,z,1,x,1);
871 }
872 #endif
873 
880 template<class Val>
881 inline void mat_xeypz( Integer len, Val *x,
882  const Val *y, const Val *z)
883 {
884  const Val* const xend=x+len;
885  for(;x!=xend;)
886  (*x++)=(*y++)+(*z++);
887 }
888 
889 #ifdef WITH_BLAS
890 
899 inline void mat_xeypz( Integer len, Real *x, const Integer incx,
900  const Real *y, const Integer incy,
901  const Real *z, const Integer incz)
902 {
903  cblas_dcopy(len,y,incy,x,incx);
904  cblas_daxpy(len,1.,z,incz,x,incx);
905 }
906 #endif
907 
917 template<class Val>
918 inline void mat_xeypz( Integer len, Val *x, const Integer incx,
919  const Val *y, const Integer incy,
920  const Val *z, const Integer incz)
921 {
922  const Val* const xend=x+len*incx;
923  for(;x!=xend;x+=incx,y+=incy,z+=incz)
924  (*x)=(*y)+(*z);
925 }
926 
927 #ifdef WITH_BLAS
928 
934 inline void mat_xeymz( Integer len, Real *x,
935  const Real *y,
936  const Real *z)
937 {
938  cblas_dcopy(len,y,1,x,1);
939  cblas_daxpy(len,-1.,z,1,x,1);
940 }
941 #endif
942 
949 template<class Val>
950 inline void mat_xeymz( Integer len, Val *x,
951  const Val *y, const Val *z)
952 {
953  const Val* const xend=x+len;
954  for(;x!=xend;)
955  (*x++)=(*y++)-(*z++);
956 }
957 
958 #ifdef WITH_BLAS
959 
968 inline void mat_xeymz( Integer len, Real *x, const Integer incx,
969  const Real *y, const Integer incy,
970  const Real *z, const Integer incz)
971 {
972  cblas_dcopy(len,y,incy,x,incx);
973  cblas_daxpy(len,-1.,z,incz,x,incx);
974 }
975 #endif
976 
986 template<class Val>
987 inline void mat_xeymz( Integer len, Val *x, const Integer incx,
988  const Val *y, const Integer incy,
989  const Val *z, const Integer incz)
990 {
991  const Val* const xend=x+len*incx;
992  for(;x!=xend;x+=incx,y+=incy,z+=incz)
993  (*x)=(*y)-(*z);
994 }
995 
996 #ifdef WITH_BLAS
997 
1005 inline void mat_xeyapzb( Integer len, Real *x,
1006  const Real *y,
1007  const Real *z, const Real a, const Real b)
1008 {
1009  cblas_dcopy(len,y,1,x,1);
1010  cblas_dscal(len,a,x,1);
1011  cblas_daxpy(len,b,z,1,x,1);
1012 }
1013 #endif
1014 
1023 template<class Val>
1024 inline void mat_xeyapzb( Integer len, Val *x,
1025  const Val *y, const Val *z,const Val a,const Val b)
1026 {
1027  const Val* const xend=x+len;
1028  for(;x!=xend;)
1029  (*x++)=a*(*y++)+b*(*z++);
1030 }
1031 
1032 #ifdef WITH_BLAS
1033 
1044 inline void mat_xeyapzb( Integer len, Real *x, const Integer incx,
1045  const Real *y, const Integer incy,
1046  const Real *z, const Integer incz, const Real a, const Real b)
1047 {
1048  cblas_dcopy(len,y,incy,x,incx);
1049  cblas_dscal(len,a,x,incx);
1050  cblas_daxpy(len,b,z,incz,x,incx);
1051 }
1052 #endif
1053 
1065 template<class Val>
1066 inline void mat_xeyapzb( Integer len, Val *x, const Integer incx,
1067  const Val *y, const Integer incy,
1068  const Val *z, const Integer incz,const Val a,const Val b)
1069 {
1070  const Val* const xend=x+len*incx;
1071  for(;x!=xend;x+=incx,y+=incy,z+=incz)
1072  (*x)=a*(*y)+b*(*z);
1073 }
1074 
1075 #ifdef WITH_BLAS
1076 
1082 inline double mat_ip( Integer len, const double *x,const double *y)
1083 {
1084  return cblas_ddot(len,x,1,y,1);
1085 }
1086 #endif
1087 
1095 template<class Val>
1096 inline Val mat_ip( Integer len, const Val *x, const Val *y,const Val* d=0)
1097 {
1098  Val sum=0;
1099  const Val* const xend=x+len;
1100  if (d==0)
1101  for(;x!=xend;)
1102  sum+=(*x++)*(*y++);
1103  else {
1104  for(;x!=xend;)
1105  sum+=(*x++)*(*y++)*(*d++);
1106  }
1107  return sum;
1108 }
1109 
1110 
1111 #ifdef WITH_BLAS
1112 
1120 inline double mat_ip( Integer len, const double *x, const Integer incx,
1121  const double *y, const Integer incy)
1122 {
1123  return cblas_ddot(len,x,incx,y,incy);
1124 }
1125 #endif
1126 
1137 template<class Val>
1138 inline Val mat_ip( Integer len, const Val *x, const Integer incx,
1139  const Val *y, const Integer incy,const Val* d=0,const Integer incd=1)
1140 {
1141  Val sum=0;
1142  const Val* const xend=x+len*incx;
1143  if (d==0)
1144  for(;x!=xend;x+=incx,y+=incy)
1145  sum+=(*x)*(*y);
1146  else
1147  for(;x!=xend;x+=incx,y+=incy,d+=incd)
1148  sum+=(*x)*(*y)*(*d);
1149  return sum;
1150 }
1151 
1165 template<class Val>
1166 inline Val mat_ip_dense_sparse(const Integer lenx, const Val *x, Integer leny, const Val *yval, const Integer* yind,const Val* d=0)
1167 {
1168  Val sum=0;
1169  const Integer* const yend=yind+leny;
1170  if (d==0)
1171  while((yind!=yend)&&(*yind<lenx))
1172  sum+=(*(x+*yind++))*(*yval++);
1173  else {
1174  while((yind!=yend)&&(*yind<lenx)){
1175  sum+=(*(x+*yind))*(*yval++)*(*(d+*yind));
1176  yind++;
1177  }
1178  }
1179  return sum;
1180 }
1181 
1197 template<class Val>
1198 inline Val mat_ip_sparse_sparse(const Integer lenx, const Val *xval, const Integer* xind,const Integer leny, const Val *yval, const Integer* yind,const Val* d)
1199 {
1200  if ((lenx<=0)||(leny<=0))
1201  return 0.;
1202  Val sum=0;
1203  const Integer* const xend=xind+lenx;
1204  const Integer* const yend=yind+leny;
1205  if (d==0) {
1206  do {
1207  if (*xind==*yind){
1208  sum+=(*xval++)*(*yval++);
1209  xind++;
1210  if (xind==xend)
1211  break;
1212  yind++;
1213  }
1214  while ((yind!=yend)&&(*yind<*xind)){
1215  yval++;
1216  yind++;
1217  }
1218  if (yind==yend)
1219  break;
1220  while ((xind!=xend)&&(*xind<*yind)){
1221  xval++;
1222  xind++;
1223  }
1224  } while(xind!=xend);
1225  }
1226  else {
1227  do {
1228  if (*xind==*yind){
1229  sum+=(*xval++)*(*yval++)*(*(d+*xind));
1230  xind++;
1231  if (xind==xend)
1232  break;
1233  yind++;
1234  }
1235  while ((yind!=yend)&&(*yind<*xind)){
1236  yval++;
1237  yind++;
1238  }
1239  if (yind==yend)
1240  break;
1241  while ((xind!=xend)&&(*xind<*yind)){
1242  xval++;
1243  xind++;
1244  }
1245  } while(xind!=xend);
1246  }
1247  return sum;
1248 }
1249 /*
1250 {
1251  Val sum=0.;
1252  if((lenx>0)&&(leny>0)){
1253  const Integer* const xend=xind+lenx;
1254  const Integer* const yend=yind+leny;
1255  while(true){
1256  Integer diff=*xind-*yind;
1257  if (diff==0){
1258  sum+=(*xval++)*(*yval++);
1259  if(++xind==xend)
1260  break;
1261  if(++yind==yend)
1262  break;
1263  }
1264  else {
1265  if (diff>0){
1266  yval++;
1267  if (++yind==yend)
1268  break;
1269  }
1270  else {
1271  xval++;
1272  if (++xind==xend)
1273  break;
1274  }
1275  }
1276  }
1277  }
1278  return sum;
1279 }
1280 */
1281  /*
1282 {
1283  Val sum=0.;
1284  if((lenx>0)&&(leny>0)){
1285  while(true){
1286  Integer diff=*xind-*yind;
1287  if (diff==0){
1288  sum+=(*xval++)*(*yval++);
1289  xind++;yind++;
1290  if ((--lenx==0)||(--leny==0))
1291  break;
1292  }
1293  else {
1294  if (diff>0){
1295  yval++;
1296  yind++;
1297  if (--leny==0)
1298  break;
1299  }
1300  else {
1301  xval++;
1302  xind++;
1303  if(--lenx==0)
1304  break;
1305  }
1306  }
1307  }
1308  }
1309  return sum;
1310 }
1311  */
1312 
1313 
1314 #ifdef WITH_BLAS
1315 
1320 inline double mat_ip( Integer len, const double *x)
1321 {
1322  return cblas_ddot(len,x,1,x,1);
1323 }
1324 #endif
1325 
1331 template<class Val>
1332 inline Val mat_ip( Integer len, const Val *x)
1333 {
1334  Val sum=0;
1335  const Val* const xend=x+len;
1336  while(x!=xend) {
1337  const Val d=(*x++);
1338  sum+=d*d;
1339  }
1340  return sum;
1341 }
1342 
1343 
1344 #ifdef WITH_BLAS
1345 
1351 inline double mat_ip( Integer len, const double *x, const Integer incx)
1352 {
1353  return cblas_ddot(len,x,incx,x,incx);
1354 }
1355 #endif
1356 
1363 template<class Val>
1364 inline Val mat_ip( Integer len, const Val *x, const Integer incx)
1365 {
1366  Val sum=0;
1367  const Val* const xend=x+len*incx;
1368  for(;x!=xend;x+=incx) {
1369  const Val d=(*x);
1370  sum+=d*d;
1371  }
1372  return sum;
1373 }
1374 
1375 
1376 
1382 template<class Val>
1383 inline Val mat_sum( Integer len, const Val *x)
1384 {
1385  Val sum=0;
1386  const Val* const xend=x+len;
1387  for(;x!=xend;)
1388  sum+=(*x++);
1389  return sum;
1390 }
1391 
1398 template<class Val>
1399 inline Val mat_sum( Integer len, const Val *x,const Integer incx)
1400 {
1401  Val sum=0;
1402  const Val* const xend=x+len*incx;
1403  for(;x!=xend;x+=incx)
1404  sum+=(*x);
1405  return sum;
1406 }
1407 
1414 template<class Val>
1415 inline bool mat_equal( Integer len, const Val *x,const Val *y)
1416 {
1417  const Val* const xend=x+len;
1418  while(x!=xend){
1419  if (*x!=*y)
1420  break;
1421  x++,y++;
1422  }
1423  return (x==xend);
1424 }
1425 
1426 
1427 #ifdef WITH_BLAS
1428 
1433 inline void mat_swap( Integer len, double *x, double *y)
1434 {
1435  cblas_dswap(len,x,1,y,1);
1436 }
1437 #endif
1438 
1444 template<class Val>
1445 inline void mat_swap( Integer len, Val *x, Val *y)
1446 {
1447  const Val* const xend=x+len;
1448  for(;x!=xend;)
1449  {Val h=*x;(*x++)=*y;(*y++)=h;}
1450 }
1451 
1452 #ifdef WITH_BLAS
1453 
1460 inline void mat_swap( Integer len, double *x, const Integer incx, double *y, const Integer incy)
1461 {
1462  cblas_dswap(len,x,incx,y,incy);
1463 }
1464 #endif
1465 
1473 template<class Val>
1474 inline void mat_swap( Integer len, Val *x, const Integer incx,
1475  Val *y, const Integer incy)
1476 {
1477  const Val* const xend=x+len*incx;
1478  for(;x!=xend;x+=incx,y+=incy)
1479  {Val h=*x;(*x)=*y;(*y)=h;}
1480 }
1481 
1483 
1484 
1487 template<class Val>
1490  mat_less_index( const Val* values ):
1491  _values(values)
1492  {}
1493 
1495  bool operator()( Integer i, Integer j) const {
1496  return _values[i] < _values[j];
1497  }
1498 
1499 private:
1501  const Val* _values;
1502 
1503 };
1504 
1505 
1508 template<class Val>
1511  mat_greater_index( const Val* values ):
1512  _values(values)
1513  {}
1514 
1516  bool operator()( Integer i, Integer j) const {
1517  return _values[i] > _values[j];
1518  }
1519 
1520 private:
1522  const Val* _values;
1523 
1524 };
1525 
1526 
1527 
1528 //-----------------------------------------------------------------------------
1529 
1534 
1535 
1536 // **************************************************************************
1537 // mat_randgen
1538 // **************************************************************************
1539 
1540 
1546  extern std::default_random_engine mat_std_randgen;
1547 extern std::mt19937 mat_mt_randgen;
1548 extern std::mt19937_64 mat_mt64_randgen;
1549 
1550 
1551 // **************************************************************************
1552 // materrout
1553 // **************************************************************************
1554 
1555 
1558 extern std::ostream* materrout;
1559 
1561 
1562 
1563 
1564 
1565 
1566 //-----------------------------------------------------------------------------
1567 
1579 
1580 // **************************************************************************
1581 // Matrix Types
1582 // **************************************************************************
1583 
1585 enum Mtype {
1592 };
1593 
1594 // **************************************************************************
1595 // "exceptions"
1596 // **************************************************************************
1597 
1599 enum MEcode {
1606 };
1607 
1610 {
1611 public:
1613  const char *message;
1615 
1617  MatrixError(MEcode c=ME_unspec,const char *mes=0,Mtype mt=MTmatrix):
1618  code(c),message(mes),mtype(mt){}
1619 };
1620 
1622 class MErange:public MatrixError
1623 {
1624 public:
1625  Integer nr;
1626  Integer r;
1627  Integer nc;
1628  Integer c;
1629 
1631  MErange(Integer inr,Integer ir,Integer inc,Integer ic,const char *mes=0,Mtype mt=MTmatrix):
1632  MatrixError(ME_range,mes,mt),nr(inr),r(ir),nc(inc),c(ic){}
1633 };
1634 
1636 class MEmem:public MatrixError
1637 {
1638 public:
1639  Integer size;
1640 
1642  MEmem(Integer s,const char *mes=0,Mtype mt=MTmatrix):
1643  MatrixError(ME_mem,mes,mt),size(s){}
1644 };
1645 
1647 class MEdim:public MatrixError
1648 {
1649 public:
1650  Integer nr1;
1651  Integer nc1;
1652  Integer nr2;
1653  Integer nc2;
1654 
1656  MEdim(Integer r1,Integer c1,Integer r2,Integer c2,const char *mes=0,Mtype mt=MTmatrix):
1657  MatrixError(ME_dim,mes,mt),nr1(r1),nc1(c1),nr2(r2),nc2(c2){}
1658 };
1659 
1661 int MEmessage(const MatrixError&);
1662 
1663 
1664 // **************************************************************************
1665 // CONICBUNDLE_DEBUG error checking templates
1666 // **************************************************************************
1667 
1668 #if (CONICBUNDLE_DEBUG>=1)
1669 
1671 template<class Mat>
1672 inline void chk_init(const Mat& A)
1673 {
1674  if (!A.get_init())
1675  MEmessage(MatrixError(ME_unspec,"Matrix not initialized",A.get_mtype()));
1676 }
1677 
1679 template<class Mat1,class Mat2>
1680 inline void chk_add(const Mat1& A,const Mat2& B)
1681 {
1682  chk_init(A); chk_init(B);
1683  if ((A.dim()==Integer(0))&&(B.dim()==Integer(0))) return;
1684  if ((A.rowdim()!=B.rowdim())||(A.coldim()!=B.coldim()))
1685  MEmessage(MEdim(A.rowdim(),A.coldim(),B.rowdim(),B.coldim(),
1686  "dimensions do not match in additive expression",
1687  A.get_mtype()));
1688 }
1689 
1691 template<class Mat1,class Mat2>
1692 inline void chk_mult(const Mat1& A,const Mat2& B,int atrans=0,int btrans=0)
1693 {
1694  chk_init(A); chk_init(B);
1695  if ((A.dim()==Integer(0))&&(B.dim()==Integer(0))) return;
1696  if (
1697  ((atrans==0)&&(btrans==0)&&(A.coldim()!=B.rowdim()))
1698  ||
1699  ((atrans==1)&&(btrans==0)&&(A.rowdim()!=B.rowdim()))
1700  ||
1701  ((atrans==0)&&(btrans==1)&&(A.coldim()!=B.coldim()))
1702  ||
1703  ((atrans==1)&&(btrans==1)&&(A.rowdim()!=B.coldim()))
1704  )
1705  MEmessage(MEdim(A.rowdim(),A.coldim(),B.rowdim(),B.coldim(),
1706  "dimensions do not match in multiplicative expression",
1707  A.get_mtype()));
1708 }
1709 
1711 template<class Mat1,class Mat2>
1712 inline void chk_rowseq(const Mat1& A,const Mat2& B)
1713 {
1714  chk_init(A); chk_init(B);
1715  if ((A.dim()==Integer(0))&&(B.dim()==Integer(0))) return;
1716 
1717  if (A.rowdim()!=B.rowdim())
1718  MEmessage(MEdim(A.rowdim(),A.coldim(),B.rowdim(),B.coldim(),
1719  "number of rows do not agree",
1720  A.get_mtype()));
1721 }
1722 
1724 template<class Mat1,class Mat2>
1725 inline void chk_colseq(const Mat1& A,const Mat2& B)
1726 {
1727  chk_init(A); chk_init(B);
1728  if ((A.dim()==Integer(0))&&(B.dim()==Integer(0))) return;
1729  if (A.coldim()!=B.coldim())
1730  MEmessage(MEdim(A.rowdim(),A.coldim(),B.rowdim(),B.coldim(),
1731  "number of columns do not agree",
1732  A.get_mtype()));
1733 }
1734 
1736 inline void chk_range(Integer r,Integer c,Integer ubr,Integer ubc)
1737 {
1738  if ((r<0)||(c<0)||
1739  ((ubr>=0)&&(r>=ubr))||
1740  ((ubc>=0)&&(c>=ubc)))
1741  MEmessage(MErange(r,ubr,c,ubc,"index out of range or negative dimension"));
1742 }
1743 
1745 inline void chk_single_range(Integer i,Integer ub)
1746 {
1747  if ((i<0)||
1748  ((ub>=0)&&(i>=ub)))
1749  MEmessage(MErange(i,ub,0,0,"index out of range or negative dimension"));
1750 }
1751 
1753 template<class Mat>
1754 inline void chk_set_init(Mat &A,bool val){A.set_init(val);}
1755 
1756 #else
1757 
1758 // if CONICBUNDLE_DEBUG is not definied all chk_functions are mapped to the null string
1759 
1761 #define chk_init(x)
1762 #define chk_add(x,y)
1764 template<class Mat1,class Mat2>
1766 inline void chk_mult(const Mat1&,const Mat2&,int =0,int =0){}
1768 #define chk_rowseq(x,y)
1769 #define chk_colseq(x,y)
1771 #define chk_range(i,j,ubi,ubj)
1773 #define chk_set_init(x,y)
1775 #define chk_single_range(i,ub)
1777 
1778 #endif
1779 
1781 
1782 }
1783 
1784 #endif
#define chk_single_range(i, ub)
CONICBUNDLE_DEBUG being undefined, the template function is removed. Otherwise it would whether index...
Definition: matop.hxx:1776
#define chk_set_init(x, y)
CONICBUNDLE_DEBUG being undefined, the template function is removed. Otherwise it would allow to set ...
Definition: matop.hxx:1774
int Integer
all integer numbers in calculations and indexing are of this type
Definition: matop.hxx:40
"less"-routine for sorting indices of value arrays by std::sort
Definition: matop.hxx:1488
void mat_xmey(Integer len, Val *x, const Val *y)
Set x[i]=-y[i] for len elements of the arrays x and y.
Definition: matop.hxx:174
void mat_xemy(Integer len, Val *x, const Val *y)
Set x[i]=-y[i] for len elements of the arrays x and y.
Definition: matop.hxx:279
void mat_xmultea(Integer len, Val *x, const Val a)
Set x[i]*=a for len elements of the array x.
Definition: matop.hxx:734
MEmem(Integer s, const char *mes=0, Mtype mt=MTmatrix)
constructor
Definition: matop.hxx:1642
error arises in a message of CH_Matrix_Classes::Indexmatrix
Definition: matop.hxx:1587
std::ostream * materrout
if not zero, this is the output stream for runtime error messages, by default it is set to &std::cout...
double Real
all real numbers in calculations are of this type
Definition: matop.hxx:50
error arises in a global function
Definition: matop.hxx:1586
void mat_xemx(Integer len, Val *x)
Set x[i]=-x[i] for an array of length len.
Definition: matop.hxx:228
mat_greater_index(const Val *values)
constructor initializing _values
Definition: matop.hxx:1511
void chk_mult(const Mat1 &, const Mat2 &, int=0, int=0)
CONICBUNDLE_DEBUG being undefined, the template function is removed. Otherwise it would check...
Definition: matop.hxx:1766
MEcode code
see MEcode for allowed error types
Definition: matop.hxx:1612
std::default_random_engine mat_std_randgen
optional fast alternative random number generator from std
void mat_xea(Integer len, Val *x, const Val a)
Set x[i]=a for len elements of the array x.
Definition: matop.hxx:77
Integer nr
number of rows
Definition: matop.hxx:1625
Mtype
serves for specifying the source (matrix class or function) of the error
Definition: matop.hxx:1585
Header declaring and (inline) implementing the random number generator CH_Tools::GB_rand.
Integer nr1
number of rows of object 1
Definition: matop.hxx:1650
Integer r
row index
Definition: matop.hxx:1626
Integer nc
number of columsn
Definition: matop.hxx:1627
void mat_xeya(Integer len, Val *x, const Val *y, const Val a)
Set x[i]=a*y[i] for len elements of the arrays x and y.
Definition: matop.hxx:340
const Real eps_Real
machine epsilon for type Real
Definition: matop.hxx:59
mat_less_index(const Val *values)
constructor initializing _values
Definition: matop.hxx:1490
Integer nr2
number of rows of object 2
Definition: matop.hxx:1652
void mat_xeypz(Integer len, Val *x, const Val *y, const Val *z)
Set x[i]=y[i]+z[i] for len elements of the arrays x, y and z.
Definition: matop.hxx:881
void mat_xeymz(Integer len, Val *x, const Val *y, const Val *z)
Set x[i]=y[i]-z[i] for len elements of the arrays x, y and z.
Definition: matop.hxx:950
void mat_xpey(Integer len, Val *x, const Val *y)
Set x[i]+=y[i] for len elements of the arrays x and y.
Definition: matop.hxx:425
void mat_xhadey(Integer len, Val *x, const Val *y)
Set x[i]*=y[i] for len elements of the arrays x and y.
Definition: matop.hxx:470
error arises in a message of CH_Matrix_Classes::Matrix
Definition: matop.hxx:1588
const Val * _values
const pointer to values
Definition: matop.hxx:1522
error due to numerical reasons (e.g. division by zero)
Definition: matop.hxx:1604
const Val * _values
const pointer to values
Definition: matop.hxx:1501
Integer nc1
number of columns of object 1
Definition: matop.hxx:1651
MErange(Integer inr, Integer ir, Integer inc, Integer ic, const char *mes=0, Mtype mt=MTmatrix)
constructor
Definition: matop.hxx:1631
const Integer min_Integer
minimal attainable value by an Integer
Definition: matop.hxx:46
error due to insufficient memory
Definition: matop.hxx:1602
Such an object is generated and passed to MEmessage(), whenever an error occurs. It holds some output...
Definition: matop.hxx:1609
int MEmessage(const MatrixError &)
displays an error message and terminates via abort() or returns 1 in case of warnings.
const char * message
the error mussage pointer must point to existing object and does not get freed
Definition: matop.hxx:1613
device independent random number generator based on long int with seed
Definition: gb_rand.hxx:28
void mat_swap(Integer len, Val *x, Val *y)
swap values x[i] and y[i] for len elements of the arrays x and y.
Definition: matop.hxx:1445
error arises in a message of CH_Matrix_Classes::Sparsemat
Definition: matop.hxx:1590
error due to range check
Definition: matop.hxx:1601
Integer nc2
number of columns of object 2
Definition: matop.hxx:1653
Mtype mtype
see Mtype for allowed error source types
Definition: matop.hxx:1614
error arises in a message of CH_Matrix_Classes::Sparsesym
Definition: matop.hxx:1591
bool mat_equal(Integer len, const Val *x, const Val *y)
returns true if the elements of the arrays x and y are exactly equal.
Definition: matop.hxx:1415
MatrixError(MEcode c=ME_unspec, const char *mes=0, Mtype mt=MTmatrix)
there is only this one constructor and no other messages
Definition: matop.hxx:1617
error due to inconsistent dimesions
Definition: matop.hxx:1603
MEcode
serves for specifying the error type.
Definition: matop.hxx:1599
Such an object is generated and passed to MEmessage() whenever a memory allocation fails...
Definition: matop.hxx:1636
Matrix Classes and Linear Algebra. See Matrix Classes (namespace CH_Matrix_Classes) for a quick intro...
Definition: PSCOracle.hxx:20
const Real min_Real
minimal attainable value by a Real
Definition: matop.hxx:56
std::mt19937 mat_mt_randgen
optional high quality alternative random number generator from std
#define chk_add(x, y)
CONICBUNDLE_DEBUG being undefined, the template function is removed. Otherwise it would check...
Definition: matop.hxx:1763
void mat_xey(Integer len, Val *x, const Val *y)
Copy an array of length len to destination x from source y.
Definition: matop.hxx:117
Val mat_ip_sparse_sparse(const Integer lenx, const Val *xval, const Integer *xind, const Integer leny, const Val *yval, const Integer *yind, const Val *d)
return sum(xval[i]*yval[j] for i,j with xind[i]==yind[j]) summing over elements of the sparse array r...
Definition: matop.hxx:1198
void mat_xinvhadey(Integer len, Val *x, const Val *y)
Set x[i]/=y[i] for len elements of the arrays x and y, no zero checking!
Definition: matop.hxx:483
void mat_xmodea(Integer len, Val *x, const Val a)
Set x[i]%=a for len elements of the array x.
Definition: matop.hxx:836
std::mt19937_64 mat_mt64_randgen
optional high quality alternative random number generator from std
const Real max_Real
maximal attainable value by a Real
Definition: matop.hxx:53
Val mat_ip(Integer len, const Val *x, const Val *y, const Val *d=0)
return sum(x[i]*y[i]) summing over len elements of the arrays x and y.
Definition: matop.hxx:1096
void mat_xpeya(Integer len, Val *x, const Val *y, const Val a)
Set x[i]+=a*y[i] for len elements of the arrays x and y.
Definition: matop.hxx:543
Such an object is generated and passed to MEmessage() whenever matrix dimensions do not agree for a d...
Definition: matop.hxx:1647
unspecified error type, not in the list below
Definition: matop.hxx:1600
Val mat_ip_dense_sparse(const Integer lenx, const Val *x, Integer leny, const Val *yval, const Integer *yind, const Val *d=0)
return sum(x[yind[j]]*yval[j]) summing over elements of the dense array x and a sparse array represen...
Definition: matop.hxx:1166
void mat_xbpeya(Integer len, Val *x, const Val *y, const Val a, const Val b)
Set x[i]=a*y[i]+b*x[i] for len elements of the arrays x and y.
Definition: matop.hxx:628
CH_Tools::GB_rand mat_randgen
common random number generator for use when a random matrix is generated.
#define chk_range(i, j, ubi, ubj)
CONICBUNDLE_DEBUG being undefined, the template function is removed. Otherwise it would check...
Definition: matop.hxx:1772
Integer c
column index
Definition: matop.hxx:1628
#define chk_colseq(x, y)
CONICBUNDLE_DEBUG being undefined, the template function is removed. Otherwise it would check...
Definition: matop.hxx:1770
error arises in a message of CH_Matrix_Classes::Symmatrix
Definition: matop.hxx:1589
bool operator()(Integer i, Integer j) const
returns (values[i] > _values[j])
Definition: matop.hxx:1516
Integer sum(const Indexmatrix &A)
returns the sum over all elements of A, i.e., (1 1 ... 1)*A*(1 1 ... 1)^T
MEdim(Integer r1, Integer c1, Integer r2, Integer c2, const char *mes=0, Mtype mt=MTmatrix)
constructor
Definition: matop.hxx:1656
#define chk_init(x)
CONICBUNDLE_DEBUG being undefined, the template function is removed. Otherwise it would check...
Definition: matop.hxx:1761
Such an object is generated and passed to MEmessage() whenever some index is out of range...
Definition: matop.hxx:1622
bool operator()(Integer i, Integer j) const
returns (values[i] < _values[j])
Definition: matop.hxx:1495
void mat_xeyapzb(Integer len, Val *x, const Val *y, const Val *z, const Val a, const Val b)
Set x[i]=a*y[i]+b*z[i] for len elements of the arrays x, y and z.
Definition: matop.hxx:1024
void mat_xpea(Integer len, Val *x, const Val a)
Set x[i]+=a for len elements of the array x.
Definition: matop.hxx:690
#define chk_rowseq(x, y)
CONICBUNDLE_DEBUG being undefined, the template function is removed. Otherwise it would check...
Definition: matop.hxx:1768
no error, but a possible source of difficulties
Definition: matop.hxx:1605
Integer size
memory size requested
Definition: matop.hxx:1639
void mat_xdivea(Integer len, Val *x, const Val a)
Set x[i]/=a for len elements of the array x.
Definition: matop.hxx:791
const Integer max_Integer
maximal attainable value by an Integer
Definition: matop.hxx:43
Val mat_sum(Integer len, const Val *x)
returns sum(x[i]) over len elements of the array x.
Definition: matop.hxx:1383
"greater"-routine for sorting indices of value arrays by std::sort
Definition: matop.hxx:1509