CVM Class Library  8.1
This C++ class library encapsulates concepts of vector and different matrices including square, band, symmetric and hermitian ones in Euclidean space of real and complex numbers.
 All Classes Files Functions Variables Typedefs Friends Macros Pages
cvm.h
Go to the documentation of this file.
1 // CVM Class Library
2 // http://cvmlib.com
3 //
4 // Copyright Sergei Nikolaev 1992-2014
5 // Distributed under the Boost Software License, Version 1.0.
6 // (See accompanying file LICENSE_1_0.txt or copy at
7 // http://www.boost.org/LICENSE_1_0.txt)
32 #ifndef _CVM_H
33 #define _CVM_H
34 
35 // 5.7 ILP64 support
36 #if defined (CVM_ILP64)
37  typedef long long int tint;
38 # define TINT_ZERO (0LL)
39 # define TINT_ONE (1LL)
40 # define CVM_TINT_FORMAT "%lld"
41 #else
42  typedef int tint;
43 # define TINT_ZERO (0)
44 # define TINT_ONE (1)
45 # define CVM_TINT_FORMAT "%d"
46 #endif
47 
48 // 5.7 0-based indexing
49 #if defined (CVM_ZERO_BASED)
50 # define CVM0 TINT_ZERO
51 #else
52 # define CVM0 TINT_ONE
53 #endif
54 
55 #if !defined (CVM_NO_MT) // as of 5.7 it's by default
56 # define CVM_MT
57 # if !defined (_PTHREADS)
58 # define _PTHREADS
59 # endif
60 #endif
61 
62 // MSVC++ 6.0 and higher settings
63 #if defined (_MSC_VER)
64 # pragma once
65 # define WIN32_LEAN_AND_MEAN // Exclude rarely-used stuff from Windows headers
66 # ifndef _WIN32_WINNT
67 # define _WIN32_WINNT 0x500 // at least Win2000 is required for InitializeCriticalSectionAndSpinCount
68 # endif
69 # include <windows.h>
70 # include <process.h>
71 // 6.0 TR1 dependency added
72 # include <memory>
73 
74 # if (_MSC_VER < 1400)
75 # error "Please use stable version 5.2 for older MSVC compilers"
76 # endif
77 # if (_MSC_VER >= 1700) // since VC11 aka MS Visual Studio 2012
78 # define CVM_STD_MUTEX
79 # include <mutex>
80 # include <thread>
81 # endif
82 # if (_MSC_VER >= 1800) // since VC12 aka MS Visual Studio 2013
83 # define CVM_USE_VARIADIC_TEMPLATES
84 # define CVM_USE_INITIALIZER_LISTS
85 # define CVM_USE_DELEGATING_CONSTRUCTORS
86 
87 # include <initializer_list>
88 # endif
89 # if (!defined(__INTEL_COMPILER) || !defined(_WIN64)) && !defined(CVM_ACML) && !(_MSC_VER >= 1500 && defined(_WIN64))
90 # define CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES
91 # endif
92 # if defined(CVM_ACML)
93 # define CVM_COMPLEX_NUMBER_RETURNED
94 # endif
95 //# if (_MSC_VER <= 1800)
96 # define noexcept throw()
97 //# endif
98 # pragma warning(disable:4290)
99 # if defined (CVM_FLOAT)
100 # pragma warning(disable:4244)
101 # endif
102 # if defined (SRC_EXPORTS) && !defined (CVM_EXPORTS)
103 # define CVM_EXPORTS
104 # endif
105 # include <unordered_map>
106 # define CVM_BLOCKS_MAP std::unordered_map
107 # ifdef CVM_STATIC
108 # define CVM_API
109 # else
110 # ifdef CVM_EXPORTS
111 # define CVM_API __declspec(dllexport)
112 # else
113 # define CVM_API __declspec(dllimport)
114 # endif
115 # endif
116 
117 # include <limits>
118 typedef __int64 CVM_LONGEST_INT;
119 
120 # if defined(_WIN64)
121 typedef unsigned long long CVM_PTR_WRAPPER;
122 # else
123 # if defined (CVM_ILP64)
124 # error "CVM_ILP64 is incompatible with 32 bit mode"
125 # endif
126 typedef unsigned long CVM_PTR_WRAPPER;
127 # endif
128 
129 # define CVM_VSNPRINTF vsnprintf_s
130 # define CVM_VSNPRINTF_S_DEFINED
131 # define CVM_STRCPY_S_DEFINED
132 
133 // GCC settings
134 #elif defined (__GNUC__)
135 # if defined (__MINGW32__)
136 # define WIN32_LEAN_AND_MEAN
137 # include <windows.h>
138 # include <process.h>
139 
140 # define CVM_USE_CRITICAL_SECTION_NOT_MUTEX
141 
142 # else
143 # ifdef __stdcall
144 # undef __stdcall
145 # endif
146 # define __stdcall
147 # include <semaphore.h> // Unix
148 # endif
149 
150 // 8.0 - since 4.7.0 we use new std::mutex features
151 # if !defined(__INTEL_COMPILER) && !defined (__MINGW32__) && (__GNUC__ > 4 || __GNUC__ == 4 && __GNUC_MINOR__ >= 7)
152 # define CVM_USE_DELEGATING_CONSTRUCTORS
153 # define CVM_USE_USER_LITERALS
154 # define CVM_STD_MUTEX
155 # include <mutex>
156 # include <thread>
157 # else
158 # define override
159 # endif
160 
161 // 8.1 - more C++11 features, see also http://gcc.gnu.org/projects/cxx0x.html
162 # if !defined (__MINGW32__) && (__GNUC__ > 4 || __GNUC__ == 4 && __GNUC_MINOR__ >= 4)
163 # define CVM_USE_VARIADIC_TEMPLATES
164 # define CVM_USE_INITIALIZER_LISTS
165 # include <initializer_list>
166 # endif
167 
168 // 6.0 TR1 dependency added
169 # include <memory>
170 
171 # define CVM_API
172 # define CVM_STDEXT stdext
173 # define CVM_BLOCKS_MAP std::map
174 
175 typedef long long CVM_LONGEST_INT;
176 
177 # if defined(__AMD64__) || defined(_WIN64) || defined(__x86_64__)
178 typedef unsigned long long CVM_PTR_WRAPPER;
179 # else
180 # if defined (CVM_ILP64)
181 # error "CVM_ILP64 is incompatible with 32 bit mode"
182 # endif
183 typedef unsigned long CVM_PTR_WRAPPER;
184 # endif
185 
186 # define CVM_VSNPRINTF vsnprintf
187 
188 // 5.7 - looking for memset
189 # include <string.h>
190 // gcc 4.6.1 fix for ptrdiff_t
191 # include <cstddef>
192 
193 //clang -dM -E -x c /dev/null
194 # if defined(__clang__) && __clang_major__ <= 5
195 # define CVM_NO_DEFAULT_RANDOM_ENGINE_SUPPORTED 1
196 # endif
197 # if defined(__clang__) && __clang_major__ >= 5
198 # define CVM_USE_VARIADIC_TEMPLATES
199 # define CVM_USE_INITIALIZER_LISTS
200 # define CVM_USE_USER_LITERALS
201 # define CVM_USE_DELEGATING_CONSTRUCTORS
202 # endif
203 
204 #else
205 # error "Unsupported compiler"
206 #endif
207 
208 #include <stdlib.h>
209 #include <stdio.h>
210 #include <stdarg.h>
211 #include <math.h>
212 #include <float.h>
213 #include <time.h>
214 #include <iostream>
215 #include <list>
216 #include <map>
217 #include <string>
218 #include <limits>
219 
220 // fix for missing __builtin_clog functions
221 #if defined (__INTEL_COMPILER)
222 # if defined (_GLIBCXX_USE_C99_COMPLEX)
223 # undef _GLIBCXX_USE_C99_COMPLEX
224 # define _GLIBCXX_USE_C99_COMPLEX 0
225 # endif
226 #endif
227 
228 #include <complex>
229 #include <algorithm>
230 #include <exception>
231 #include <utility>
232 #include <new>
233 
234 #if !defined(CVM_NO_DEFAULT_RANDOM_ENGINE_SUPPORTED)
235 # include <random>
236 #endif
237 
238 #if defined (STLPORT)
239 # define CVM_USES_STLPORT
240 #endif
241 
242 #if defined (_DEBUG) || defined (DEBUG)
243 # define CVM_DEBUG
244 # include <assert.h>
245 # define CVM_ASSERT(p,n) _cvm_assert(p,n);
246 #else
247 # define CVM_ASSERT(p,n)
248 #endif
249 
250 #define CVM_OK 0
251 #define CVM_OUTOFMEMORY 1
252 #define CVM_WRONGSIZE 5
253 #define CVM_SIZESMISMATCH 6
254 #define CVM_WRONGMKLARG 7
255 #define CVM_WRONGMKLARG2 8
256 #define CVM_SINGULARMATRIX 9
257 #define CVM_NOTPOSITIVEDEFINITE 10
258 #define CVM_WRONGCHOLESKYFACTOR 11
259 #define CVM_WRONGBUNCHKAUFMANFACTOR 12
260 #define CVM_NOTPOSITIVEDIAG 13
261 #define CVM_CONVERGENCE_ERROR 14
262 #define CVM_DIVISIONBYZERO 15
263 #define CVM_SEMAPHOREERROR 16
264 #define CVM_READ_ONLY_ACCESS 17
265 #define CVM_SUBMATRIXACCESSERROR 18
266 #define CVM_SUBMATRIXNOTAVAILABLE 19
267 #define CVM_MATRIXNOTSYMMETRIC 20
268 #define CVM_MATRIXNOTHERMITIAN 21
269 #define CVM_BREAKS_HERMITIANITY 22
270 #define CVM_METHODNOTAVAILABLE 23
271 #define CVM_NOTIMPLEMENTED 24
272 #define CVM_CANT_RESIZE_SHARED_MEM 25
273 #define CVM_NOT_CONJUGATED 26
274 // 8.1
275 #define CVM_WRONGSIZE_LT 27
276 #define CVM_WRONGSIZE_LE 28
277 #define CVM_INDEX_GT 29
278 #define CVM_INDEX_GE 30
279 #define CVM_OUTOFRANGE_LTGT 31
280 #define CVM_OUTOFRANGE_LTGE 32
281 #define CVM_OUTOFRANGE_LTGE1 33
282 #define CVM_OUTOFRANGE_LTGE2 34
283 #define CVM_SIZESMISMATCH_GT 35
284 #define CVM_SIZESMISMATCH_LT 36
285 
286 #define CVM_THE_LAST_ERROR_CODE 37
287 #define CVM_MATRIX_ELEMENT_SEPARATOR " "
288 #define CVM_EOL std::endl
289 
290 typedef unsigned char tbyte;
291 
292 #ifdef CVM_NO_NAMESPACE
293 # define CVM_NAMESPACE_BEG
294 # define CVM_NAMESPACE_END
295 #else
296 # define CVM_NAMESPACE_BEG namespace cvm {
297 # define CVM_NAMESPACE_END }
298 #endif
299 
300 
302 
304 template<typename T>
306 public:
307  void operator () (T* d) const {
308  if (d != nullptr) {
309  ::delete[] d;
310  }
311  }
312 };
313 
314 template<typename TR, typename T> class basic_array;
315 template<typename TR, typename TC> class Matrix;
316 template<typename TR, typename TC> class SqMatrix;
317 template<typename TR, typename TC> class BandMatrix;
318 template<typename TR> class basic_rvector;
319 template<typename TR> class basic_rmatrix;
320 template<typename TR> class basic_srmatrix;
321 template<typename TR, typename TC> class basic_cvector;
322 template<typename TR, typename TC> class basic_cmatrix;
323 template<typename TR, typename TC> class basic_scmatrix;
324 template<typename TR> class basic_srbmatrix;
325 template<typename TR, typename TC> class basic_scbmatrix;
326 template<typename TR> class basic_srsmatrix;
327 template<typename TR, typename TC> class basic_schmatrix;
328 template<typename T, typename TR> class type_proxy;
329 
332 {
333  // message string maps
334  typedef std::map<int, std::string, std::less<int> > map_Msg;
335  typedef map_Msg::iterator itr_Msg;
336  typedef map_Msg::const_iterator citr_Msg;
337  typedef std::pair<int, std::string> pair_Msg;
338 
339 private:
340  std::string msUnknown;
341  map_Msg mmMsg;
342  CVM_API ErrMessages();
343 
344 public:
345  CVM_API const std::string& _get(int nException);
346  CVM_API bool _add(int nNewCause, const char* szNewMessage);
347 
348  static CVM_API ErrMessages& ErrMessagesInstance();
350 };
351 
352 
358 class cvmexception : public std::exception
359 {
360 protected:
361  int mnCause;
362  mutable char mszMsg[256];
363 
365  virtual const char* _get_message(int nCause) const {
366  return ErrMessages::ErrMessagesInstance()._get(nCause).c_str();
367  }
369 
370 public:
377  : mnCause(CVM_OK)
378  {
379  mszMsg[0] = '\0';
380  }
381 
389  CVM_API explicit cvmexception(int nCause, ...); // TODO
390 
394  CVM_API cvmexception(const cvmexception& e);
395 
399  virtual ~cvmexception() throw() { }
400 
409  int cause() const {
410  return mnCause;
411  }
412 
420  const char* what() const throw() override {
421  return mszMsg;
422  }
423 
431  static int getNextCause() {
433  }
434 
444  static bool add(int nNewCause, const char* szNewMessage) {
445  return ErrMessages::ErrMessagesInstance()._add(nNewCause, szNewMessage);
446  }
447 };
448 
450 template<typename T>
451 CVM_API void __copy(tint nSize, const T* pFrom, tint nFromIncr, T* pTo, tint nToIncr);
452 template<typename T>
453 CVM_API void __swap(tint nSize, T* p1, tint n1Incr, T* p2, tint n2Incr);
454 
455 template<typename TC, typename TR>
456 CVM_API TR _real(const TC& mT);
457 template<typename TC, typename TR>
458 CVM_API TR _imag(const TC& mT);
459 
460 template<typename TR>
461 CVM_API TR __dot(const TR* pDm, tint mnSize, tint mnIncr, const TR* pd, tint nIncr);
462 template<typename TC>
463 CVM_API TC __dotu(const TC* pDm, tint mnSize, tint mnIncr, const TC* pd, tint nIncr);
464 template<typename TC>
465 CVM_API TC __dotc(const TC* pDm, tint mnSize, tint mnIncr, const TC* pd, tint nIncr);
466 
467 template<typename TR, typename TC>
468 CVM_API TR __norm(const TC* pd, tint nSize, tint nIncr);
469 
470 template<typename TC>
471 CVM_API tint __idamax(const TC* pd, tint nSize, tint nIncr);
472 template<typename TC>
473 CVM_API tint __idamin(const TC* pd, tint nSize, tint nIncr);
474 
475 template<typename TC>
476 CVM_API void __add(TC* pDm, tint mnSize, tint mnIncr, const TC* pv, tint nIncr);
477 template<typename TC>
478 CVM_API void __subtract(TC* pDm, tint mnSize, tint mnIncr, const TC* pv, tint nIncr);
480 template<typename TR, typename TC>
481 CVM_API void __scal(TC* pDm, tint mnSize, tint mnIncr, TR dScal);
482 template<typename TC>
483 CVM_API void __conj(TC* pDm, tint mnSize, tint mnIncr);
484 
485 template<typename TR, typename TC>
486 CVM_API void __copy_real(TC* pDm, tint mnSize, tint mnIncr, const TR* pRe, tint nReIncr);
487 template<typename TR, typename TC>
488 CVM_API void __copy_imag(TC* pDm, tint mnSize, tint mnIncr, const TR* pRe, tint nReIncr);
489 template<typename TR, typename TC>
490 CVM_API void __copy2(TC* pDm, tint mnSize, tint mnIncr, const TR* pRe, const TR* pIm, tint nReIncr = 1, tint nImIncr = 1);
491 
492 template<typename TC, typename TM, typename TV>
493 CVM_API void __gemv(bool bLeft, const TM& m, TC dAlpha, const TV& v, TC dBeta, TV& vRes);
494 template<typename TC, typename TM, typename TV>
495 CVM_API void __gbmv(bool bLeft, const TM& m, TC dAlpha, const TV& v, TC dBeta, TV& vRes);
496 template<typename TR, typename TM, typename TV>
497 CVM_API void __symv(const TM& m, TR dAlpha, const TV& v, TR dBeta, TV& vRes);
498 template<typename TC, typename TM, typename TV>
499 CVM_API void __shmv(const TM& m, TC dAlpha, const TV& v, TC dBeta, TV& vRes);
500 template<typename TC, typename TM>
501 CVM_API void __gemm(const TM& ml, bool bTrans1, const TM& mr, bool bTrans2, TC dAlpha, TM& mRes, TC dBeta);
502 template<typename TR, typename TSM, typename TM>
503 CVM_API void __symm(bool bLeft, const TSM& ml, const TM& mr, TR dAlpha, TM& mRes, TR dBeta);
504 template<typename TC, typename TSM, typename TM>
505 CVM_API void __hemm(bool bLeft, const TSM& ml, const TM& mr, TC dAlpha, TM& mRes, TC dBeta);
506 
507 template<typename TC, typename TV>
508 CVM_API void __polynom(TC* pDm, tint ldP, tint mnM, const TC* pd, tint ldA, const TV& v);
509 template<typename T>
510 CVM_API void __inv(T& m, const T& mArg) throw(cvmexception);
511 template<typename T, typename TR>
512 CVM_API void __exp(T& m, const T& mArg, TR tol) throw(cvmexception);
513 template<typename T, typename TR>
514 CVM_API void __exp_symm(T& m, const T& mArg, TR tol) throw(cvmexception);
515 template<typename TR, typename TC, typename TRM>
516 CVM_API void __solve(const TRM& m, tint nrhs, const TC* pB, tint ldB, TC* pX, tint ldX, TR& dErr,
517  const TC* pLU, const tint* pPivots, int transp_mode) throw(cvmexception);
518 template<typename TC, typename TM, typename TSM>
519 CVM_API void __svd(TC* pd, tint nSize, tint nIncr, const TM& mArg, TSM* mU, TSM* mVH) throw(cvmexception);
520 template<typename TR, typename TM, typename TX>
521 CVM_API void __pinv(TX& mX, const TM& mArg, TR threshold) throw(cvmexception);
522 template<typename TV, typename TSM, typename TSCM>
523 CVM_API void __eig(TV& vRes, const TSM& mArg, TSCM* mEigVect, bool bRightVect) throw(cvmexception);
524 template<typename TR, typename TM>
525 CVM_API void __cond_num(const TM& mArg, TR& dCondNum) throw(cvmexception);
526 template<typename TM>
527 void __low_up(TM& m, tint* nPivots) throw(cvmexception);
528 
529 template<typename TR>
530 CVM_API void __randomize(TR* pDm, tint mnSize, tint mnIncr, TR dFrom, TR dTo);
531 template<typename TC, typename TR>
532 CVM_API void __randomize_real(TC* pDm, tint mnSize, tint mnIncr, TR dFrom, TR dTo);
533 template<typename TC, typename TR>
534 CVM_API void __randomize_imag(TC* pDm, tint mnSize, tint mnIncr, TR dFrom, TR dTo);
535 
536 template<typename TR, typename TM, typename TV>
537 CVM_API void __ger(TM& m, const TV& vCol, const TV& vRow, TR dAlpha);
538 template<typename TC, typename TM, typename TV>
539 CVM_API void __geru(TM& m, const TV& vCol, const TV& vRow, TC cAlpha);
540 template<typename TC, typename TM, typename TV>
541 CVM_API void __gerc(TM& m, const TV& vCol, const TV& vRow, TC cAlpha);
542 
543 template<typename TC, typename TSM>
544 CVM_API void __syrk(bool bTransp, TC alpha, tint k, const TC* pA, tint ldA, TC beta, TSM& m);
545 template<typename TC, typename TSM>
546 CVM_API void __syr2k(bool bTransp, TC alpha, tint k, const TC* pA, tint ldA, const TC* pB, tint ldB, TC beta, TSM& m);
547 template<typename TR, typename TC, typename TSM>
548 CVM_API void __herk(bool bTransp, TR alpha, tint k, const TC* pA, tint ldA, TR beta, TSM& m);
549 template<typename TR, typename TC, typename TSM>
550 CVM_API void __her2k(bool bTransp, TC alpha, tint k, const TC* pA, tint ldA, const TC* pB, tint ldB, TR beta, TSM& m);
551 
552 template<typename TM>
553 CVM_API tint __cholesky(TM& m);
554 template<typename TM>
555 CVM_API void __bunch_kaufman(TM& m, tint* nPivots) throw(cvmexception);
556 template<typename TR, typename TM, typename TV>
557 CVM_API void __poequ(const TM& m, TV& vScalings, TR& dCond, TR& dMax);
558 
559 template<typename TM, typename TSM>
560 CVM_API void __qre(const TM& mA, TM& mQ, TSM& mR) throw(cvmexception);
561 template<typename TM, typename TSM>
562 CVM_API void __qrf(const TM& mA, TSM& mQ, TM& mR) throw(cvmexception);
563 template<typename TM, typename TSM>
564 CVM_API void __rqe(const TM& mA, TSM& mR, TM& mQ) throw(cvmexception);
565 template<typename TM, typename TSM>
566 CVM_API void __rqf(const TM& mA, TM& mR, TSM& mQ) throw(cvmexception);
567 
568 template<typename TM, typename TSM>
569 CVM_API void __lqe(const TM& mA, TSM& mL, TM& mQ) throw(cvmexception);
570 template<typename TM, typename TSM>
571 CVM_API void __lqf(const TM& mA, TM& mL, TSM& mQ) throw(cvmexception);
572 template<typename TM, typename TSM>
573 CVM_API void __qle(const TM& mA, TM& mQ, TSM& mL) throw(cvmexception);
574 template<typename TM, typename TSM>
575 CVM_API void __qlf(const TM& mA, TSM& mQ, TM& mL) throw(cvmexception);
576 
577 // Linear Least Squares problems
578 template<typename TM, typename TV>
579 CVM_API void __gels(bool transpose, TM& mA, const TM& mB, TM& mX, TV& vErr) throw(cvmexception);
580 template<typename TR, typename TM>
581 CVM_API void __gelsy(TM& mA, const TM& mB, TM& mX, TR tol, tint& rank) throw(cvmexception);
582 template<typename TR, typename TV, typename TM>
583 CVM_API void __gelss(TM& mA, const TM& mB, TM& mX, TR rcond, TV& vSV, tint& rank) throw(cvmexception);
584 template<typename TR, typename TV, typename TM>
585 CVM_API void __gelsd(TM& mA, const TM& mB, TM& mX, TR rcond, TV& vSV, tint& rank) throw(cvmexception);
586 
587 // 8.1 Generalized Eigenvalue Problem
588 template<typename TSRM, typename TSCM, typename TRV, typename TCV>
589 CVM_API void __ggev(TSRM& mA, TSRM& mB, TCV& vAlpha, TRV& vBeta,
590  TSCM* mEigVectLeft, TSCM* mEigVectRight) throw (cvmexception);
591 
592 
593 template<typename T>
594 inline const T& _cvm_min(const T& x, const T& y) {
595  return x < y ? x : y;
596 }
597 
598 template<typename T>
599 inline const T& _cvm_max(const T& x, const T& y) {
600  return x > y ? x : y;
601 }
602 
603 template<class TR>
604 inline const TR* __get_real_p(const std::complex<TR>* c)
605 {
606 #if defined (CVM_USES_STLPORT)
607  return &c->_M_re;
608 #else
609  return reinterpret_cast<const TR*>(c);
610 #endif
611 }
612 
613 template<class TR>
614 inline const TR* __get_imag_p(const std::complex<TR>* c)
615 {
616 #if defined (CVM_USES_STLPORT)
617  return &c->_M_im;
618 #else
619  return reinterpret_cast<const TR*>(c) + 1;
620 #endif
621 }
622 
623 template<class TR>
624 inline TR* __get_real_p(std::complex<TR>* c)
625 {
626 #if defined (CVM_USES_STLPORT)
627  return &c->_M_re;
628 #else
629  return reinterpret_cast<TR*>(c);
630 #endif
631 }
632 
633 template<class TR>
634 inline TR* __get_imag_p(std::complex<TR>* c)
635 {
636 #if defined (CVM_USES_STLPORT)
637  return &c->_M_im;
638 #else
639  return reinterpret_cast<TR*>(c) + 1;
640 #endif
641 }
642 
643 template<class TR> inline const TR* __get_real_p(const TR* d) { return d; }
644 template<class TR> inline const TR* __get_imag_p(const TR* d) { return d; }
645 template<class TR> inline TR* __get_real_p(TR* d) { return d; }
646 template<class TR> inline TR* __get_imag_p(TR* d) { return d; }
647 
649 class Chars
650 {
651 protected:
652  static char mchars[15];
653 
654 public:
655  static const char* pT() {return mchars;}
656  static const char* pN() {return mchars + 1;}
657  static const char* pU() {return mchars + 2;}
658  static const char* pL() {return mchars + 3;}
659  static const char* pP() {return mchars + 4;}
660  static const char* pQ() {return mchars + 5;}
661  static const char* pB() {return mchars + 6;}
662  static const char* pE() {return mchars + 7;}
663  static const char* pR() {return mchars + 8;}
664  static const char* pA() {return mchars + 9;}
665  static const char* pS() {return mchars + 10;}
666  static const char* pV() {return mchars + 11;}
667  static const char* pO() {return mchars + 12;}
668  static const char* pI() {return mchars + 13;}
669  static const char* pC() {return mchars + 14;}
670 };
671 
672 template<class TR>
673 inline TR basic_cvmMachMin() {
674  static const TR _min = (std::numeric_limits<TR>::min)();
675  return _min;
676 }
677 
678 template<class TR>
679 inline TR basic_cvmMachSp() {
680  static const TR _eps = (std::numeric_limits<TR>::epsilon)();
681  return _eps;
682 }
683 
684 template<typename TR>
685 inline TR basic_cvmMachMax() {
686  static const TR _max((std::numeric_limits<TR>::max)());
687  return _max;
688 }
689 
690 template<typename TR>
691 inline TR basic_cvmLogMachMax() {
692  static const TR _log_max(::log(basic_cvmMachMax<TR>()));
693  return _log_max;
694 }
695 
696 // range and index checkers
697 template<typename T>
698 void _check_negative(int err_code, T val) throw(cvmexception)
699 {
700  static T zero = T(0);
701  if (val < zero) {
702  throw cvmexception(err_code, val);
703  }
704 }
705 
706 template<typename T>
707 void _check_positive(int err_code, T val) throw(cvmexception)
708 {
709  static T zero = T(0);
710  if (val > zero) {
711  throw cvmexception(err_code, val);
712  }
713 }
714 
715 template<typename T>
716 void _check_positive(int err_code, T val, const char* func, const char* file, int line) throw(cvmexception)
717 {
718  static T zero = T(0);
719  if (val > zero) {
720  throw cvmexception(err_code, func, file, line);
721  }
722 }
723 
724 template<typename T>
725 void _check_lt(int err_code, T idx, T limit) throw(cvmexception)
726 {
727  if (idx < limit) {
728  throw cvmexception(err_code, idx, limit);
729  }
730 }
731 
732 template<typename T>
733 void _check_le(int err_code, T idx, T limit) throw(cvmexception)
734 {
735  if (idx <= limit) {
736  throw cvmexception(err_code, idx, limit);
737  }
738 }
739 
740 template<typename T>
741 void _check_ge(int err_code, T idx, T limit) throw(cvmexception)
742 {
743  if (idx >= limit) {
744  throw cvmexception(err_code, idx, limit);
745  }
746 }
747 
748 template<typename T>
749 void _check_gt(int err_code, T idx, T limit) throw(cvmexception)
750 {
751  if (idx > limit) {
752  throw cvmexception(err_code, idx, limit);
753  }
754 }
755 
756 template<typename T>
757 void _check_ne(int err_code, T idx, T limit) throw(cvmexception)
758 {
759  if (idx != limit) {
760  throw cvmexception(err_code, idx, limit);
761  }
762 }
763 
764 template<typename T>
765 void _check_lt_gt(int err_code, T idx, T low_limit, T up_limit) throw(cvmexception)
766 {
767  if (idx < low_limit || idx > up_limit) {
768  throw cvmexception(err_code, idx, low_limit, up_limit);
769  }
770 }
771 
772 template<typename T>
773 void _check_lt_ge(int err_code, T idx, T low_limit, T up_limit) throw(cvmexception)
774 {
775  if (idx < low_limit || idx >= up_limit) {
776  throw cvmexception(err_code, idx, low_limit, up_limit);
777  }
778 }
780 
781 
783 #ifdef CVM_USE_POOL_MANAGER
784 
786 class CVM_API MemoryBlocks
787 {
788  struct BlockProperty
789  {
790  size_t mnSize;
791  tint mnRefCount;
792  BlockProperty(size_t nSize, tint nRefCount) : mnSize(nSize), mnRefCount(nRefCount) { }
793  };
794 
795  typedef std::map<tbyte*, BlockProperty, std::less<tbyte*> > map_Blocks;
796  typedef map_Blocks::iterator itr_Blocks;
797 
798  typedef std::multimap<size_t, tbyte*> map_FreeBs;
799  typedef map_FreeBs::iterator itr_FreeBs;
800 
801  typedef std::map<tbyte*, itr_FreeBs, std::less<tbyte*> > map_FreeIt;
802  typedef map_FreeIt::iterator itr_FreeIt;
803 
804  map_FreeBs mFreeBs;
805  map_FreeIt mFreeIt;
806 
807  map_Blocks mBlocks;
808 
809 
810 public:
811  void AddBlock(tbyte* pBlock, size_t nBytes, bool bOccupied);
812  tbyte* GetFreeBlock(size_t nBytes);
813  void AddPair(tbyte* pBlock, size_t nBytes, size_t nRest);
814 #ifdef CVM_DEBUG
815  void Assert(const void* pvBlock, size_t nBytes);
816 #endif
817  void AddNew(tbyte* pBlock, size_t nBytes);
818  tbyte* AddRef(const tbyte* pBlock);
819  tint FreeBlock(tbyte* pBlock);
820 };
821 
822 
824 class CVM_API MemoryPool
825 {
826  typedef std::list<tbyte*> list_blocks;
827 
828  struct DeletePtr {
829  template<class T>
830  void operator () (T* p) const {
831  ::delete[] p;
832  }
833  };
834 
835  list_blocks mOutBlocks;
836  MemoryBlocks mMemoryBlocks;
837 
838 public:
839  MemoryPool();
840  ~MemoryPool();
841 
842  tbyte* Malloc(size_t nBytes) throw(cvmexception);
843  tbyte* AddRef(const tbyte* pd);
844  tint Free(tbyte*& pToFree) throw(cvmexception);
845 
846 #ifdef CVM_DEBUG
847  void Assert(const void* pvBlock, size_t nBytes) {
848  mMemoryBlocks.Assert(pvBlock, nBytes);
849  }
850 #endif
851  void Clear();
852 };
853 
854 CVM_API tbyte* _cvmMalloc(size_t nBytes) throw(cvmexception);
855 CVM_API tbyte* _cvmAddRef(const tbyte* pd);
856 CVM_API tint _cvmFree(tbyte*& pd);
857 
864 template<typename T>
865 inline T* cvmMalloc(size_t nEls) throw(cvmexception) {
866  return reinterpret_cast<T*>(_cvmMalloc(nEls * sizeof(T)));
867 }
868 
875 template<typename T>
876 inline T* cvmAddRef(const T* pd) {
877  return reinterpret_cast<T*>(_cvmAddRef(reinterpret_cast<const tbyte*>(pd)));
878 }
879 
887 template<typename T>
888 inline tint cvmFree(T*& pd) {
889  return _cvmFree(reinterpret_cast<tbyte*&>(pd));
890 }
891 
892 #else // CVM_USE_POOL_MANAGER
893 
901 template<typename T>
902 inline T* cvmMalloc(tint nEls) throw(cvmexception)
903 {
904  _check_lt(CVM_WRONGSIZE_LT, nEls, TINT_ZERO);
905  return nEls > 0 ? ::new T[nEls] : nullptr;
906 }
907 
914 template<typename T>
915 inline tint cvmFree(T*& pd)
916 {
917  if (pd != nullptr) {
918  ::delete[] pd;
919  pd = nullptr;
920  }
921  return 0;
922 }
923 
924 #endif // !CVM_USE_POOL_MANAGER
925 
932 CVM_API void cvmExit();
933 
942 template<typename T>
943 inline void cvmZeroMemory(T* p, tint nEls) {
944  memset(p, 0, nEls * sizeof(T));
945 }
946 
948 CVM_API void _cvm_assert(const void* pvBlock, size_t nBytes);
949 
950 inline float _abs(const float& v) {
951  return static_cast<float>(fabs(v));
952 }
953 
954 inline double _abs(const double& v) {
955  return fabs(v);
956 }
957 
958 inline long double _abs(const long double& v) {
959  return fabs(v);
960 }
961 
962 inline float _abs(const std::complex<float>& v) {
963  return static_cast<float>(sqrt(v.real() * v.real() + v.imag() * v.imag()));
964 }
965 
966 inline double _abs(const std::complex<double>& v) {
967  return sqrt(v.real() * v.real() + v.imag() * v.imag());
968 }
969 
970 inline long double _abs(const std::complex<long double>& v) {
971  return sqrt(v.real() * v.real() + v.imag() * v.imag());
972 }
973 
974 inline tint _abs(const tint& v) {
975  return v < 0 ? -v : v;
976 }
977 
978 inline float _sqrt(const float& v) {
979  return static_cast<float>(sqrt(v));
980 }
981 
982 inline double _sqrt(const double& v) {
983  return sqrt(v);
984 }
985 
986 inline std::complex<float> _conjugate(const std::complex<float>& v) {
987  return std::complex<float>(v.real(), -v.imag());
988 }
989 
990 inline std::complex<double> _conjugate(const std::complex<double>& v) {
991  return std::complex<double>(v.real(), -v.imag());
992 }
993 
994 template<typename TR>
995 inline bool _conjugated(const std::complex<TR>& v1, const std::complex<TR>& v2, TR tol)
996 {
997  return _abs(v1.real() - v2.real()) <= tol &&
998  _abs(v1.imag() + v2.imag()) <= tol;
999 }
1001 
1002 
1004 template<typename T, typename TR>
1005 inline std::ostream& operator << (std::ostream& os, const type_proxy<T,TR>& mOut)
1006 {
1007 #if defined (_MSC_VER) && !defined (CVM_USES_STLPORT)
1008  os.imbue(std::locale::empty());
1009 #endif
1010  os << static_cast<T>(mOut);
1011  return os;
1012 }
1013 
1015 template<typename T, typename TR>
1016 inline std::istream& operator >> (std::istream& is, type_proxy<T,TR>& v)
1017 {
1018 #if defined (_MSC_VER) && !defined (CVM_USES_STLPORT)
1019  is.imbue(std::locale::empty());
1020 #endif
1021  T t;
1022  is >> t;
1023  v = t;
1024  return is;
1025 }
1026 
1060 template<typename TR, typename TC>
1061 std::istream& operator >> (std::istream& is, basic_array<TR,TC>& aIn)
1062 {
1063 #if !defined (CVM_USES_STLPORT) && defined (_MSC_VER)
1064  is.imbue(std::locale::empty());
1065 #endif
1066  const tint nSize = aIn.size() * aIn.incr();
1067  CVM_ASSERT(aIn.get(), ((aIn.size() - 1) * aIn.incr() + 1) * sizeof(TC))
1068  for (tint i = 0; i < nSize && is.good(); i += aIn.incr()) {
1069  is >> aIn.get()[i];
1070  }
1071  return is;
1072 }
1073 
1097 template<typename TR, typename TC>
1098 std::ostream& operator << (std::ostream& os, const basic_array<TR,TC>& aOut)
1099 {
1100 #if !defined (CVM_USES_STLPORT) && defined (_MSC_VER)
1101  os.imbue(std::locale::empty());
1102 #endif
1103  const tint nSize = aOut.size() * aOut.incr();
1104  CVM_ASSERT(aOut.get(), ((aOut.size() - 1) * aOut.incr() + 1) * sizeof(TC))
1105  for (tint i = 0; i < nSize && os.good(); i += aOut.incr()) {
1106  os << aOut.get()[i] << CVM_MATRIX_ELEMENT_SEPARATOR;
1107  }
1108  os << CVM_EOL;
1109  return os;
1110 }
1111 
1150 template<typename TR, typename TC>
1151 std::istream& operator >> (std::istream& is, Matrix<TR,TC>& mIn)
1152 {
1153 #if defined (_MSC_VER) && !defined (CVM_USES_STLPORT)
1154  is.imbue(std::locale::empty());
1155 #endif
1156  TC v;
1157  for (tint i = 0; i < mIn.msize(); ++i) {
1158  for (tint j = 0; j < mIn.nsize() && is.good(); ++j) {
1159  is >> v;
1160  mIn._ij_proxy_val(i, j) = v;
1161  }
1162  }
1163  return is;
1164 }
1165 
1193 template<typename TR, typename TC>
1194 std::ostream& operator << (std::ostream& os, const Matrix<TR,TC>& mOut)
1195 {
1196 #if defined (_MSC_VER) && !defined (CVM_USES_STLPORT)
1197  os.imbue(std::locale::empty());
1198 #endif
1199  for (tint i = 0; i < mOut.msize(); ++i) {
1200  for (tint j = 0; j < mOut.nsize() && os.good(); ++j) {
1201  os << mOut._ij_val(i, j) << CVM_MATRIX_ELEMENT_SEPARATOR;
1202  }
1203  os << CVM_EOL;
1204  }
1205  return os;
1206 }
1207 
1209 template<typename TR, typename TC, typename RM, typename RBM>
1210 inline void _copy_b_matrix(RM& m, RBM& mb, bool bLeftToRight)
1211 {
1212  const tint nM = mb.msize();
1213  const tint nN = mb.nsize();
1214  const tint nKL = mb.lsize();
1215  const tint nKU = mb.usize();
1216  const tint nCol = 1 + nKL + nKU;
1217  tint nS, nShiftL, nShiftR;
1218  TC* pL;
1219  TC* pR;
1220 
1221  for (tint i = 0; i < nN; ++i) {
1222  nS = nCol;
1223  nShiftL = 0;
1224  nShiftR = 0;
1225  if (i < nKU) {
1226  nShiftR = nKU - i;
1227  nS -= nShiftR;
1228  } else {
1229  nShiftL = i - nKU;
1230  }
1231  if (nN - i <= nKL) {
1232  nS -= nKL + 1 - (nN - i);
1233  }
1234 
1235  pL = m.get() + i * nM + nShiftL;
1236  pR = mb._pb() + i * nCol + nShiftR;
1237 
1238  __copy<TC>(nS,
1239  bLeftToRight ? pL : pR,
1240  1,
1241  bLeftToRight ? pR : pL,
1242  1);
1243  }
1244 }
1245 
1246 // pd = p1 + p2
1247 template<typename TR, typename TC>
1248 inline void _sum(TC* pd, tint nSize, tint nIncr,
1249  const TC* p1, tint nIncr1, const TC* p2, tint nIncr2)
1250 {
1251  if (pd == p1) {
1252  if (pd == p2) {
1253  static const TR two(2.);
1254  __scal<TR,TC>(pd, nSize, nIncr, two);
1255  } else {
1256  __add<TC>(pd, nSize, nIncr, p2, nIncr2);
1257  }
1258  } else {
1259  if (pd == p2) {
1260  __add<TC>(pd, nSize, nIncr, p1, nIncr1);
1261  } else {
1262  __copy<TC>(nSize, p1, nIncr1, pd, nIncr);
1263  __add<TC>(pd, nSize, nIncr, p2, nIncr2);
1264  }
1265  }
1266 }
1267 
1268 // pd = p1 - p2
1269 template<typename TR, typename TC>
1270 inline void _diff(TC* pd, tint nSize, tint nIncr,
1271  const TC* p1, tint nIncr1, const TC* p2, tint nIncr2)
1272 {
1273  if (pd == p1) {
1274  if (pd == p2) {
1275  static const TR zero(0.);
1276  __scal<TR,TC>(pd, nSize, nIncr, zero);
1277  } else {
1278  __subtract<TC>(pd, nSize, nIncr, p2, nIncr2);
1279  }
1280  } else {
1281  if (pd == p2) {
1282  static const TR mone(-1.);
1283  __subtract<TC>(pd, nSize, nIncr, p1, nIncr1);
1284  __scal<TR,TC>(pd, nSize, nIncr, mone);
1285  } else {
1286  __copy<TC>(nSize, p1, nIncr1, pd, nIncr);
1287  __subtract<TC>(pd, nSize, nIncr, p2, nIncr2);
1288  }
1289  }
1290 }
1291 
1292 // pDm = pDm + pd
1293 template<typename TR, typename TC>
1294 inline void _incr(TC* pDm, tint mnSize, tint mnIncr, const TC* pd, tint nIncr)
1295 {
1296  if (pDm == pd) {
1297  static const TR two(2.);
1298  __scal<TR,TC>(pDm, mnSize, mnIncr, two);
1299  } else {
1300  __add<TC>(pDm, mnSize, mnIncr, pd, nIncr);
1301  }
1302 }
1303 
1304 // pDm = pDm - pd
1305 template<typename TR, typename TC>
1306 inline void _decr(TC* pDm, tint mnSize, tint mnIncr, const TC* pd, tint nIncr)
1307 {
1308  if (pDm == pd) {
1309  static const TR zero(0.);
1310  __scal<TR,TC>(pDm, mnSize, mnIncr, zero);
1311  } else {
1312  __subtract<TC>(pDm, mnSize, mnIncr, pd, nIncr);
1313  }
1314 }
1315 
1316 // fills real part
1317 template<typename TR, typename TC>
1318 void _set_real(TC* pDm, tint mnSize, tint mnIncr, TR d)
1319 {
1320  const tint nIncr2 = 2 * mnIncr;
1321  const tint nSize = mnSize * nIncr2;
1322  TR* pd = __get_real_p<TR>(pDm);
1323  for (tint i = 0; i < nSize; i += nIncr2) {
1324  CVM_ASSERT(pd, (i + 1) * sizeof(TR))
1325  pd[i] = d;
1326  }
1327 }
1328 
1329 // fills imaginary part
1330 template<typename TR, typename TC>
1331 void _set_imag(TC* pDm, tint mnSize, tint mnIncr, TR d)
1332 {
1333  const tint nIncr2 = 2 * mnIncr;
1334  const tint nSize = mnSize * nIncr2;
1335  TR* pd = __get_imag_p<TR>(pDm);
1336  for (tint i = 0; i < nSize; i += nIncr2) {
1337  CVM_ASSERT(pd, (i + 1) * sizeof(TR))
1338  pd[i] = d;
1339  }
1340 }
1342 
1343 
1349 template<typename T, typename TR>
1351 {
1352 protected:
1353  T& mT;
1354  bool mbReadOnly;
1355 
1356 public:
1364  type_proxy(T& ref, bool read_only)
1365  : mT(ref),
1366  mbReadOnly(read_only)
1367  {
1368  }
1369 
1377  type_proxy(const T& ref, bool read_only = true)
1378  : mT(const_cast<T&>(ref)),
1379  mbReadOnly(read_only)
1380  {
1381  }
1382 
1389  type_proxy(const type_proxy& p)
1390  : mT(p.mT),
1391  mbReadOnly(p.mbReadOnly)
1392  {
1393  }
1394 
1401  operator T() const {
1402  return mT;
1403  }
1404 
1405 /*
1406  doesn't work - masks the operator above
1407  operator const T& () const
1408  {
1409  return mT;
1410  }
1411 
1412  operator T& ()
1413  {
1414  if (mbReadOnly) throw cvmexception (CVM_READ_ONLY_ACCESS);
1415  return mT;
1416  }
1417 */
1418 
1425  T val() const {
1426  return mT;
1427  }
1428 
1435  T& get() throw(cvmexception) {
1436  if (mbReadOnly) throw cvmexception(CVM_READ_ONLY_ACCESS);
1437  return mT;
1438  }
1439 
1446  const T& get() const {
1447  return mT;
1448  }
1449 
1457  type_proxy& operator = (const type_proxy& p)
1458  {
1459  mT = p.mT;
1460  mbReadOnly = p.mbReadOnly;
1461  return *this;
1462  }
1463 
1471  type_proxy& operator = (const T& v) throw(cvmexception)
1472  {
1473  if (mbReadOnly) throw cvmexception(CVM_READ_ONLY_ACCESS);
1474  mT = v;
1475  return *this;
1476  }
1477 
1484  T* operator &() throw(cvmexception)
1485  {
1486  if (mbReadOnly) throw cvmexception(CVM_READ_ONLY_ACCESS);
1487  return &mT;
1488  }
1489 
1496  const T* operator &() const {
1497  return &mT;
1498  }
1499 
1507  template<typename U>
1508  T& operator += (const U& u) throw(cvmexception)
1509  {
1510  if (mbReadOnly) throw cvmexception(CVM_READ_ONLY_ACCESS);
1511  mT += T(u);
1512  return mT;
1513  }
1514 
1522  template<typename U>
1523  T& operator -= (const U& u) throw(cvmexception)
1524  {
1525  if (mbReadOnly) throw cvmexception(CVM_READ_ONLY_ACCESS);
1526  mT -= T(u);
1527  return mT;
1528  }
1529 
1537  template<typename U>
1538  T& operator *= (const U& u) throw(cvmexception)
1539  {
1540  if (mbReadOnly) throw cvmexception(CVM_READ_ONLY_ACCESS);
1541  mT *= T(u);
1542  return mT;
1543  }
1544 
1552  template<typename U>
1553  T& operator /= (const U& u) throw(cvmexception)
1554  {
1555  if (mbReadOnly) throw cvmexception(CVM_READ_ONLY_ACCESS);
1556  mT /= T(u);
1557  return mT;
1558  }
1559 
1567  template<typename U>
1568  T operator + (const U& u) const {
1569  return mT + T(u);
1570  }
1571 
1579  template<typename U>
1580  T operator - (const U& u) const {
1581  return mT - T(u);
1582  }
1583 
1591  template<typename U>
1592  T operator * (const U& u) const {
1593  return mT * T(u);
1594  }
1595 
1603  template<typename U>
1604  T operator / (const U& u) const {
1605  return mT / T(u);
1606  }
1607 
1615  T operator + (const type_proxy<T,TR>& u) const {
1616  return mT + u.mT;
1617  }
1618 
1626  T operator - (const type_proxy<T,TR>& u) const {
1627  return mT - u.mT;
1628  }
1629 
1637  T operator * (const type_proxy<T,TR>& u) const {
1638  return mT * u.mT;
1639  }
1640 
1648  T operator / (const type_proxy<T,TR>& u) const {
1649  return mT / u.mT;
1650  }
1651 
1658  T operator - () const {
1659  return -mT;
1660  }
1661 
1669  TR real() const {
1670  return _real<T,TR>(mT);
1671  }
1672 
1680  TR imag() const {
1681  return _imag<T,TR>(mT);
1682  }
1683 };
1684 
1686 template<typename T, typename TR>
1687 inline std::complex<TR> operator + (const std::complex<TR>& u, const type_proxy<T,TR>& p) {
1688  return u + p.val();
1689 }
1690 
1692 template<typename T, typename TR>
1693 inline std::complex<TR> operator - (const std::complex<TR>& u, const type_proxy<T,TR>& p) {
1694  return u - p.val();
1695 }
1696 
1698 template<typename T, typename TR>
1699 inline std::complex<TR> operator * (const std::complex<TR>& u, const type_proxy<T,TR>& p) {
1700  return u * p.val();
1701 }
1702 
1704 template<typename T, typename TR>
1705 inline std::complex<TR> operator / (const std::complex<TR>& u, const type_proxy<T,TR>& p) {
1706  return u / p.val();
1707 }
1708 
1710 template<typename U, typename T, typename TR>
1711 inline T operator + (U u, const type_proxy<T,TR>& p) {
1712  return T(u) + p.val();
1713 }
1714 
1716 template<typename U, typename T, typename TR>
1717 inline T operator - (U u, const type_proxy<T,TR>& p) {
1718  return T(u) - p.val();
1719 }
1720 
1722 template<typename U, typename T, typename TR>
1723 inline T operator * (U u, const type_proxy<T,TR>& p) {
1724  return T(u) * p.val();
1725 }
1726 
1728 template<typename U, typename T, typename TR>
1729 inline T operator / (U u, const type_proxy<T,TR>& p) {
1730  return T(u) / p.val();
1731 }
1732 
1733 
1735 template<typename T>
1737 {
1738  T mMax;
1739 #if !defined(CVM_NO_DEFAULT_RANDOM_ENGINE_SUPPORTED)
1740  std::random_device mre;
1741  std::uniform_int_distribution<int> mDist;
1742 #endif
1743 
1748  Randomizer() noexcept
1749 #if defined(CVM_NO_DEFAULT_RANDOM_ENGINE_SUPPORTED)
1750  : mMax(static_cast<T>(RAND_MAX))
1751  {
1752  srand(static_cast<unsigned int>(time(nullptr)));
1753  }
1754 #else
1755  : mMax(static_cast<T>((std::numeric_limits<int>::max)())),
1756  mre(),
1757  mDist(0, (std::numeric_limits<int>::max)())
1758  {
1759  }
1760 #endif
1761 
1766  T _get(T dFrom, T dTo) noexcept
1767  {
1768  const T dMin = _cvm_min<T>(dFrom, dTo);
1769  const T dMax = _cvm_max<T>(dFrom, dTo);
1770 #if defined(CVM_NO_DEFAULT_RANDOM_ENGINE_SUPPORTED)
1771  unsigned int nr = static_cast<unsigned int>(rand());
1772  return dMin + static_cast<T>(nr) * (dMax - dMin) / mMax;
1773 #else
1774  return dMin + static_cast<T>(mDist(mre)) * (dMax - dMin) / mMax;
1775 #endif
1776  }
1777 
1778 public:
1783  ~Randomizer() noexcept
1784  {
1785  }
1786 
1795  static T get(T dFrom, T dTo) noexcept
1796  {
1797  static Randomizer r;
1798  return r._get(dFrom, dTo);
1799  }
1800 };
1801 
1802 
1817 template<typename TR, typename TC>
1819 {
1820 protected:
1823 
1824 #ifdef CVM_USE_POOL_MANAGER
1825  TC* mpd;
1826 #else
1827  std::shared_ptr<TC> mp;
1828  TC* mpf;
1829 #endif
1830 
1831 public:
1832  typedef TC value_type;
1833  typedef value_type* pointer;
1835  typedef const value_type* const_iterator;
1836  typedef const value_type* const_pointer;
1838  typedef const value_type& const_reference;
1839  typedef size_t size_type;
1840  typedef ptrdiff_t difference_type;
1841  typedef std::reverse_iterator<const_iterator> const_reverse_iterator;
1842  typedef std::reverse_iterator<iterator> reverse_iterator;
1843 
1864  : msz(0),
1865  mincr(0),
1866 #ifdef CVM_USE_POOL_MANAGER
1867  mpd(nullptr)
1868 #else
1869  mp(),
1870  mpf(nullptr)
1871 #endif
1872  {
1873  }
1874 
1896  explicit basic_array(tint nSize, bool bZeroMemory = true)
1897  : msz(nSize), mincr(1),
1898 #ifdef CVM_USE_POOL_MANAGER
1899  mpd(cvmMalloc<TC>(static_cast<size_t>(msz)))
1900 #else
1901  mp(cvmMalloc<TC>(msz), ArrayDeleter<TC>()),
1902  mpf(nullptr)
1903 #endif
1904  {
1905  CVM_ASSERT(this->get(), msz * sizeof(TC))
1906  if (bZeroMemory) cvmZeroMemory<TC>(this->get(), msz);
1907  }
1908 
1936  basic_array(TC* pd, tint nSize, tint nIncr = 1)
1937  : msz(nSize), mincr(nIncr),
1938 #ifdef CVM_USE_POOL_MANAGER
1939  mpd(cvmAddRef<TC>(pd))
1940 #else
1941  mp(),
1942  mpf(pd)
1943 #endif
1944  {
1945  _check_le(CVM_WRONGSIZE_LE, msz, TINT_ZERO);
1946  CVM_ASSERT(this->get(), msz * sizeof(TC))
1947  }
1948 
1973  basic_array(const TC* pd, tint nSize, tint nIncr = 1)
1974  : msz(nSize), mincr(1),
1975 #ifdef CVM_USE_POOL_MANAGER
1976  mpd(cvmMalloc<TC>(static_cast<size_t>(msz)))
1977 #else
1978  mp(cvmMalloc<TC>(msz), ArrayDeleter<TC>()),
1979  mpf(nullptr)
1980 #endif
1981  {
1982  _check_le(CVM_WRONGSIZE_LE, msz, TINT_ZERO);
1983  CVM_ASSERT(this->get(), msz * sizeof(TC))
1984  __copy<TC>(msz, pd, nIncr, this->get(), this->incr());
1985  }
1986 
2010  basic_array(const TC* begin, const TC* end)
2011  : msz(tint(end - begin)), mincr(1),
2012 #ifdef CVM_USE_POOL_MANAGER
2013  mpd(cvmMalloc<TC>(static_cast<size_t>(msz)))
2014 #else
2015  mp(cvmMalloc<TC>(msz), ArrayDeleter<TC>()),
2016  mpf(nullptr)
2017 #endif
2018  {
2019  CVM_ASSERT(this->get(), msz * sizeof(TC))
2020  __copy<TC>(msz, begin, 1, this->get(), this->incr());
2021  }
2022 
2043  : msz(a.msz), mincr(1),
2044 #ifdef CVM_USE_POOL_MANAGER
2045  mpd(cvmMalloc<TC>(msz))
2046 #else
2047  mp(cvmMalloc<TC>(msz), ArrayDeleter<TC>()),
2048  mpf(nullptr)
2049 #endif
2050  {
2051  this->_assign(a.get(), a.incr());
2052  CVM_ASSERT(this->get(), msz * sizeof(TC))
2053  }
2054 
2072  basic_array(basic_array&& a) noexcept
2073  : msz(a.msz), mincr(a.mincr),
2074 #ifdef CVM_USE_POOL_MANAGER
2075  mpd(nullptr)
2076 #else
2077  mp(),
2078  mpf(nullptr)
2079 #endif
2080  {
2081  _move(std::move(a));
2082  CVM_ASSERT(this->get(), msz * sizeof(TC))
2083  }
2084 
2107  {
2108  _check_ne(CVM_SIZESMISMATCH, a.size(), this->size());
2109  this->_assign(a.get(), a.incr());
2110  return *this;
2111  }
2112 
2127  {
2128  _check_ne(CVM_SIZESMISMATCH, a.size(), this->size());
2129  _move(std::move(a));
2130  CVM_ASSERT(this->get(), msz * sizeof(TC))
2131  return *this;
2132  }
2133 
2139  virtual ~basic_array()
2140  {
2141 #ifdef CVM_USE_POOL_MANAGER
2142  cvmFree<TC>(mpd);
2143 #endif
2144  }
2145 
2164  tint size() const {
2165  return msz;
2166  }
2167 
2200  TC* get()
2201  {
2202 #ifdef CVM_USE_POOL_MANAGER
2203  return this->mpd;
2204 #else
2205  // if foreign pointer is set, this is it, shared one otherwise
2206  return this->mpf == nullptr ? this->mp.get() : this->mpf;
2207 #endif
2208  }
2209 
2242  const TC* get() const
2243  {
2244 #ifdef CVM_USE_POOL_MANAGER
2245  return this->mpd;
2246 #else
2247  // if foreign pointer is set, this is it, shared one otherwise
2248  return this->mpf == nullptr ? this->mp.get() : this->mpf;
2249 #endif
2250  }
2251 
2285  operator TC* () {
2286  return this->get();
2287  }
2288 
2322  operator const TC* () const {
2323  return this->get();
2324  }
2325 
2357  return this->at(static_cast<size_type>(n - CVM0));
2358  }
2359 
2390  TC operator () (tint n) const throw(cvmexception) {
2391  return this->at(static_cast<size_type>(n - CVM0));
2392  }
2393 
2425  return this->at(static_cast<size_type>(n - CVM0));
2426  }
2427 
2458  TC operator [] (tint n) const throw(cvmexception) {
2459  return this->at(static_cast<size_type>(n - CVM0));
2460  }
2461 
2484  basic_array& assign(const TC* p)
2485  {
2486  this->_assign(p, 1);
2487  return *this;
2488  }
2489 
2509  basic_array& set(TC x)
2510  {
2511  this->_set(x);
2512  return *this;
2513  }
2514 
2546  {
2547  this->_resize(nNewSize);
2548  return *this;
2549  }
2550 
2551 
2575  tint incr() const {
2576  return mincr;
2577  }
2578 
2605  tint indofmax() const {
2606  return this->_indofmax();
2607  }
2608 
2635  tint indofmin() const {
2636  return this->_indofmin();
2637  }
2638 
2674  virtual TR norm() const {
2675  return __norm<TR,TC>(this->get(), this->size(), this->incr());
2676  }
2677 
2711  virtual TR norminf() const
2712  {
2713  CVM_ASSERT(this->get(), ((this->_indofmax() - CVM0) * this->incr() + CVM0) * sizeof(TC))
2714  return _abs(this->get()[(this->_indofmax() - CVM0) * this->incr()]);
2715  }
2716 
2752  virtual TR norm1() const
2753  {
2754  TR dNorm(0);
2755  const tint nSize = this->size() * this->incr();
2756  for (tint i = 0; i < nSize; i += this->incr()) {
2757  dNorm += _abs(this->get()[i]);
2758  }
2759  return dNorm;
2760  }
2761 
2801  virtual TR norm2() const
2802  {
2803  return this->norm();
2804  }
2805 
2807  const tint* _pincr() const {
2808  return &this->mincr;
2809  }
2810 
2811  const tint* _psize() const {
2812  return &this->msz;
2813  }
2814 
2815  // this = this / d for real only
2816  void _div(TR d) throw(cvmexception)
2817  {
2818  static const TR one(1.);
2819  if (_abs(d) <= basic_cvmMachMin<TR>()) {
2821  }
2822  this->_scalr(one / d);
2823  }
2824 
2825  // "native pointer": to be redefined in classes with non-traditional storage, like band matrices etc.
2826  virtual TC* _pd() {
2827  return this->get();
2828  }
2829 
2830  // "native pointer": to be redefined in classes with non-traditional storage, like band matrices etc.
2831  virtual const TC* _pd() const {
2832  return this->get();
2833  }
2835 
2837  iterator begin() { return this->get();}
2839  const_iterator begin() const { return this->get();}
2841  iterator end() { return this->get() + this->size();}
2843  const_iterator end() const { return this->get() + this->size();}
2844 
2853 
2855  size_type max_size() const { return size_type(-1) / sizeof(TC);}
2857  size_type capacity() const { return static_cast<size_type>(this->size());}
2859  bool empty() const { return this->size() > 0;}
2860 
2862  reference front() { return *begin();}
2864  const_reference front() const { return *begin();}
2866  reference back() { return *(end() - 1);}
2868  const_reference back() const { return *(end() - 1);}
2869 
2870 // deprecated since 5.4.2 due to its incosistency with STL counterpart
2871 //
2872 // void reserve (size_type n) throw (cvmexception)
2873 // {
2874 // this->_resize(tint(n));
2875 // }
2876 
2878  void assign(size_type n, const TC& val) throw(cvmexception)
2879  {
2880  _check_gt(CVM_INDEX_GT, n, static_cast<size_type>(this->size()));
2881  CVM_ASSERT(this->get(), this->size() * sizeof(TC))
2882  for (tint i = 0; i < (tint)n; ++i) {
2883  this->get()[i] = val;
2884  }
2885  }
2886 
2889  {
2890  const tint n = end - begin;
2891  _check_gt(CVM_INDEX_GT, n, this->size());
2892  CVM_ASSERT(this->get(), this->size() * sizeof(TC))
2893  for (tint i = 0; i < n; ++i) {
2894  this->get()[i] = *(begin + i);
2895  }
2896  }
2897 
2899  void clear() {
2900  this->_resize(0);
2901  }
2902 
2905  {
2906  _check_ne(CVM_SIZESMISMATCH, v.size(), this->size());
2907  __swap<TC>(this->size(), this->get(), this->incr(), v.get(), v.incr());
2908  }
2909 
2912  {
2913  _check_ge(CVM_INDEX_GE, n, static_cast<size_type>(this->size()));
2914  CVM_ASSERT(this->get(), (n + 1) * sizeof(TC))
2915  return this->get()[n * this->incr()];
2916  }
2917 
2920  {
2921  _check_ge(CVM_INDEX_GE, n, static_cast<size_type>(this->size()));
2922  CVM_ASSERT(this->get(), (n + 1) * sizeof(TC))
2923  return this->get()[n * this->incr()];
2924  }
2925 
2953  void push_back(const TC& x) throw(cvmexception)
2954  {
2955  this->_resize(this->size() + 1);
2956  this->get()[this->size() - 1] = x;
2957  }
2958 
2985  void pop_back() throw(cvmexception) {
2986  if (this->size() > 0) this->_resize(this->size() - 1);
2987  }
2988 
3023  iterator insert(iterator position, const TC& x) throw(cvmexception)
3024  {
3025  const tint n = tint(position - this->begin());
3026  _check_lt_gt(CVM_OUTOFRANGE_LTGT, n, TINT_ZERO, this->size());
3027  CVM_ASSERT(this->get(), this->size() * sizeof(TC))
3028  this->_resize(this->size() + 1);
3029  for (tint i = this->size() - 1; i > n; --i) {
3030  this->get()[i] = this->get()[i - 1];
3031  }
3032  this->get()[n] = x;
3033  return iterator(this->get() + n);
3034  }
3035 
3070  {
3071  const tint n = tint(position - this->begin());
3072  _check_lt_gt(CVM_OUTOFRANGE_LTGT, n, TINT_ZERO, this->size());
3073 
3074  CVM_ASSERT(this->get(), this->size() * sizeof(TC))
3075  for (tint i = n; i < this->size() - 1; ++i) {
3076  this->get()[i] = this->get()[i + 1];
3077  }
3078  if (this->size() > 0) this->_resize(this->size() - 1);
3079  return iterator(this->get() + n);
3080  }
3081 
3082 protected:
3084  // internal protected constructor to be used in Matrix
3085  basic_array(tint sz, tint incr) noexcept
3086  : msz(sz), mincr(incr),
3087 #ifdef CVM_USE_POOL_MANAGER
3088  mpd(nullptr)
3089 #else
3090  mp(),
3091  mpf(nullptr)
3092 #endif
3093  {
3094  }
3095 
3097  void _move(basic_array&& a) noexcept
3098  {
3099 #ifdef CVM_USE_POOL_MANAGER
3100  mpd = a.mpd;
3101  a.mpd = nullptr;
3102 #else
3103  if (a.mpf == nullptr) {
3104  if (mpf == nullptr) {
3105  mp = std::move(a.mp);
3106  a.msz = 0;
3107  } else {
3108  // *this has foreign array inside, no place to move to
3109  // so, this is the case where we don't even touch a, just copy its content
3110  this->_assign(a.mp.get(), a.incr());
3111  }
3112  } else {
3113  // ... and here we can't rely on foreign poiter life cycle, thus copying
3114  this->_resize(a.size());
3115  this->_assign(a.mpf, a.incr());
3116  //a.mpf = nullptr;
3117  //a.msz = 0;
3118  }
3119 #endif
3120  }
3121 
3122  // compares array elements
3123  bool _equals(const basic_array& a) const
3124  {
3125  bool bRes = false;
3126  if (this->size() == a.size()) {
3127  bRes = true;
3128  if (this->get() != a.get()) {
3129  for (tint i = 0; i < this->size(); ++i) {
3130  if (_abs(this->get()[i * this->incr()] - a.get()[i * a.incr()]) > basic_cvmMachMin<TR>()) {
3131  bRes = false;
3132  break;
3133  }
3134  }
3135  }
3136  }
3137  return bRes;
3138  }
3139 
3140  // array normalizing
3141  void _normalize()
3142  {
3143  const TR dNorm = this->norm();
3144  if (dNorm > basic_cvmMachMin<TR>()) {
3145  static const TR one(1.);
3146  this->_scalr(one / dNorm);
3147  }
3148  }
3149 
3150  // this = a (no matter what)
3151  void _replace(const basic_array& a) throw(cvmexception)
3152  {
3153 #ifdef CVM_USE_POOL_MANAGER
3154  cvmFree<TC>(this->mpd);
3155  this->mpd = cvmMalloc<TC>(a.size());
3156 #else
3157  this->mp.reset(cvmMalloc<TC>(a.size()), ArrayDeleter<TC>());
3158  this->mpf = nullptr;
3159 #endif
3160  this->msz = a.size();
3161  this->mincr = 1;
3162  CVM_ASSERT(this->get(), ((this->size() - 1) * this->incr() + 1) * sizeof(TC))
3163  }
3164 
3165  void _resize(tint nNewSize) throw(cvmexception)
3166  {
3167  _check_lt(CVM_WRONGSIZE_LT, nNewSize, TINT_ZERO);
3168  const bool is_empty = this->_is_empty();
3169  if (nNewSize != this->size() || is_empty) {
3170  TC* pd = cvmMalloc<TC>(nNewSize);
3171  if (nNewSize > this->size()) cvmZeroMemory<TC>(pd, nNewSize);
3172  const tint nMinSize = _cvm_min<tint>(nNewSize, this->size());
3173  if (nMinSize > 0 && !is_empty) {
3174  __copy<TC>(nMinSize, this->get(), this->incr(), pd, 1);
3175  }
3176 #ifdef CVM_USE_POOL_MANAGER
3177  cvmFree<TC>(this->mpd);
3178  this->mpd = pd;
3179 #else
3180  this->mp.reset(pd, ArrayDeleter<TC>());
3181  this->mpf = nullptr;
3182 #endif
3183  this->msz = nNewSize;
3184  this->mincr = 1;
3185  CVM_ASSERT(this->get(), ((this->size() - 1) * this->incr() + 1) * sizeof(TC))
3186  }
3187  }
3188 
3189  bool _is_empty() const
3190  {
3191 #ifdef CVM_USE_POOL_MANAGER
3192  return this->mpd == nullptr;
3193 #else
3194  return this->mpf == nullptr && this->mp.get() == nullptr;
3195 #endif
3196  }
3197 
3198  // virtual methods
3199  virtual void _scalr(TR d) {
3200  __scal<TR,TC>(this->get(), this->size(), this->incr(), d);
3201  }
3202 
3203  // index of max module, undefined for matrices, CVM0-based
3204  virtual tint _indofmax() const {
3205  return __idamax<TC>(this->get(), this->size(), this->incr());
3206  }
3207 
3208  // index of min module, undefined for matrices, CVM0-based
3209  virtual tint _indofmin() const {
3210  return __idamin<TC>(this->get(), this->size(), this->incr());
3211  }
3212 
3213  // fills the content
3214  virtual void _set(TC d)
3215  {
3216  CVM_ASSERT(this->get(), ((this->size() - 1) * this->incr() + 1) * sizeof(TC))
3217  const tint nSize = this->size() * this->incr();
3218  for (tint i = 0; i < nSize; i += this->incr()) {
3219  this->get()[i] = d;
3220  }
3221  }
3222 
3223  virtual void _assign(const TC* pd, tint nIncr)
3224  {
3225  if (this->get() != pd) {
3226  __copy<TC>(this->size(), pd, nIncr, this->get(), this->incr());
3227  }
3228  }
3229 
3230  virtual void _assign_shifted(TC* pDshifted, const TC* pd, tint nSize, tint nIncr, tint)
3231  {
3232  if (pDshifted != pd) {
3233  __copy<TC>(nSize, pd, nIncr, pDshifted, this->incr());
3234  }
3235  }
3237 
3271  friend std::istream& operator >> <> (std::istream& is, basic_array& aIn);
3272 
3296  friend std::ostream& operator << <> (std::ostream& os, const basic_array& aOut);
3297 };
3298 
3299 
3306 template<typename TR>
3307 class basic_rvector : public basic_array<TR,TR>
3308 {
3309  typedef std::complex<TR> TC;
3310  typedef basic_array<TR,TR> BaseArray;
3311 
3312 public:
3338  {
3339  }
3340 
3361  explicit basic_rvector(tint nSize)
3362  : BaseArray(nSize)
3363  {
3364  }
3365 
3366 #if defined(CVM_USE_INITIALIZER_LISTS)
3367 
3388  basic_rvector(std::initializer_list<TR> list)
3389  : BaseArray(static_cast<tint>(list.size()))
3390  {
3391  tint i = 0;
3392  TR* p = this->mp.get();
3393  for (auto it = list.begin(); it != list.end(); ++it) {
3394  p[i++] = *it;
3395  }
3396  }
3397 #endif
3398 
3420  basic_rvector(tint nSize, TR d)
3421  : BaseArray(nSize, false)
3422  {
3423  this->_set(d);
3424  }
3425 
3490  basic_rvector(TR* pd, tint nSize, tint nIncr = 1)
3491  : BaseArray(pd, nSize, nIncr)
3492  {
3493  }
3494 
3536  basic_rvector(const TR* pd, tint nSize, tint nIncr = 1)
3537  : BaseArray(pd, nSize, nIncr)
3538  {
3539  }
3540 
3541 
3571  : BaseArray(v.size(), false)
3572  {
3573  __copy<TR>(this->size(), v, v.incr(), this->get(), this->incr());
3574  }
3575 
3580  : BaseArray(std::move(a))
3581  {
3582  }
3583 
3617  {
3618  _check_ne(CVM_SIZESMISMATCH, v.size(), this->size());
3619  this->_assign(v, v.incr());
3620  return *this;
3621  }
3622 
3627  {
3628  // size check is in BasicArray
3629  BaseArray::operator = (std::move(a));
3630  return *this;
3631  }
3632 
3664  basic_rvector& assign(const TR* pd, tint nIncr = 1)
3665  {
3666  this->_assign(pd, nIncr);
3667  return *this;
3668  }
3669 
3701  basic_rvector& assign(tint n, const TR* pd, tint nIncr = 1)
3702  {
3703  _check_lt_ge(CVM_OUTOFRANGE_LTGE, n, CVM0, this->size() + CVM0);
3704  n -= CVM0;
3705  this->_assign_shifted(this->get() + this->incr() * n, pd, this->size() - n, nIncr, 0);
3706  return *this;
3707  }
3708 
3742  basic_rvector& assign(tint n, const TR* pd, tint nSize, tint nIncr)
3743  {
3744  _check_lt_ge(CVM_OUTOFRANGE_LTGE, n, CVM0, this->size() + CVM0);
3745  n -= CVM0;
3746  this->_assign_shifted(this->get() + this->incr() * n, pd, _cvm_min<tint>(nSize, this->size() - n), nIncr, 0);
3747  return *this;
3748  }
3749 
3781  {
3782  _check_gt(CVM_SIZESMISMATCH_GT, v.size() + n, this->size());
3783  return assign(n, v, v.size(), v.incr());
3784  }
3785 
3809  basic_rvector& set(TR d)
3810  {
3811  this->_set(d);
3812  return *this;
3813  }
3814 
3853  {
3854  this->_resize(nNewSize);
3855  return *this;
3856  }
3857 
3887  bool operator == (const basic_rvector& v) const {
3888  return this->_equals(v);
3889  }
3890 
3920  bool operator != (const basic_rvector& v) const {
3921  return !(this->operator == (v));
3922  }
3923 
3962  {
3963  this->_replace(v);
3964  __copy<TR>(this->size(), v.get(), v.incr(), this->get(), this->incr());
3965  return *this;
3966  }
3967 
4006  {
4007  _check_ne(CVM_SIZESMISMATCH, v.size(), this->size());
4008  basic_rvector vSum(*this);
4009  __add<TR>(vSum.get(), vSum.size(), vSum.incr(), v._pd(), v.incr());
4010  return vSum;
4011  }
4012 
4050  {
4051  _check_ne(CVM_SIZESMISMATCH, v.size(), this->size());
4052  basic_rvector vDiff(*this);
4053  __subtract<TR>(vDiff.get(), vDiff.size(), vDiff.incr(), v._pd(), v.incr());
4054  return vDiff;
4055  }
4056 
4095  {
4096  _check_ne(CVM_SIZESMISMATCH, v1.size(), this->size());
4097  _check_ne(CVM_SIZESMISMATCH, v2.size(), this->size());
4098  _sum<TR,TR>(this->get(), this->size(), this->incr(), v1, v1.incr(), v2, v2.incr());
4099  return *this;
4100  }
4101 
4140  {
4141  _check_ne(CVM_SIZESMISMATCH, v1.size(), this->size());
4142  _check_ne(CVM_SIZESMISMATCH, v2.size(), this->size());
4143  _diff<TR,TR>(this->get(), this->size(), this->incr(), v1, v1.incr(), v2, v2.incr());
4144  return *this;
4145  }
4146 
4187  {
4188  _check_ne(CVM_SIZESMISMATCH, v.size(), this->size());
4189  _incr<TR,TR>(this->get(), this->size(), this->incr(), v, v.incr());
4190  return *this;
4191  }
4192 
4233  {
4234  _check_ne(CVM_SIZESMISMATCH, v.size(), this->size());
4235  _decr<TR,TR>(this->get(), this->size(), this->incr(), v, v.incr());
4236  return *this;
4237  }
4238 
4263  {
4264  static const TR mone(-1.);
4265  basic_rvector vRes(*this);
4266  vRes._scalr(mone);
4267  return vRes;
4268  }
4269 
4295  basic_rvector operator * (TR dMult) const throw(cvmexception)
4296  {
4297  basic_rvector vRes(*this);
4298  vRes._scalr(dMult);
4299  return vRes;
4300  }
4301 
4335  {
4336  basic_rvector vRes(*this);
4337  vRes._div(dDiv);
4338  return vRes;
4339  }
4340 
4366  {
4367  this->_scalr(dMult);
4368  return *this;
4369  }
4370 
4403  {
4404  this->_div(dDiv);
4405  return *this;
4406  }
4407 
4437  {
4438  this->_normalize();
4439  return *this;
4440  }
4441 
4473  TR operator * (const basic_rvector& v) const throw(cvmexception)
4474  {
4475  _check_ne(CVM_SIZESMISMATCH, v.size(), this->size());
4476  return __dot<TR>(this->get(), this->size(), this->incr(), v.get(), v.incr());
4477  }
4478 
4513  {
4514  _check_ne(CVM_SIZESMISMATCH, m.msize(), this->size());
4515  basic_rvector vRes(m.nsize());
4516  m._multiply(vRes, *this, true);
4517  return vRes;
4518  }
4519 
4552  {
4553  _check_ne(CVM_SIZESMISMATCH, m.nsize(), this->size());
4554  _check_ne(CVM_SIZESMISMATCH, m.msize(), v.size());
4555  m._multiply(*this, v, true);
4556  return *this;
4557  }
4558 
4591  {
4592  _check_ne(CVM_SIZESMISMATCH, m.msize(), this->size());
4593  _check_ne(CVM_SIZESMISMATCH, m.nsize(), v.size());
4594  m._multiply(*this, v, false);
4595  return *this;
4596  }
4597 
4642  {
4643  basic_rmatrix<TR> mRes(this->size(), v.size());
4644  static const TR one(1.);
4645  __ger<TR,basic_rmatrix<TR>, basic_rvector>(mRes, *this, v, one);
4646  return mRes;
4647  }
4648 
4695  basic_rvector& solve(const basic_srmatrix<TR>& mA, const basic_rvector& vB, TR& dErr) throw(cvmexception) {
4696  return _solve_helper(mA, vB, dErr, 0);
4697  }
4698 
4699 // 6.1: transpose added
4746  basic_rvector& solve_tran(const basic_srmatrix<TR>& mA, const basic_rvector& vB, TR& dErr) throw(cvmexception) {
4747  return _solve_helper(mA, vB, dErr, 1);
4748  }
4749 
4792  {
4793  static TR dErr(0.);
4794  return this->solve(mA, vB, dErr);
4795  }
4796 
4797 // 6.1: transpose added
4841  {
4842  static TR dErr(0.);
4843  return this->solve_tran(mA, vB, dErr);
4844  }
4845 
4846 // 6.1: MATLAB-style operator B/A returning solution of X*A=B equation
4886  return _operator_solve_helper(mA, 1);
4887  }
4888 
4889  // 6.1: similar to operator / this one returns solution of A*X=B equation
4926  return _operator_solve_helper(mA, 0);
4927  }
4928 
4991  basic_rvector& solve_lu(const basic_srmatrix<TR>& mA, const basic_srmatrix<TR>& mLU, const tint* pPivots,
4992  const basic_rvector& vB, TR& dErr) throw(cvmexception)
4993  {
4994  _check_ne(CVM_SIZESMISMATCH, mA.msize(), this->size());
4995  _check_ne(CVM_SIZESMISMATCH, mA.msize(), vB.size());
4996  _check_ne(CVM_SIZESMISMATCH, mA.msize(), mLU.msize());
4997  mA._solve(vB, *this, dErr, mLU, pPivots, 0);
4998  return *this;
4999  }
5000 
5056  basic_rvector& solve_lu(const basic_srmatrix<TR>& mA, const basic_srmatrix<TR>& mLU, const tint* pPivots,
5057  const basic_rvector& vB) throw(cvmexception)
5058  {
5059  static TR dErr(0.);
5060  return this->solve_lu(mA, mLU, pPivots, vB, dErr);
5061  }
5062 
5117  basic_rvector& gels(bool transpose, const basic_rmatrix<TR>& mA, const basic_rvector& vB,
5118  TR& dErr) throw(cvmexception)
5119  {
5120  _check_ne(CVM_SIZESMISMATCH, vB.size(), transpose ? mA.nsize() : mA.msize());
5121  _check_ne(CVM_SIZESMISMATCH, this->size(), transpose ? mA.msize() : mA.nsize());
5122  basic_rmatrix<TR> mA2(mA); // this algorithm overrides A
5123  basic_rmatrix<TR> mB(vB, vB.size(), 1);
5124  basic_rmatrix<TR> mX;
5125  basic_rvector vErr(1);
5126  __gels(transpose, mA2, mB, mX, vErr);
5127  dErr = vErr(CVM0);
5128  *this = mX(CVM0);
5129  return *this;
5130  }
5131 
5179  TR tol = basic_cvmMachSp<TR>()) throw(cvmexception)
5180  {
5181  _check_ne(CVM_SIZESMISMATCH, vB.size(), mA.msize());
5182  _check_ne(CVM_SIZESMISMATCH, this->size(), mA.nsize());
5183  basic_rmatrix<TR> mA2(mA); // this algorithm overrides A
5184  basic_rmatrix<TR> mB(vB, vB.size(), 1);
5185  basic_rmatrix<TR> mX;
5186  __gelsy(mA2, mB, mX, tol, rank);
5187  *this = mX(CVM0);
5188  return *this;
5189  }
5190 
5244  tint& rank, TR tol = basic_cvmMachSp<TR>()) throw(cvmexception)
5245  {
5246  _gels_sd(true, mA, vB, sv, rank, tol);
5247  return *this;
5248  }
5249 
5304  tint& rank, TR tol = basic_cvmMachSp<TR>()) throw(cvmexception)
5305  {
5306  _gels_sd(false, mA, vB, sv, rank, tol);
5307  return *this;
5308  }
5309 
5383  {
5384  mA._svd(*this, nullptr, nullptr);
5385  return *this;
5386  }
5387 
5461  {
5462  mA._svd(*this, nullptr, nullptr);
5463  return *this;
5464  }
5465 
5544  basic_srmatrix<TR>& mVH) throw(cvmexception)
5545  {
5546  mA._svd(*this, &mU, &mVH);
5547  return *this;
5548  }
5549 
5628  basic_scmatrix<TR,TC>& mVH) throw(cvmexception)
5629  {
5630  mA._svd(*this, &mU, &mVH);
5631  return *this;
5632  }
5633 
5634 // we don't use _eig here since this is the special case - symmetric matrix
5707  {
5708  _check_ne(CVM_SIZESMISMATCH, this->size(), mA.nsize());
5709  __eig<basic_rvector, basic_srsmatrix<TR>, basic_srmatrix<TR> >(*this, mA, nullptr, true);
5710  return *this;
5711  }
5712 
5787  {
5788  _check_ne(CVM_SIZESMISMATCH, this->size(), mA.nsize());
5789  _check_ne(CVM_SIZESMISMATCH, this->size(), mEigVect.nsize());
5790  __eig<basic_rvector, basic_srsmatrix<TR>, basic_srmatrix<TR> >(*this, mA, &mEigVect, true);
5791  return *this;
5792  }
5793 
5794 // we don't use _eig here since this is the special case - hermitian matrix
5867  {
5868  _check_ne(CVM_SIZESMISMATCH, this->size(), mA.nsize());
5869  __eig<basic_rvector, basic_schmatrix<TR,TC>, basic_scmatrix<TR,TC> >(*this, mA, nullptr, true);
5870  return *this;
5871  }
5872 
5947  {
5948  _check_ne(CVM_SIZESMISMATCH, this->size(), mA.nsize());
5949  _check_ne(CVM_SIZESMISMATCH, this->size(), mEigVect.nsize());
5950  __eig<basic_rvector, basic_schmatrix<TR,TC>, basic_scmatrix<TR,TC> >(*this, mA, &mEigVect, true);
5951  return *this;
5952  }
5953 
5954 // ?gemv routines perform a matrix-vector operation defined as
5955 // this = alpha*m*v + beta * this,
6008  basic_rvector& gemv(bool bLeft, const basic_rmatrix<TR>& m, TR dAlpha, const basic_rvector& v, TR dBeta) throw(cvmexception)
6009  {
6010  _check_ne(CVM_SIZESMISMATCH, this->size(), bLeft ? m.nsize() : m.msize());
6011  _check_ne(CVM_SIZESMISMATCH, v.size(), bLeft ? m.msize() : m.nsize());
6012  m._gemv(bLeft, dAlpha, v, dBeta, *this);
6013  return *this;
6014  }
6015 
6068  basic_rvector& gbmv(bool bLeft, const basic_srbmatrix<TR>& m, TR dAlpha, const basic_rvector& v, TR dBeta) throw(cvmexception)
6069  {
6070  _check_ne(CVM_SIZESMISMATCH, this->size(), bLeft ? m.nsize() : m.msize());
6071  _check_ne(CVM_SIZESMISMATCH, v.size(), bLeft ? m.msize() : m.nsize());
6072  m._gbmv(bLeft, dAlpha, v, dBeta, *this);
6073  return *this;
6074  }
6075 
6099  basic_rvector& randomize(TR dFrom, TR dTo)
6100  {
6101  __randomize<TR>(this->get(), this->size(), this->incr(), dFrom, dTo);
6102  return *this;
6103  }
6104 
6105 private:
6107  basic_rvector& _solve_helper(const basic_srmatrix<TR>& mA, const basic_rvector& vB,
6108  TR& dErr, int transp_mode) throw(cvmexception)
6109  {
6110  _check_ne(CVM_SIZESMISMATCH, this->size(), mA.msize());
6111  _check_ne(CVM_SIZESMISMATCH, vB.size(), mA.msize());
6112  mA._solve(vB, *this, dErr, nullptr, nullptr, transp_mode);
6113  return *this;
6114  }
6115 
6116  basic_rvector _operator_solve_helper(const basic_srmatrix<TR>& mA, int transp_mode) const throw(cvmexception)
6117  {
6118  _check_ne(CVM_SIZESMISMATCH, this->size(), mA.msize());
6119  static TR dErr(0.);
6120  basic_rvector vX(this->size());
6121  mA._solve(*this, vX, dErr, nullptr, nullptr, transp_mode);
6122  return vX;
6123  }
6124 
6125  // helper for svd and divide&conquer methods
6126  void _gels_sd(bool svd, const basic_rmatrix<TR>& mA, const basic_rvector& vB,
6127  basic_rvector<TR>& sv, tint& rank, TR tol) throw(cvmexception)
6128  {
6129  _check_ne(CVM_SIZESMISMATCH, this->size(), mA.nsize());
6130  _check_ne(CVM_SIZESMISMATCH, vB.size(), mA.msize());
6131  _check_ne(CVM_SIZESMISMATCH, sv.size(), _cvm_min<tint>(mA.msize(), mA.nsize()));
6132  basic_rmatrix<TR> mA2(mA); // this algorithm overrides A
6133  basic_rmatrix<TR> mB(vB, vB.size(), 1);
6134  basic_rmatrix<TR> mX;
6135  basic_rvector<TR> sv1(sv.size()); // to ensure that incr=1
6136  if (svd) {
6137  __gelss(mA2, mB, mX, tol, sv1, rank);
6138  } else {
6139  __gelsd(mA2, mB, mX, tol, sv1, rank);
6140  }
6141  sv = sv1;
6142  *this = mX(CVM0);
6143  }
6145 };
6146 
6147 
6154 template<typename TR, typename TC>
6155 class basic_cvector : public basic_array<TR,TC>
6156 {
6157  typedef basic_array<TR,TC> BaseArray;
6158 
6159 public:
6185  {
6186  }
6187 
6210  explicit basic_cvector(tint nSize)
6211  : BaseArray(nSize)
6212  {
6213  }
6214 
6215 #if defined(CVM_USE_INITIALIZER_LISTS)
6216 
6242  basic_cvector(std::initializer_list<TC> list)
6243  : BaseArray(static_cast<tint>(list.size()))
6244  {
6245  tint i = 0;
6246  TC* p = this->mp.get();
6247  for (auto it = list.begin(); it != list.end(); ++it) {
6248  p[i++] = *it;
6249  }
6250  }
6251 #endif
6252 
6274  basic_cvector(tint nSize, TC c)
6275  : BaseArray(nSize, false)
6276  {
6277  this->_set(c);
6278  }
6279 
6343  basic_cvector(TC* pd, tint nSize, tint nIncr = 1)
6344  : BaseArray(pd, nSize, nIncr)
6345  {
6346  }
6347 
6388  basic_cvector(const TC* pd, tint nSize, tint nIncr = 1)
6389  : BaseArray(pd, nSize, nIncr)
6390  {
6391  }
6392 
6422  : BaseArray(v.size(), false)
6423  {
6424  __copy<TC>(this->size(), v, v.incr(), this->get(), this->incr());
6425  }
6426 
6445  : BaseArray(std::move(a))
6446  {
6447  }
6448 
6488  basic_cvector(const TR* pRe, const TR* pIm, tint nSize, tint nIncrRe = 1, tint nIncrIm = 1)
6489  : BaseArray(nSize, pRe == nullptr || pIm == nullptr)
6490  {
6491  __copy2<TR,TC>(this->get(), this->size(), this->incr(), pRe, pIm, nIncrRe, nIncrIm);
6492  }
6493 
6528  : BaseArray(vRe.size(), false)
6529  {
6530  _check_ne(CVM_SIZESMISMATCH, vRe.size(), vIm.size());
6531  __copy2<TR,TC>(this->get(), this->size(), this->incr(), vRe, vIm, vRe.incr(), vIm.incr());
6532  }
6533 
6566  basic_cvector(const TR* pA, tint nSize, bool bRealPart = true, tint nIncr = 1)
6567  : BaseArray(nSize)
6568  {
6569  if (bRealPart) __copy2<TR,TC>(this->get(), this->size(), this->incr(), pA, nullptr, nIncr, 0);
6570  else __copy2<TR,TC>(this->get(), this->size(), this->incr(), nullptr, pA, 0, nIncr);
6571  }
6572 
6603  explicit basic_cvector(const basic_rvector<TR>& v, bool bRealPart = true)
6604  : BaseArray(v.size())
6605  {
6606  if (bRealPart) __copy2<TR,TC>(this->get(), this->size(), this->incr(), v, nullptr, v.incr(), 0);
6607  else __copy2<TR,TC>(this->get(), this->size(), this->incr(), nullptr, v, 0, v.incr());
6608  }
6609 
6638  return basic_rvector<TR>(__get_real_p<TR>(this->get()), this->size(), this->incr() * 2);
6639  }
6640 
6669  return basic_rvector<TR>(__get_imag_p<TR>(this->get()), this->size(), this->incr() * 2);
6670  }
6671 
6701  {
6702  _check_ne(CVM_SIZESMISMATCH, this->size(), v.size());
6703  this->_assign(v, v.incr());
6704  return *this;
6705  }
6706 
6721  {
6722  // size check is in BasicArray
6723  BaseArray::operator = (std::move(a));
6724  return *this;
6725  }
6726 
6757  basic_cvector& assign(const TC* pd, tint nIncr = 1)
6758  {
6759  this->_assign(pd, nIncr);
6760  return *this;
6761  }
6762 
6792  basic_cvector& assign(tint n, const TC* pd, tint nIncr = 1)
6793  {
6794  _check_lt_ge(CVM_OUTOFRANGE_LTGE, n, CVM0, this->size() + CVM0);
6795  n -= CVM0;
6796  this->_assign_shifted(this->get() + this->incr() * n, pd, this->size() - n, nIncr, 0);
6797  return *this;
6798  }
6799 
6827  basic_cvector& assign(tint n, const TC* pd, tint nSize, tint nIncr)
6828  {
6829  _check_lt_ge(CVM_OUTOFRANGE_LTGE, n, CVM0, this->size() + CVM0);
6830  n -= CVM0;
6831  this->_assign_shifted(this->get() + this->incr() * n, pd, _cvm_min<tint>(nSize, this->size() - n), nIncr, 0);
6832  return *this;
6833  }
6834 
6863  {
6864  _check_gt(CVM_SIZESMISMATCH_GT, v.size() + n, this->size());
6865  return assign(n, v, v.size(), v.incr());
6866  }
6867 
6891  basic_cvector& set(TC c)
6892  {
6893  this->_set(c);
6894  return *this;
6895  }
6896 
6924  {
6925  _check_ne(CVM_SIZESMISMATCH, vRe.size(), this->size());
6926  __copy_real<TR,TC>(this->get(), this->size(), this->incr(), vRe, vRe.incr());
6927  return *this;
6928  }
6929 
6957  {
6958  _check_ne(CVM_SIZESMISMATCH, vIm.size(), this->size());
6959  __copy_imag<TR,TC>(this->get(), this->size(), this->incr(), vIm, vIm.incr());
6960  return *this;
6961  }
6962 
6986  {
6987  this->_set_real_number(d);
6988  return *this;
6989  }
6990 
7014  {
7015  this->_set_imag_number(d);
7016  return *this;
7017  }
7018 
7055  {
7056  this->_resize(nNewSize);
7057  return *this;
7058  }
7059 
7089  bool operator == (const basic_cvector& v) const {
7090  return this->_equals(v);
7091  }
7092 
7119  bool operator != (const basic_cvector& v) const {
7120  return !(this->operator == (v));
7121  }
7122 
7161  {
7162  this->_replace(v);
7163  __copy<TC>(this->size(), v.get(), v.incr(), this->get(), this->incr());
7164  return *this;
7165  }
7166 
7204  {
7205  _check_ne(CVM_SIZESMISMATCH, v.size(), this->size());
7206  basic_cvector vSum(*this);
7207  __add<TC>(vSum.get(), vSum.size(), vSum.incr(), v._pd(), v.incr());
7208  return vSum;
7209  }
7210 
7247  {
7248  _check_ne(CVM_SIZESMISMATCH, v.size(), this->size());
7249  basic_cvector vDiff(*this);
7250  __subtract<TC>(vDiff.get(), vDiff.size(), vDiff.incr(), v._pd(), v.incr());
7251  return vDiff;
7252  }
7253 
7292  {
7293  _check_ne(CVM_SIZESMISMATCH, v1.size(), this->size());
7294  _check_ne(CVM_SIZESMISMATCH, v2.size(), this->size());
7295  _sum<TR,TC>(this->get(), this->size(), this->incr(), v1, v1.incr(), v2, v2.incr());
7296  return *this;
7297  }
7298 
7337  {
7338  _check_ne(CVM_SIZESMISMATCH, v1.size(), this->size());
7339  _check_ne(CVM_SIZESMISMATCH, v2.size(), this->size());
7340  _diff<TR,TC>(this->get(), this->size(), this->incr(), v1, v1.incr(), v2, v2.incr());
7341  return *this;
7342  }
7343 
7384  {
7385  _check_ne(CVM_SIZESMISMATCH, v.size(), this->size());
7386  _incr<TR,TC>(this->get(), this->size(), this->incr(), v, v.incr());
7387  return *this;
7388  }
7389 
7430  {
7431  _check_ne(CVM_SIZESMISMATCH, v.size(), this->size());
7432  _decr<TR,TC>(this->get(), this->size(), this->incr(), v, v.incr());
7433  return *this;
7434  }
7435 
7460  {
7461  static const TR mone(-1.);
7462  basic_cvector vRes(*this);
7463  vRes._scalr(mone);
7464  return vRes;
7465  }
7466 
7493  basic_cvector operator * (TR dMult) const
7494  {
7495  basic_cvector vRes(*this);
7496  vRes._scalr(dMult);
7497  return vRes;
7498  }
7499 
7529  {
7530  basic_cvector vRes(*this);
7531  vRes._div(dDiv);
7532  return vRes;
7533  }
7534 
7561  basic_cvector operator * (TC cMult) const
7562  {
7563  basic_cvector vRes(*this);
7564  vRes._scalc(cMult);
7565  return vRes;
7566  }
7567 
7597  {
7598  basic_cvector vRes(*this);
7599  vRes._div(cDiv);
7600  return vRes;
7601  }
7602 
7628  {
7629  this->_scalr(dMult);
7630  return *this;
7631  }
7632 
7660  {
7661  this->_div(dDiv);
7662  return *this;
7663  }
7664 
7691  {
7692  this->_scalc(cMult);
7693  return *this;
7694  }
7695 
7724  {
7725  this->_div(cDiv);
7726  return *this;
7727  }
7728 
7756  {
7757  this->_normalize();
7758  return *this;
7759  }
7760 
7790  {
7791  basic_cvector vRes(*this);
7792  __conj<TC>(vRes.get(), vRes.size(), vRes.incr());
7793  return vRes;
7794  }
7795 
7825  {
7826  _check_ne(CVM_SIZESMISMATCH, v.size(), this->size());
7827  this->_assign(v, v.incr());
7828  __conj<TC>(this->get(), this->size(), this->incr());
7829  return *this;
7830  }
7831 
7860  {
7861  __conj<TC>(this->get(), this->size(), this->incr());
7862  return *this;
7863  }
7864 
7890  TC operator * (const basic_cvector& v) const throw(cvmexception)
7891  {
7892  _check_ne(CVM_SIZESMISMATCH, v.size(), this->size());
7893  return __dotu<TC>(this->get(), this->size(), this->incr(), v.get(), v.incr());
7894  }
7895 
7923  TC operator % (const basic_cvector& v) const throw(cvmexception)
7924  {
7925  _check_ne(CVM_SIZESMISMATCH, v.size(), this->size());
7926  return __dotc<TC>(this->get(), this->size(), this->incr(), v.get(), v.incr());
7927  }
7928 
7964  {
7965  _check_ne(CVM_SIZESMISMATCH, m.msize(), this->size());
7966  basic_cvector vRes(m.nsize());
7967  m._multiply(vRes, *this, true);
7968  return vRes;
7969  }
7970 
8013  {
8014  _check_ne(CVM_SIZESMISMATCH, m.nsize(), this->size());
8015  _check_ne(CVM_SIZESMISMATCH, m.msize(), v.size());
8016  m._multiply(*this, v, true);
8017  return *this;
8018  }
8019 
8061  {
8062  _check_ne(CVM_SIZESMISMATCH, m.msize(), this->size());
8063  _check_ne(CVM_SIZESMISMATCH, v.size(), m.nsize());
8064  m._multiply(*this, v, false);
8065  return *this;
8066  }
8067 
8112  {
8113  basic_cmatrix<TR,TC> mRes(this->size(), v.size());
8114  static const TC one(1., 0.);
8115  __geru<TC, basic_cmatrix<TR,TC>, basic_cvector>(mRes, *this, v, one);
8116  return mRes;
8117  }
8118 
8169  {
8170  basic_cmatrix<TR,TC> mRes(this->size(), v.size());
8171  static const TC one(1., 0.);
8172  __gerc<TC, basic_cmatrix<TR,TC>, basic_cvector>(mRes, *this, v, one);
8173  return mRes;
8174  }
8175 
8221  basic_cvector& solve(const basic_scmatrix<TR,TC>& mA, const basic_cvector& vB, TR& dErr) throw(cvmexception) {
8222  return _solve_helper(mA, vB, dErr, 0);
8223  }
8224 
8225 // 6.1: transpose added
8274  return _solve_helper(mA, vB, dErr, 1);
8275  }
8276 
8277 // 6.1: conjugate added
8324  return _solve_helper(mA, vB, dErr, 2);
8325  }
8326 
8367  {
8368  static TR dErr(0.);
8369  return this->solve(mA, vB, dErr);
8370  }
8371 
8372 // 6.1: transpose added
8416  {
8417  static TR dErr(0.);
8418  return this->solve_tran(mA, vB, dErr);
8419  }
8420 
8421 // 6.1: conjugate added
8463  {
8464  static TR dErr(0.);
8465  return this->solve_conj(mA, vB, dErr);
8466  }
8467 
8468 // 6.1: MATLAB-style operator B/A returning solution of X*A=B equation
8508  return _operator_solve_helper(mA, 1);
8509  }
8510 
8511 // 6.1: similar to operator / this one returns solution of A*X=B equation
8548  return _operator_solve_helper(mA, 0);
8549  }
8550 
8615  const basic_cvector& vB, TR& dErr) throw(cvmexception)
8616  {
8617  _check_ne(CVM_SIZESMISMATCH, mA.msize(), this->size());
8618  _check_ne(CVM_SIZESMISMATCH, vB.size(), mA.msize());
8619  _check_ne(CVM_SIZESMISMATCH, mLU.msize(), mA.msize());
8620  mA._solve(vB, *this, dErr, mLU, pPivots, 0);
8621  return *this;
8622  }
8623 
8680  const basic_cvector& vB) throw(cvmexception)
8681  {
8682  static TR dErr(0.);
8683  return this->solve_lu(mA, mLU, pPivots, vB, dErr);
8684  }
8685 
8750  basic_cvector& gels(bool conjugate, const basic_cmatrix<TR,TC>& mA, const basic_cvector& vB,
8751  TC& cErr) throw(cvmexception)
8752  {
8753  _check_ne(CVM_SIZESMISMATCH, vB.size(), conjugate ? mA.nsize() : mA.msize());
8754  _check_ne(CVM_SIZESMISMATCH, this->size(), conjugate ? mA.msize() : mA.nsize());
8755  basic_cmatrix<TR,TC> mA2(mA); // this algorithm overrides A
8756  basic_cmatrix<TR,TC> mB(vB, vB.size(), 1);
8758  basic_cvector vErr(1);
8759  __gels(conjugate, mA2, mB, mX, vErr);
8760  cErr = vErr(CVM0);
8761  *this = mX(CVM0);
8762  return *this;
8763  }
8764 
8814  TR tol = basic_cvmMachSp<TR>()) throw(cvmexception)
8815  {
8816  _check_ne(CVM_SIZESMISMATCH, vB.size(), mA.msize());
8817  _check_ne(CVM_SIZESMISMATCH, this->size(), mA.nsize());
8818  basic_cmatrix<TR,TC> mA2(mA); // this algorithm overrides A
8819  basic_cmatrix<TR,TC> mB(vB, vB.size(), 1);
8821  __gelsy(mA2, mB, mX, tol, rank);
8822  *this = mX(CVM0);
8823  return *this;
8824  }
8825 
8881  tint& rank, TR tol = basic_cvmMachSp<TR>()) throw(cvmexception)
8882  {
8883  _gels_sd(true, mA, vB, sv, rank, tol);
8884  return *this;
8885  }
8886 
8943  tint& rank, TR tol = basic_cvmMachSp<TR>()) throw(cvmexception)
8944  {
8945  _gels_sd(false, mA, vB, sv, rank, tol);
8946  return *this;
8947  }
8948 
8995  {
8996  _check_ne(CVM_SIZESMISMATCH, mA.nsize(), this->size());
8997  mA._eig(*this, nullptr, true);
8998  return *this;
8999  }
9000 
9074  basic_cvector& eig(const basic_srmatrix<TR>& mA, basic_scmatrix<TR,TC>& mEigVect, bool bRightVect = true) throw(cvmexception)
9075  {
9076  _check_ne(CVM_SIZESMISMATCH, mA.nsize(), this->size());
9077  _check_ne(CVM_SIZESMISMATCH, mEigVect.nsize(), this->size());
9078  mA._eig(*this, &mEigVect, bRightVect);
9079  return *this;
9080  }
9081 
9121  {
9122  _check_ne(CVM_SIZESMISMATCH, mA.nsize(), this->size());
9123  mA._eig(*this, nullptr, true);
9124  return *this;
9125  }
9126 
9192  basic_cvector& eig(const basic_scmatrix<TR,TC>& mA, basic_scmatrix<TR,TC>& mEigVect, bool bRightVect = true) throw(cvmexception)
9193  {
9194  _check_ne(CVM_SIZESMISMATCH, mA.nsize(), this->size());
9195  mA._eig(*this, &mEigVect, bRightVect);
9196  return *this;
9197  }
9198 
9245  {
9246  _check_ne(CVM_SIZESMISMATCH, mA.nsize(), this->size());
9247  _check_ne(CVM_SIZESMISMATCH, mB.nsize(), this->size());
9248  _check_ne(CVM_SIZESMISMATCH, vBeta.size(), this->size());
9249  basic_srmatrix<TR> ma(mA); // overriden below. also, we need to care about band ones
9250  basic_srmatrix<TR> mb(mB);
9251  __ggev<basic_srmatrix<TR>, basic_scmatrix<TR,TC>, basic_rvector<TR>, basic_cvector>(ma, mb, *this, vBeta, nullptr, nullptr);
9252  return *this;
9253  }
9254 
9329  basic_scmatrix<TR,TC>& mEigVectLeft, basic_scmatrix<TR,TC>& mEigVectRight) throw(cvmexception)
9330  {
9331  _check_ne(CVM_SIZESMISMATCH, mA.nsize(), this->size());
9332  _check_ne(CVM_SIZESMISMATCH, mB.nsize(), this->size());
9333  _check_ne(CVM_SIZESMISMATCH, vBeta.size(), this->size());
9334  _check_ne(CVM_SIZESMISMATCH, mEigVectLeft.nsize(), this->size());
9335  _check_ne(CVM_SIZESMISMATCH, mEigVectRight.nsize(), this->size());
9336  basic_srmatrix<TR> ma(mA); // overriden below. also, we need to care about band ones
9337  basic_srmatrix<TR> mb(mB);
9338  __ggev(ma, mb, *this, vBeta, &mEigVectLeft, &mEigVectRight);
9339  return *this;
9340  }
9341 
9419  basic_scmatrix<TR,TC>& mEigVect, bool bRightVect = true) throw(cvmexception)
9420  {
9421  _check_ne(CVM_SIZESMISMATCH, mA.nsize(), this->size());
9422  _check_ne(CVM_SIZESMISMATCH, mB.nsize(), this->size());
9423  _check_ne(CVM_SIZESMISMATCH, vBeta.size(), this->size());
9424  _check_ne(CVM_SIZESMISMATCH, mEigVect.nsize(), this->size());
9425  basic_srmatrix<TR> ma(mA); // overriden below. also, we need to care about band ones
9426  basic_srmatrix<TR> mb(mB);
9427  __ggev(ma, mb, *this, vBeta, bRightVect ? nullptr : &mEigVect, bRightVect ? &mEigVect : nullptr);
9428  return *this;
9429  }
9430 
9480  {
9481  _check_ne(CVM_SIZESMISMATCH, mA.nsize(), this->size());
9482  _check_ne(CVM_SIZESMISMATCH, mB.nsize(), this->size());
9483  _check_ne(CVM_SIZESMISMATCH, vBeta.size(), this->size());
9484  basic_scmatrix<TR,TC> ma(mA); // overriden below. also, we need to care about band ones
9485  basic_scmatrix<TR,TC> mb(mB);
9486  __ggev<basic_scmatrix<TR,TC>, basic_scmatrix<TR,TC>, basic_cvector, basic_cvector>(ma, mb, *this, vBeta, nullptr, nullptr);
9487  return *this;
9488  }
9489 
9567  basic_scmatrix<TR,TC>& mEigVectLeft, basic_scmatrix<TR,TC>& mEigVectRight) throw(cvmexception)
9568  {
9569  _check_ne(CVM_SIZESMISMATCH, mA.nsize(), this->size());
9570  _check_ne(CVM_SIZESMISMATCH, mB.nsize(), this->size());
9571  _check_ne(CVM_SIZESMISMATCH, vBeta.size(), this->size());
9572  _check_ne(CVM_SIZESMISMATCH, mEigVectLeft.nsize(), this->size());
9573  _check_ne(CVM_SIZESMISMATCH, mEigVectRight.nsize(), this->size());
9574  basic_scmatrix<TR,TC> ma(mA); // overriden below. also, we need to care about band ones
9575  basic_scmatrix<TR,TC> mb(mB);
9576  __ggev(ma, mb, *this, vBeta, &mEigVectLeft, &mEigVectRight);
9577  return *this;
9578  }
9579 
9660  basic_scmatrix<TR,TC>& mEigVect, bool bRightVect = true) throw(cvmexception)
9661  {
9662  _check_ne(CVM_SIZESMISMATCH, mA.nsize(), this->size());
9663  _check_ne(CVM_SIZESMISMATCH, mB.nsize(), this->size());
9664  _check_ne(CVM_SIZESMISMATCH, vBeta.size(), this->size());
9665  _check_ne(CVM_SIZESMISMATCH, mEigVect.nsize(), this->size());
9666  basic_scmatrix<TR,TC> ma(mA); // overriden below. also, we need to care about band ones
9667  basic_scmatrix<TR,TC> mb(mB);
9668  __ggev(ma, mb, *this, vBeta, bRightVect ? nullptr : &mEigVect, bRightVect ? &mEigVect : nullptr);
9669  return *this;
9670  }
9671 
9725  basic_cvector& gemv(bool bLeft, const basic_cmatrix<TR,TC>& m, TC dAlpha, const basic_cvector& v, TC dBeta) throw(cvmexception)
9726  {
9727  _check_ne(CVM_SIZESMISMATCH, v.size(), bLeft ? m.msize() : m.nsize());
9728  _check_ne(CVM_SIZESMISMATCH, this->size(), bLeft ? m.nsize() : m.msize());
9729  m._gemv(bLeft, dAlpha, v, dBeta, *this);
9730  return *this;
9731  }
9732 
9786  basic_cvector& gbmv(bool bLeft, const basic_scbmatrix<TR,TC>& m, TC dAlpha, const basic_cvector& v, TC dBeta) throw(cvmexception)
9787  {
9788  _check_ne(CVM_SIZESMISMATCH, v.size(), bLeft ? m.msize() : m.nsize());
9789  _check_ne(CVM_SIZESMISMATCH, this->size(), bLeft ? m.nsize() : m.msize());
9790  m._gbmv(bLeft, dAlpha, v, dBeta, *this);
9791  return *this;
9792  }
9793 
9817  basic_cvector& randomize_real(TR dFrom, TR dTo)
9818  {
9819  __randomize_real<TC, TR>(this->get(), this->size(), this->incr(), dFrom, dTo);
9820  return *this;
9821  }
9822 
9846  basic_cvector& randomize_imag(TR dFrom, TR dTo)
9847  {
9848  __randomize_imag<TC, TR>(this->get(), this->size(), this->incr(), dFrom, dTo);
9849  return *this;
9850  }
9851 
9852 protected:
9854  void _div(TC d) throw(cvmexception)
9855  {
9856  if (_abs(d) <= basic_cvmMachMin<TR>()) {
9858  }
9859  static const TC one(1., 0.);
9860  this->_scalc(one / d);
9861  }
9862 
9863  void _scalc(TC d)
9864  {
9865  __scal<TC, TC>(this->get(), this->size(), this->incr(), d);
9866  }
9867 
9868  void _set_real_number(TR d)
9869  {
9870  _set_real<TR,TC>(this->get(), this->size(), this->incr(), d);
9871  }
9872 
9873  void _set_imag_number(TR d)
9874  {
9875  _set_imag<TR,TC>(this->get(), this->size(), this->incr(), d);
9876  }
9877 
9878 private:
9879  basic_cvector& _solve_helper(const basic_scmatrix<TR,TC>& mA, const basic_cvector& vB,
9880  TR& dErr, int transp_mode) throw(cvmexception)
9881  {
9882  _check_ne(CVM_SIZESMISMATCH, mA.msize(), this->size());
9883  _check_ne(CVM_SIZESMISMATCH, mA.msize(), vB.size());
9884  mA._solve(vB, *this, dErr, nullptr, nullptr, transp_mode);
9885  return *this;
9886  }
9887 
9888  basic_cvector _operator_solve_helper(const basic_scmatrix<TR,TC>& mA, int transp_mode) const throw(cvmexception)
9889  {
9890  _check_ne(CVM_SIZESMISMATCH, mA.msize(), this->size());
9891  static TR dErr(0.);
9892  basic_cvector vX(this->size());
9893  mA._solve(*this, vX, dErr, nullptr, nullptr, transp_mode);
9894  return vX;
9895  }
9896 
9897  // helper for svd and divide&conquer methods
9898  void _gels_sd(bool svd, const basic_cmatrix<TR,TC>& mA, const basic_cvector& vB,
9899  basic_rvector<TR>& sv, tint& rank, TR tol = basic_cvmMachSp<TR>()) throw(cvmexception)
9900  {
9901  _check_ne(CVM_SIZESMISMATCH, mA.nsize(), this->size());
9902  _check_ne(CVM_SIZESMISMATCH, mA.msize(), vB.size());
9903  _check_ne(CVM_SIZESMISMATCH, sv.size(), _cvm_min<tint>(mA.msize(), mA.nsize()));
9904  basic_cmatrix<TR,TC> mA2(mA); // this algorithm overrides A
9905  basic_cmatrix<TR,TC> mB(vB, vB.size(), 1);
9907  basic_rvector<TR> sv1(sv.size()); // to ensure that incr=1
9908  if (svd) {
9909  __gelss(mA2, mB, mX, tol, sv1, rank);
9910  } else {
9911  __gelsd(mA2, mB, mX, tol, sv1, rank);
9912  }
9913  sv = sv1;
9914  *this = mX(CVM0);
9915  }
9917 };
9918 
9926 template<typename TR, typename TC>
9927 class Matrix : public basic_array<TR,TC>
9928 {
9929  typedef basic_array<TR,TC> BaseArray;
9930 
9931 protected:
9935 
9942  : mm(0),
9943  mn(0),
9944  mld(0)
9945  {
9946  }
9947 
9957  Matrix(tint nM, tint nN, tint nLD, bool bZeroMemory)
9958  : BaseArray(nLD * nN, bZeroMemory),
9959  mm(nM),
9960  mn(nN),
9961  mld(nLD)
9962  {
9963  }
9964 
9977  Matrix(TC* pd, tint nM, tint nN, tint nLD, tint nSize)
9978  : BaseArray(pd, nSize, 1),
9979  mm(nM),
9980  mn(nN),
9981  mld(nLD)
9982  {
9983  }
9984 
9998  Matrix(const TC* pd, tint nM, tint nN, tint nLD, tint nSize)
9999  : BaseArray(pd, nSize, 1),
10000  mm(nM),
10001  mn(nN),
10002  mld(nLD)
10003  {
10004  }
10005 
10014  Matrix(const BaseArray& v, bool bBeColumn)
10015  : BaseArray(v.size()),
10016  mm(bBeColumn ? v.size() : 1),
10017  mn(bBeColumn ? 1 : v.size()),
10018  mld(mm)
10019  {
10020  __copy<TC>(this->size(), v, v.incr(), this->get(), this->incr());
10021  }
10022 
10040  Matrix(Matrix&& m) noexcept
10041 #if defined(CVM_USE_DELEGATING_CONSTRUCTORS)
10042  : Matrix(m.size(), m.incr(), m.msize(), m.nsize(), m.ld())
10043 #else
10044  : BaseArray(m.size(), m.incr()),
10045  mm(m.msize()),
10046  mn(m.nsize()),
10047  mld(m.ld())
10048 #endif
10049  {
10050  _mmove(std::move(m));
10051  }
10052 
10066  Matrix& operator = (Matrix&& m) throw(cvmexception)
10067  {
10068  _check_ne(CVM_SIZESMISMATCH, m.msize(), this->msize());
10069  _check_ne(CVM_SIZESMISMATCH, m.nsize(), this->nsize());
10070  _mmove(std::move(m));
10071  return *this;
10072  }
10073 
10074 public:
10092  tint msize() const {
10093  return mm;
10094  }
10095 
10113  tint nsize() const {
10114  return mn;
10115  }
10116 
10138  tint ld() const {
10139  return mld;
10140  }
10141 
10165  tint rowofmax() const {
10166  return (this->_indofmax() - CVM0) % mm + CVM0;
10167  }
10168 
10192  tint rowofmin() const {
10193  return (this->_indofmin() - CVM0) % mm + CVM0;
10194  }
10195 
10219  tint colofmax() const {
10220  return (this->_indofmax() - CVM0) / mm + CVM0;
10221  }
10222 
10246  tint colofmin() const {
10247  return (this->_indofmin() - CVM0) / mm + CVM0;
10248  }
10249 
10250  TR norm1() const override
10251  {
10252  tint i, j, k;
10253  TR rSum, rNorm(0.);
10254 
10255  for (j = 0; j < this->nsize(); ++j) {
10256  rSum = TR(0.);
10257 
10258  k = j * this->ld();
10259  for (i = 0; i < this->msize(); ++i) {
10260  CVM_ASSERT(this->get(), (k + i + 1) * sizeof(TC))
10261  rSum += _abs(this->get()[k + i]);
10262  }
10263 
10264  if (rSum > rNorm) {
10265  rNorm = rSum;
10266  }
10267  }
10268  return rNorm;
10269  }
10270 
10271  TR norminf() const override
10272  {
10273  tint i, j;
10274  TR rSum, rNorm(0.);
10275 
10276  for (i = 0; i < this->msize(); ++i) {
10277  rSum = TR(0.);
10278 
10279  for (j = 0; j < this->nsize(); ++j) {
10280  CVM_ASSERT(this->get(), (j * this->ld() + i + 1) * sizeof(TC))
10281  rSum += _abs(this->get()[j * this->ld() + i]);
10282  }
10283 
10284  if (rSum > rNorm) {
10285  rNorm = rSum;
10286  }
10287  }
10288  return rNorm;
10289  }
10290 
10292  const tint* _pm() const {
10293  return &mm;
10294  }
10295 
10296  const tint* _pn() const {
10297  return &mn;
10298  }
10299 
10300  const tint* _pld() const {
10301  return &mld;
10302  }
10303 
10304  TC* _sub_pointer_nocheck(tint row, tint col) {
10305  return this->_pd() + ((col - CVM0) * this->ld() + row - CVM0);
10306  }
10307 
10308  TC* _sub_pointer(tint row, tint col, tint height, tint width) throw(cvmexception)
10309  {
10310  _check_lt(CVM_SIZESMISMATCH_LT, row, CVM0);
10311  _check_lt(CVM_SIZESMISMATCH_LT, col, CVM0);
10312  _check_lt(CVM_SIZESMISMATCH_LT, height, CVM0);
10313  _check_lt(CVM_SIZESMISMATCH_LT, width, CVM0);
10314  _check_gt(CVM_SIZESMISMATCH_GT, row + height, this->msize() + CVM0);
10315  _check_gt(CVM_SIZESMISMATCH_GT, col + width, this->nsize() + CVM0);
10316  return _sub_pointer_nocheck(row, col);
10317  }
10318 
10319  virtual tint _ldm() const {
10320  return this->ld();
10321  }
10322 
10323  virtual const tint* _pldm() const {
10324  return this->_pld();
10325  }
10326 
10327  virtual bool _continuous() const {
10328  return this->msize() == this->ld();
10329  }
10330 
10331  void _check_ld() const
10332  {
10333  if (!this->_continuous()) {
10334  throw cvmexception(CVM_SUBMATRIXACCESSERROR);
10335  }
10336  }
10337 
10338  void _scalr(TR d) override
10339  {
10340  if (this->_continuous()) {
10341  __scal<TR,TC>(this->get(), this->size(), this->incr(), d);
10342  } else for (tint i = 0; i < this->nsize(); ++i) {
10343  __scal<TR,TC>(this->get() + this->ld() * i, this->msize(), this->incr(), d);
10344  }
10345  }
10346 
10347 protected:
10348  Matrix(tint size, tint incr, tint msize, tint nsize, tint ld) noexcept
10349  : BaseArray(size, incr),
10350  mm(msize),
10351  mn(nsize),
10352  mld(ld)
10353  {
10354  }
10355 
10356  // move semantics implementation for matrices
10357  void _mmove(Matrix&& m) noexcept
10358  {
10359  // no size checks here
10360 #ifdef CVM_USE_POOL_MANAGER
10361  mpd = m.mpd;
10362  m.mpd = nullptr;
10363 #else
10364  if (m.mpf == nullptr) {
10365  if (this->mpf == nullptr && this->_continuous() && m._continuous()) {
10366  this->mp = std::move(m.mp);
10367  m.msz = m.mm = m.mn = m.mld = 0;
10368  } else {
10369  // *this has foreign array inside or not continuous; no place to move to
10370  // i.e. this is the case where we don't even touch m, just copy its content
10371  this->_massign(m);
10372  }
10373  } else {
10374  // ... and here we can't rely on foreign poiter life cycle, thus copying
10375  this->_resize2(m.msize(), m.nsize());
10376  this->_massign(m);
10377  }
10378 #endif
10379  }
10380 
10381  void _diag_helper(tint nDiag, tint& nShift, tint& nSize) const throw(cvmexception)
10382  {
10383  if (nDiag >= 0) {
10384  _check_ge(CVM_INDEX_GE, nDiag, this->nsize());
10385  nShift = nDiag * this->ld();
10386  nSize = this->nsize() > this->msize() ? (nDiag > this->nsize() - this->msize() ? this->nsize() - nDiag : this->msize()) :
10387  this->nsize() - nDiag;
10388  } else {
10389  nShift = -nDiag;
10390  _check_ge(CVM_INDEX_GE, nShift, this->msize());
10391  nSize = this->msize() > this->nsize() ? (nShift > this->msize() - this->nsize() ? this->msize() - nShift : this->nsize()) :
10392  this->msize() - nShift;
10393  }
10394  }
10395 
10396  tint _indofmax() const override
10397  {
10398  this->_check_ld();
10399  return __idamax<TC>(this->get(), this->size(), this->incr());
10400  }
10401 
10402  tint _indofmin() const override
10403  {
10404  this->_check_ld();
10405  return __idamin<TC>(this->get(), this->size(), this->incr());
10406  }
10407 
10408  void _assign(const TC* pd, tint nIncr) override
10409  {
10410  if (this->get() != pd) {
10411  if (this->_continuous()) {
10412  __copy<TC>(this->size(), pd, nIncr, this->get(), this->incr());
10413  } else for (tint i = 0; i < this->nsize(); ++i) {
10414  __copy<TC>(this->msize(), pd + this->msize() * i * nIncr, nIncr, this->get() + this->ld() * i, this->incr());
10415  }
10416  }
10417  }
10418 
10419  void _assign_shifted(TC* pDshifted, const TC* pd, tint nRows, tint nCols, tint nLD) override // reusing nSise and nIncr parameter
10420  {
10421  if (pDshifted != pd) {
10422  for (tint i = 0; i < nCols; ++i) {
10423  __copy<TC>(nRows, pd + nLD * i, 1, pDshifted + this->ld() * i, this->incr());
10424  }
10425  }
10426  }
10427 
10428  void _set(TC d) override
10429  {
10430  CVM_ASSERT(this->get(), this->size() * sizeof(TC))
10431  tint i, j, k;
10432  for (j = 0; j < this->nsize(); ++j) {
10433  k = j * this->ld();
10434  for (i = 0; i < this->msize(); ++i) {
10435  this->get()[k + i] = d;
10436  }
10437  }
10438  }
10439 
10440  virtual void _massign(const Matrix& m)
10441  {
10442  if (this->get() != m.get()) {
10443  if (this->_continuous() && m._continuous()) {
10444  __copy<TC>(this->size(), m._pd(), m.incr(), this->get(), this->incr());
10445  } else {
10446  const TC* p = m._pd();
10447  const tint nLD = m._ldm();
10448  for (tint i = 0; i < this->nsize(); ++i) {
10449  __copy<TC>(this->msize(), p + nLD * i, m.incr(), this->get() + this->ld() * i, this->incr());
10450  }
10451  }
10452  }
10453  }
10454 
10455  virtual void _resize2(tint nNewM, tint nNewN) throw(cvmexception)
10456  {
10457  _check_lt(CVM_WRONGSIZE_LT, nNewM, TINT_ZERO);
10458  _check_lt(CVM_WRONGSIZE_LT, nNewN, TINT_ZERO);
10459  const bool is_empty = this->_is_empty();
10460  if (nNewM != this->msize() || nNewN != this->nsize() || is_empty) {
10461  if (!is_empty) {
10462  this->_check_ld();
10463  }
10464  const tint nNewSize = nNewM * nNewN;
10465  TC* pd = cvmMalloc<TC>(nNewSize);
10466  cvmZeroMemory<TC>(pd, nNewSize);
10467  const tint nMinM = _cvm_min<tint>(nNewM, this->msize());
10468  const tint nMinN = _cvm_min<tint>(nNewN, this->nsize());
10469  if (nNewSize > 0 && !is_empty) {
10470  for (tint i = 0; i < nMinN; ++i) {
10471  __copy<TC>(nMinM, this->get() + i * this->msize(), this->incr(), pd + i * nNewM, 1);
10472  }
10473  }
10474 #ifdef CVM_USE_POOL_MANAGER
10475  cvmFree<TC>(this->mpd);
10476  this->mpd = pd;
10477 #else
10478  this->mp.reset(pd, ArrayDeleter<TC>());
10479  this->mpf = nullptr;
10480 #endif
10481  this->msz = nNewSize;
10482  CVM_ASSERT(this->get(), this->size() * sizeof(TC))
10483  this->mm = nNewM;
10484  this->mn = nNewN;
10485  this->mld = nNewM;
10486  this->mincr = 1;
10487  }
10488  }
10489 
10490  virtual tint _ld_for_replace() const {
10491  return this->mld;
10492  }
10493 
10494  virtual tint _size_for_replace() const {
10495  return this->size();
10496  }
10497 
10498  void _replace(const Matrix& m) throw(cvmexception)
10499  {
10500  // submatrix replacement is obviously not possible
10501  this->_check_ld();
10502 #ifdef CVM_USE_POOL_MANAGER
10503  cvmFree<TC>(this->mpd);
10504  this->msz = m._size_for_replace();
10505  this->mpd = cvmMalloc<TC>(this->size());
10506 #else
10507  this->msz = m._size_for_replace();
10508  this->mp.reset(cvmMalloc<TC>(this->size()), ArrayDeleter<TC>());
10509  this->mpf = nullptr;
10510 #endif
10511  this->mincr = 1;
10512  CVM_ASSERT(this->get(), (this->size() * sizeof(TC)))
10513  this->mm = m.msize();
10514  this->mn = m.nsize();
10515  this->mld = m._ld_for_replace();
10516  }
10517 
10518  void _transp_m(const Matrix& m)
10519  {
10520  tint i;
10521  if (this->msize() > this->nsize()) for (i = 0; i < this->nsize(); ++i) {
10522  __copy<TC>(m.nsize(), m.get() + i, m.ld(), this->get() + i * this->ld(), 1);
10523  }
10524  else for (i = 0; i < this->msize(); ++i) {
10525  __copy<TC>(m.msize(), m.get() + i * m.ld(), 1, this->get() + i, this->ld());
10526  }
10527  }
10528 
10529  virtual type_proxy<TC, TR> _ij_proxy_val(tint i, tint j) { // always zero based
10530  CVM_ASSERT(this->get(), (this->ld() * j + i + 1) * sizeof(TC))
10531  return type_proxy<TC, TR>(this->get()[this->ld() * j + i], false);
10532  }
10533 
10534  virtual TC _ij_val(tint i, tint j) const { // always zero based
10535  CVM_ASSERT(this->get(), (this->ld() * j + 1) * sizeof(TC))
10536  return this->get()[this->ld() * j + i];
10537  }
10538 
10539  virtual void _swap_rows(tint n1, tint n2) throw(cvmexception)
10540  {
10541  _check_lt_ge(CVM_OUTOFRANGE_LTGE1, n1, CVM0, this->msize() + CVM0);
10542  _check_lt_ge(CVM_OUTOFRANGE_LTGE2, n2, CVM0, this->msize() + CVM0);
10543  if (n1 != n2) {
10544  __swap<TC>(this->nsize(), this->get() + n1 - CVM0, this->ld(), this->get() + n2 - CVM0, this->ld());
10545  }
10546  }
10547 
10548  virtual void _swap_cols(tint n1, tint n2) throw(cvmexception)
10549  {
10550  _check_lt_ge(CVM_OUTOFRANGE_LTGE1, n1, CVM0, this->nsize() + CVM0);
10551  _check_lt_ge(CVM_OUTOFRANGE_LTGE1, n2, CVM0, this->nsize() + CVM0);
10552  if (n1 != n2) {
10553  __swap<TC>(this->msize(), this->get() + (n1 - CVM0) * this->ld(), 1, this->get() + (n2 - CVM0) * this->ld(), 1);
10554  }
10555  }
10556 
10557  virtual const TC* _pp(const Matrix& m) const {
10558  return m._pd();
10559  }
10560 
10561  // matrix cleaning (we ALWAYS have mincr = 1 for matrices)
10562  virtual void _vanish()
10563  {
10564  CVM_ASSERT(this->get(), this->size() * sizeof(TC))
10565  if (this->_continuous()) {
10566  memset(this->get(), 0, this->size() * sizeof(TC));
10567  } else for (tint i = 0; i < this->nsize(); ++i) {
10568  memset(this->get() + this->ld() * i, 0, this->msize() * sizeof(TC));
10569  }
10570  }
10571 
10572  void _msum(const Matrix& m1, const Matrix& m2)
10573  {
10574  if (this->_continuous() && m1._continuous() && m2._continuous()) {
10575  _sum<TR,TC>(this->get(), this->size(), this->incr(), _pp(m1), m1.incr(), _pp(m2), m2.incr());
10576  } else for (tint i = 0; i < this->nsize(); ++i) {
10577  _sum<TR,TC>(this->get() + this->ld() * i, this->msize(), this->incr(), _pp(m1) + m1._ldm() * i, m1.incr(), _pp(m2) + m2._ldm() * i, m2.incr());
10578  }
10579  }
10580 
10581  void _mdiff(const Matrix& m1, const Matrix& m2)
10582  {
10583  if (this->_continuous() && m1._continuous() && m2._continuous()) {
10584  _diff<TR,TC>(this->get(), this->size(), this->incr(), _pp(m1), m1.incr(), _pp(m2), m2.incr());
10585  } else for (tint i = 0; i < this->nsize(); ++i) {
10586  _diff<TR,TC>(this->get() + this->ld() * i, this->msize(), this->incr(), _pp(m1) + m1._ldm() * i, m1.incr(), _pp(m2) + m2._ldm() * i, m2.incr());
10587  }
10588  }
10589 
10590  void _mincr(const Matrix& m)
10591  {
10592  if (this->_continuous() && m._continuous()) {
10593  _incr<TR,TC>(this->get(), this->size(), this->incr(), _pp(m), m.incr());
10594  } else for (tint i = 0; i < this->nsize(); ++i) {
10595  _incr<TR,TC>(this->get() + this->ld() * i, this->msize(), this->incr(), _pp(m) + m._ldm() * i, m.incr());
10596  }
10597  }
10598 
10599  void _mdecr(const Matrix& m)
10600  {
10601  if (this->_continuous() && m._continuous()) {
10602  _decr<TR,TC>(this->get(), this->size(), this->incr(), _pp(m), m.incr());
10603  } else for (tint i = 0; i < this->nsize(); ++i) {
10604  _decr<TR,TC>(this->get() + this->ld() * i, this->msize(), this->incr(), _pp(m) + m._ldm() * i, m.incr());
10605  }
10606  }
10608 
10609  public:
10610  friend std::ostream& operator << <> (std::ostream& os, const Matrix<TR,TC>& mOut);
10611  friend std::istream& operator >> <> (std::istream& is, Matrix<TR,TC>& mIn);
10612 };
10613 
10614 
10622 template<typename TR, typename TC>
10624 {
10625  typedef Matrix<TR,TC> BaseMatrix;
10626 
10627 protected:
10630  {
10631  }
10632 
10633  virtual ~SqMatrix()
10634  {
10635  }
10636 
10638  virtual tint _size() const = 0;
10639  virtual tint _msize() const = 0;
10640  virtual tint _nsize() const = 0;
10641  virtual tint _ld() const = 0;
10642 
10643  // it's the same as get, but let's keep get not virtual for performance sake
10644  virtual const TC* _pv() const = 0;
10645  virtual TC* _pv() = 0;
10646 
10647  // it differs from Matrix::_transp_m because in this case we can do it in-place.
10648  void _sq_transp()
10649  {
10650  const tint mm = this->_msize();
10651  const tint mld = this->_ld();
10652  TC* pd = this->_pv();
10653  if (mm > 1) {
10654  const tint nM1 = mld + 1, nM1m = mld - 1, nM2m = mm - 1;
10655  tint i = 1, j = 1, m;
10656  for (;;) {
10657  m = mm - i;
10658  __swap<TC>(m, pd + j, 1, pd + j + nM1m, mld);
10659  if (i >= nM2m) {
10660  break;
10661  }
10662  ++i;
10663  j += nM1;
10664  }
10665  }
10666  }
10667 
10668  // plus identity
10669  void _sq_plus_plus()
10670  {
10671  TC* pd = this->_pv();
10672  static const TC one(1.);
10673  const tint nSize = this->_size();
10674  const tint nNext = this->_msize() + 1;
10675  CVM_ASSERT(pd, nSize * sizeof(TC))
10676  for (tint i = 0; i < nSize; i += nNext) {
10677  pd[i] += one;
10678  }
10679  }
10680 
10681  // minus identity
10682  void _sq_minus_minus()
10683  {
10684  TC* pd = this->_pv();
10685  static const TC one(1.);
10686  const tint nSize = this->_size();
10687  const tint nNext = this->_msize() + 1;
10688  CVM_ASSERT(pd, nSize * sizeof(TC))
10689  for (tint i = 0; i < nSize; i += nNext) {
10690  pd[i] -= one;
10691  }
10692  }
10693 
10694 public:
10695  void _clean_low_triangle()
10696  {
10697  const tint mm = this->_msize();
10698  const tint mld = this->_ld();
10699  TC* pd = this->_pv();
10700  tint n = 1;
10701  static const TR zero(0.);
10702  for (tint i = 1; i < mm; ++i) {
10703  __scal<TR,TC>(pd + n, mm - i, 1, zero); // column by column
10704  n += mld + 1;
10705  }
10706  }
10708 };
10709 
10716 template<typename TR>
10717 class basic_rmatrix : public Matrix<TR,TR>
10718 {
10719  typedef std::complex<TR> TC;
10720  typedef basic_array<TR,TR> BaseArray;
10721  typedef Matrix<TR,TR> BaseMatrix;
10722  typedef basic_rvector<TR> RVector;
10723 
10724  friend class basic_rvector<TR>; // _multiply
10725 
10726 public:
10750  {
10751  }
10752 
10780  : BaseMatrix(nM, nN, nM, true)
10781  {
10782  }
10783 
10816  basic_rmatrix(TR* pd, tint nM, tint nN)
10817  : BaseMatrix(pd, nM, nN, nM, nM * nN)
10818  {
10819  }
10820 
10855  basic_rmatrix(const TR* pd, tint nM, tint nN)
10856  : BaseMatrix(pd, nM, nN, nM, nM * nN)
10857  {
10858  }
10859 
10888  : BaseMatrix(m.msize(), m.nsize(), m.msize(), false)
10889  {
10890  this->_massign(m);
10891  }
10892 
10911  : BaseMatrix(std::move(m))
10912  {
10913  }
10914 
10949  explicit basic_rmatrix(const RVector& v, bool bBeColumn = true)
10950  : BaseMatrix(v, bBeColumn)
10951  {
10952  }
10953 
10984  basic_rmatrix(basic_rmatrix& m, tint nRow, tint nCol, tint nHeight, tint nWidth)
10985  : BaseMatrix(m._sub_pointer(nRow, nCol, nHeight, nWidth), nHeight, nWidth, m.ld(), nHeight * nWidth)
10986  {
10987  m._check_submatrix();
10988  }
10989 
11028  type_proxy<TR,TR> operator () (tint nRow, tint nCol) throw(cvmexception)
11029  {
11030  _check_lt_ge(CVM_OUTOFRANGE_LTGE1, nRow, CVM0, this->msize() + CVM0);
11031  _check_lt_ge(CVM_OUTOFRANGE_LTGE2, nCol, CVM0, this->nsize() + CVM0);
11032  return this->_ij_proxy_val(nRow - CVM0, nCol - CVM0);
11033  }
11034 
11065  TR operator () (tint nRow, tint nCol) const throw(cvmexception)
11066  {
11067  _check_lt_ge(CVM_OUTOFRANGE_LTGE1, nRow, CVM0, this->msize() + CVM0);
11068  _check_lt_ge(CVM_OUTOFRANGE_LTGE2, nCol, CVM0, this->nsize() + CVM0);
11069  return this->_ij_val(nRow - CVM0, nCol - CVM0);
11070  }
11071 
11107  RVector operator () (tint nCol) throw(cvmexception)
11108  {
11109  _check_lt_ge(CVM_OUTOFRANGE_LTGE, nCol, CVM0, this->nsize() + CVM0);
11110  return this->_col(nCol - CVM0);
11111  }
11112 
11149  RVector operator [] (tint nRow) throw(cvmexception)
11150  {
11151  _check_lt_ge(CVM_OUTOFRANGE_LTGE, nRow, CVM0, this->msize() + CVM0);
11152  return this->_row(nRow - CVM0);
11153  }
11154 
11181  const RVector operator () (tint nCol) const throw(cvmexception)
11182  {
11183  _check_lt_ge(CVM_OUTOFRANGE_LTGE, nCol, CVM0, this->nsize() + CVM0);
11184  return this->_col(nCol - CVM0);
11185  }
11186 
11213  const RVector operator [] (tint nRow) const throw(cvmexception)
11214  {
11215  _check_lt_ge(CVM_OUTOFRANGE_LTGE, nRow, CVM0, this->msize() + CVM0);
11216  return this->_row(nRow - CVM0);
11217  }
11218 
11252  RVector diag(tint nDiag) throw(cvmexception) {
11253  return this->_diag(nDiag);
11254  }
11255 
11288  const RVector diag(tint nDiag) const throw(cvmexception) {
11289  return this->_diag(nDiag);
11290  }
11291 
11325  basic_rmatrix& operator = (const basic_rmatrix& m) throw(cvmexception)
11326  {
11327  _check_ne(CVM_SIZESMISMATCH, this->msize(), m.msize());
11328  _check_ne(CVM_SIZESMISMATCH, this->nsize(), m.nsize());
11329  this->_massign(m);
11330  return *this;
11331  }
11332 
11346  basic_rmatrix& operator = (basic_rmatrix&& m) throw(cvmexception)
11347  {
11348  // size check is in BaseMatrix
11349  BaseMatrix::operator = (std::move(m));
11350  return *this;
11351  }
11352 
11382  basic_rmatrix& assign(const RVector& v) throw(cvmexception)
11383  {
11384  _check_gt(CVM_SIZESMISMATCH_GT, this->size(), v.size());
11385  this->_assign(v, v.incr());
11386  return *this;
11387  }
11388 
11416  basic_rmatrix& assign(const TR* pd)
11417  {
11418  this->_assign(pd, 1);
11419  return *this;
11420  }
11421 
11452  basic_rmatrix& assign(tint nRow, tint nCol, const basic_rmatrix& m) throw(cvmexception)
11453  {
11454  _check_lt_ge(CVM_OUTOFRANGE_LTGE1, nRow, CVM0, this->msize() + CVM0);
11455  _check_lt_ge(CVM_OUTOFRANGE_LTGE2, nCol, CVM0, this->nsize() + CVM0);
11456  _check_gt(CVM_SIZESMISMATCH_GT, m.msize() + nRow - CVM0, this->msize());
11457  _check_gt(CVM_SIZESMISMATCH_GT, m.nsize() + nCol - CVM0, this->nsize());
11458  this->_assign_shifted(this->_sub_pointer_nocheck(nRow, nCol), m._pd(), m.msize(), m.nsize(), m.ld());
11459  return *this;
11460  }
11461 
11484  basic_rmatrix& set(TR d)
11485  {
11486  this->_set(d);
11487  return *this;
11488  }
11489 
11533  basic_rmatrix& resize(tint nNewM, tint nNewN) throw(cvmexception)
11534  {
11535  this->_resize2(nNewM, nNewN);
11536  return *this;
11537  }
11538 
11565  bool operator == (const basic_rmatrix& m) const {
11566  return this->msize() == m.msize() && this->nsize() == m.nsize() && this->_mequals(m);
11567  }
11568 
11595  bool operator != (const basic_rmatrix& m) const {
11596  return !operator == (m);
11597  }
11598 
11639  basic_rmatrix& operator << (const basic_rmatrix& m) throw(cvmexception)
11640  {
11641  this->_replace(m);
11642  this->_massign(m);
11643  return *this;
11644  }
11645 
11682  basic_rmatrix operator + (const basic_rmatrix& m) const throw(cvmexception)
11683  {
11684  _check_ne(CVM_SIZESMISMATCH, this->msize(), m.msize());
11685  _check_ne(CVM_SIZESMISMATCH, this->nsize(), m.nsize());
11686  basic_rmatrix mSum(*this);
11687  mSum._mincr(m);
11688  return mSum;
11689  }
11690 
11726  basic_rmatrix operator - (const basic_rmatrix& m) const throw(cvmexception)
11727  {
11728  _check_ne(CVM_SIZESMISMATCH, this->msize(), m.msize());
11729  _check_ne(CVM_SIZESMISMATCH, this->nsize(), m.nsize());
11730  basic_rmatrix mDiff(*this);
11731  mDiff._mdecr(m);
11732  return mDiff;
11733  }
11734 
11771  basic_rmatrix& sum(const basic_rmatrix& m1, const basic_rmatrix& m2) throw(cvmexception)
11772  {
11773  _check_ne(CVM_SIZESMISMATCH, this->msize(), m1.msize());
11774  _check_ne(CVM_SIZESMISMATCH, this->nsize(), m1.nsize());
11775  _check_ne(CVM_SIZESMISMATCH, this->msize(), m2.msize());
11776  _check_ne(CVM_SIZESMISMATCH, this->nsize(), m2.nsize());
11777  this->_msum(m1, m2);
11778  return *this;
11779  }
11780 
11817  basic_rmatrix& diff(const basic_rmatrix& m1, const basic_rmatrix& m2) throw(cvmexception)
11818  {
11819  _check_ne(CVM_SIZESMISMATCH, this->msize(), m1.msize());
11820  _check_ne(CVM_SIZESMISMATCH, this->nsize(), m1.nsize());
11821  _check_ne(CVM_SIZESMISMATCH, this->msize(), m2.msize());
11822  _check_ne(CVM_SIZESMISMATCH, this->nsize(), m2.nsize());
11823  this->_mdiff(m1, m2);
11824  return *this;
11825  }
11826 
11866  basic_rmatrix& operator += (const basic_rmatrix& m) throw(cvmexception)
11867  {
11868  _check_ne(CVM_SIZESMISMATCH, this->msize(), m.msize());
11869  _check_ne(CVM_SIZESMISMATCH, this->nsize(), m.nsize());
11870  this->_mincr(m);
11871  return *this;
11872  }
11873 
11913  basic_rmatrix& operator -= (const basic_rmatrix& m) throw(cvmexception)
11914  {
11915  _check_ne(CVM_SIZESMISMATCH, this->msize(), m.msize());
11916  _check_ne(CVM_SIZESMISMATCH, this->nsize(), m.nsize());
11917  this->_mdecr(m);
11918  return *this;
11919  }
11920 
11942  static const TR mone(-1.);
11943  basic_rmatrix mRes(*this);
11944  mRes._scalr(mone);
11945  return mRes;
11946  }
11947 
11970  basic_rmatrix operator * (TR dMult) const {
11971  basic_rmatrix mRes(*this);
11972  mRes._scalr(dMult);
11973  return mRes;
11974  }
11975 
12008  basic_rmatrix operator / (TR dDiv) const throw(cvmexception)
12009  {
12010  basic_rmatrix mRes(*this);
12011  mRes._div(dDiv);
12012  return mRes;
12013  }
12014 
12038  {
12039  this->_scalr(dMult);
12040  return *this;
12041  }
12042 
12075  basic_rmatrix& operator /= (TR dDiv) throw(cvmexception)
12076  {
12077  this->_div(dDiv);
12078  return *this;
12079  }
12080 
12107  {
12108  this->_normalize();
12109  return *this;
12110  }
12111 
12139  basic_rmatrix operator ~() const throw(cvmexception)
12140  {
12141  basic_rmatrix mRes(this->nsize(), this->msize());
12142  mRes._transp_m(*this);
12143  return mRes;
12144  }
12145 
12178  basic_rmatrix& transpose(const basic_rmatrix& m) throw(cvmexception)
12179  {
12180  _check_ne(CVM_SIZESMISMATCH, this->msize(), m.nsize());
12181  _check_ne(CVM_SIZESMISMATCH, this->nsize(), m.msize());
12182  if (this->get() == m.get()) {
12183  basic_rmatrix mTmp(m);
12184  this->_transp_m(mTmp);
12185  } else {
12186  this->_transp_m(m);
12187  }
12188  return *this;
12189  }
12190 
12220  basic_rmatrix& transpose() throw(cvmexception)
12221  {
12222  basic_rmatrix mTmp(*this);
12223  this->_resize2(this->nsize(), this->msize());
12224  this->_transp_m(mTmp);
12225  return *this;
12226  }
12227 
12257  RVector operator * (const RVector& v) const throw(cvmexception)
12258  {
12259  _check_ne(CVM_SIZESMISMATCH, this->nsize(), v.size());
12260  RVector vRes(this->msize());
12261  this->_multiply(vRes, v, false);
12262  return vRes;
12263  }
12264 
12294  basic_rmatrix operator * (const basic_rmatrix& m) const throw(cvmexception)
12295  {
12296  _check_ne(CVM_SIZESMISMATCH, this->nsize(), m.msize());
12297  basic_rmatrix mRes(this->msize(), m.nsize());
12298  mRes.mult(*this, m);
12299  return mRes;
12300  }
12301 
12336  basic_rmatrix& mult(const basic_rmatrix& m1, const basic_rmatrix& m2) throw(cvmexception)
12337  {
12338  this->_mult(m1, m2);
12339  return *this;
12340  }
12341 
12394  basic_rmatrix& rank1update(const RVector& vCol, const RVector& vRow) throw(cvmexception)
12395  {
12396  static const TR one(1.);
12397  this->_check_rank1update();
12398  _check_ne(CVM_SIZESMISMATCH, this->msize(), vCol.size());
12399  _check_ne(CVM_SIZESMISMATCH, this->nsize(), vRow.size());
12400  this->_vanish();
12401  __ger<TR,basic_rmatrix, RVector>(*this, vCol, vRow, one);
12402  return *this;
12403  }
12404 
12438  basic_rmatrix& swap_rows(tint n1, tint n2) throw(cvmexception)
12439  {
12440  this->_swap_rows(n1, n2);
12441  return *this;
12442  }
12443 
12475  basic_rmatrix& swap_cols(tint n1, tint n2) throw(cvmexception)
12476  {
12477  this->_swap_cols(n1, n2);
12478  return *this;
12479  }
12480 
12522  basic_rmatrix& solve(const basic_srmatrix<TR>& mA, const basic_rmatrix& mB, TR& dErr) throw(cvmexception) {
12523  return _solve_helper(mA, mB, dErr, 0);
12524  }
12525 
12526 // 6.1: transpose added
12574  basic_rmatrix& solve_tran(const basic_srmatrix<TR>& mA, const basic_rmatrix& mB, TR& dErr) throw(cvmexception) {
12575  return _solve_helper(mA, mB, dErr, 1);
12576  }
12577 
12615  basic_rmatrix& solve(const basic_srmatrix<TR>& mA, const basic_rmatrix& mB) throw(cvmexception)
12616  {
12617  static TR dErr(0.);
12618  return this->solve(mA, mB, dErr);
12619  }
12620 
12621 // 6.1: transpose added
12666  basic_rmatrix& solve_tran(const basic_srmatrix<TR>& mA, const basic_rmatrix& mB) throw(cvmexception)
12667  {
12668  static TR dErr(0.);
12669  return this->solve_tran(mA, mB, dErr);
12670  }
12671 
12741  basic_rmatrix& solve_lu(const basic_srmatrix<TR>& mA, const basic_srmatrix<TR>& mLU, const tint* pPivots,
12742  const basic_rmatrix& mB, TR& dErr) throw(cvmexception)
12743  {
12744  _check_ne(CVM_SIZESMISMATCH, this->msize(), mA.msize());
12745  _check_ne(CVM_SIZESMISMATCH, this->nsize(), mB.nsize());
12746  _check_ne(CVM_SIZESMISMATCH, mA.msize(), mB.msize());
12747  _check_ne(CVM_SIZESMISMATCH, mA.msize(), mLU.msize());
12748  mA._solve(mB, *this, dErr, mLU, pPivots, 0);
12749  return *this;
12750  }
12751 
12818  basic_rmatrix& solve_lu(const basic_srmatrix<TR>& mA, const basic_srmatrix<TR>& mLU, const tint* pPivots,
12819  const basic_rmatrix& mB) throw(cvmexception)
12820  {
12821  static TR dErr(0.);
12822  return this->solve_lu(mA, mLU, pPivots, mB, dErr);
12823  }
12824 
12867  RVector svd() const throw(cvmexception)
12868  {
12869  RVector vRes(_cvm_min<tint>(this->msize(), this->nsize()));
12870  this->_svd(vRes, nullptr, nullptr);
12871  return vRes;
12872  }
12873 
12949  RVector svd(basic_srmatrix<TR>& mU, basic_srmatrix<TR>& mVH) const throw(cvmexception)
12950  {
12951  RVector vRes(_cvm_min<tint>(this->msize(), this->nsize()));
12952  this->_svd(vRes, &mU, &mVH);
12953  return vRes;
12954  }
12955 
13009  basic_rmatrix pinv(TR threshold = basic_cvmMachSp<TR>()) const throw(cvmexception)
13010  {
13011  basic_rmatrix mAx(this->nsize(), this->msize());
13012  this->_pinv(mAx, threshold);
13013  return mAx;
13014  }
13015 
13073  basic_rmatrix& pinv(const basic_rmatrix& mA, TR threshold = basic_cvmMachSp<TR>()) throw(cvmexception)
13074  {
13075  _check_ne(CVM_SIZESMISMATCH, this->nsize(), mA.msize());
13076  _check_ne(CVM_SIZESMISMATCH, this->msize(), mA.nsize());
13077  mA._pinv(*this, threshold);
13078  return *this;
13079  }
13080 
13146  basic_rmatrix gels(bool transpose, const basic_rmatrix& mB, basic_rvector<TR>& vErr) const throw(cvmexception)
13147  {
13148  _check_ne(CVM_SIZESMISMATCH, mB.nsize(), vErr.size());
13149  _check_ne(CVM_SIZESMISMATCH, mB.msize(), transpose ? this->nsize() : this->msize());
13150  basic_rmatrix mA(*this); // this algorithm overrides A
13151  basic_rmatrix mX;
13152  __gels(transpose, mA, mB, mX, vErr);
13153  return mX;
13154  }
13155 
13222  basic_rmatrix& gels(bool transpose, const basic_rmatrix& mA, const basic_rmatrix& mB,
13223  basic_rvector<TR>& vErr) throw(cvmexception)
13224  {
13225  _check_ne(CVM_SIZESMISMATCH, this->nsize(), mB.nsize());
13226  _check_ne(CVM_SIZESMISMATCH, mB.nsize(), vErr.size());
13227  _check_ne(CVM_SIZESMISMATCH, this->msize(), transpose ? mA.msize() : mA.nsize());
13228  _check_ne(CVM_SIZESMISMATCH, mB.msize(), transpose ? mA.nsize() : mA.msize());
13229  basic_rmatrix mA2(mA); // this algorithm overrides A
13230  __gels(transpose, mA2, mB, *this, vErr);
13231  return *this;
13232  }
13233 
13295  basic_rvector<TR> gels(bool transpose, const basic_rvector<TR>& vB, TR& dErr) const throw(cvmexception)
13296  {
13297  _check_ne(CVM_SIZESMISMATCH, vB.size(), transpose ? this->nsize() : this->msize());
13298  basic_rmatrix mA(*this); // this algorithm overrides A
13299  basic_rmatrix mB(vB, vB.size(), 1);
13300  basic_rmatrix mX;
13301  basic_rvector<TR> vErr(1);
13302  __gels(transpose, mA, mB, mX, vErr);
13303  dErr = vErr(CVM0);
13304  const TR* pResult = mX;
13305  return basic_rvector<TR>(pResult, mX.msize());
13306  }
13307 
13361  TR tol = basic_cvmMachSp<TR>()) const throw(cvmexception)
13362  {
13363  _check_ne(CVM_SIZESMISMATCH, this->msize(), mB.msize());
13364  basic_rmatrix mA(*this); // this algorithm overrides A
13365  basic_rmatrix mX;
13366  __gelsy(mA, mB, mX, tol, rank);
13367  return mX;
13368  }
13369 
13425  tint& rank, TR tol = basic_cvmMachSp<TR>()) throw(cvmexception)
13426  {
13427  _check_ne(CVM_SIZESMISMATCH, this->msize(), mA.nsize());
13428  _check_ne(CVM_SIZESMISMATCH, this->nsize(), mB.nsize());
13429  _check_ne(CVM_SIZESMISMATCH, mA.msize(), mB.msize());
13430  basic_rmatrix mA2(mA); // this algorithm overrides A
13431  __gelsy(mA2, mB, *this, tol, rank);
13432  return *this;
13433  }
13434 
13488  TR tol = basic_cvmMachSp<TR>()) const throw(cvmexception)
13489  {
13490  _check_ne(CVM_SIZESMISMATCH, this->msize(), vB.size());
13491  basic_rmatrix mA(*this); // this algorithm overrides A
13492  basic_rmatrix mB(vB, vB.size(), 1);
13493  basic_rmatrix mX;
13494  __gelsy(mA, mB, mX, tol, rank);
13495  const TR* pResult = mX;
13496  return basic_rvector<TR>(pResult, mX.msize());
13497  }
13498 
13558  TR tol = basic_cvmMachSp<TR>()) const throw(cvmexception)
13559  {
13560  return _gels_sd(true, mB, sv, rank, tol);
13561  }
13562 
13624  tint& rank, TR tol = basic_cvmMachSp<TR>()) throw(cvmexception)
13625  {
13626  _gels_sd(true, mA, mB, sv, rank, tol);
13627  return *this;
13628  }
13629 
13689  tint& rank, TR tol = basic_cvmMachSp<TR>()) const throw(cvmexception)
13690  {
13691  return _gels_sd(true, vB, sv, rank, tol);
13692  }
13693 
13753  TR tol = basic_cvmMachSp<TR>()) const throw(cvmexception)
13754  {
13755  return _gels_sd(false, mB, sv, rank, tol);
13756  }
13757 
13819  tint& rank, TR tol = basic_cvmMachSp<TR>()) throw(cvmexception)
13820  {
13821  _gels_sd(false, mA, mB, sv, rank, tol);
13822  return *this;
13823  }
13824 
13885  tint& rank, TR tol = basic_cvmMachSp<TR>()) const throw(cvmexception)
13886  {
13887  return _gels_sd(false, vB, sv, rank, tol);
13888  }
13889 
13927  tint rank(TR tol = basic_cvmMachSp<TR>()) const throw(cvmexception)
13928  {
13929  tint nRank = 0;
13930  RVector vS(_cvm_min<tint>(this->msize(), this->nsize()));
13931  this->_svd(vS, nullptr, nullptr);
13932  vS.normalize();
13933  for (; nRank < vS.size(); ++nRank) {
13934  if (vS[nRank * this->incr() + CVM0] < tol) break;
13935  }
13936  return nRank;
13937  }
13938 
13939 // QR factorization
13940 // Case 1: "economy" mode, A is (m x n) and Q is (m x n) and R is (n x n)
13989  void qr(basic_rmatrix<TR>& mQ, basic_srmatrix<TR>& mR) const throw(cvmexception) {
13990  this->_qr_rs(mQ, mR);
13991  }
13992 
13993 // Case 2: full mode, A is (m x n) and Q is (m x m) and R is (m x n)
14041  void qr(basic_srmatrix<TR>& mQ, basic_rmatrix<TR>& mR) const throw(cvmexception) {
14042  this->_qr_sr(mQ, mR);
14043  }
14044 
14045 // LQ factorization
14046 // Case 1: "economy" mode, A is (m x n) and L is (m x m) and Q is (m x n)
14093  void lq(basic_srmatrix<TR>& mL, basic_rmatrix<TR>& mQ) const throw(cvmexception) {
14094  this->_lq_sr(mL, mQ);
14095  }
14096 
14097 // Case 2: full mode, A is (m x n) and L is (m x n) and Q is (n x n)
14144  void lq(basic_rmatrix<TR>& mL, basic_srmatrix<TR>& mQ) const throw(cvmexception) {
14145  this->_lq_rs(mL, mQ);
14146  }
14147 
14148 // RQ factorization
14149 // Case 1: "economy" mode, A is (m x n) and R is (m x m) and Q is (m x n)
14191  void rq(basic_srmatrix<TR>& mR, basic_rmatrix<TR>& mQ) const throw(cvmexception) {
14192  this->_rq_sr(mR, mQ);
14193  }
14194 
14195 // Case 2: full mode, A is (m x n) and R is (m x n) and Q is (n x n)
14237  void rq(basic_rmatrix<TR>& mR, basic_srmatrix<TR>& mQ) const throw(cvmexception) {
14238  this->_rq_rs(mR, mQ);
14239  }
14240 
14241 // QL factorization
14242 // Case 1: "economy" mode, A is (m x n) and Q is (m x n) and L is (n x n)
14283  void ql(basic_rmatrix<TR>& mQ, basic_srmatrix<TR>& mL) const throw(cvmexception) {
14284  this->_ql_rs(mQ, mL);
14285  }
14286 
14287 // Case 2: full mode, A is (m x n) and Q is (m x m) and L is (m x n)
14328  void ql(basic_srmatrix<TR>& mQ, basic_rmatrix<TR>& mL) const throw(cvmexception) {
14329  this->_ql_sr(mQ, mL);
14330  }
14331 
14332 // this += alpha * v_col * v_row (rank-1 update)
14392  basic_rmatrix& ger(TR alpha, const RVector& vCol, const RVector& vRow) throw(cvmexception)
14393  {
14394  this->_check_ger();
14395  _check_ne(CVM_SIZESMISMATCH, this->msize(), vCol.size());
14396  _check_ne(CVM_SIZESMISMATCH, this->nsize(), vRow.size());
14397  __ger<TR,basic_rmatrix, RVector>(*this, vCol, vRow, alpha);
14398  return *this;
14399  }
14400 
14401 // this = alpha * m1 * m2 + beta * this
14458  basic_rmatrix& gemm(const basic_rmatrix& m1, bool bTrans1, const basic_rmatrix& m2,
14459  bool bTrans2, TR dAlpha, TR dBeta) throw(cvmexception)
14460  {
14461  this->_check_gemm();
14462  _check_ne(CVM_SIZESMISMATCH, this->msize(), bTrans1 ? m1.nsize() : m1.msize());
14463  _check_ne(CVM_SIZESMISMATCH, this->nsize(), bTrans2 ? m2.msize() : m2.nsize());
14464  _check_ne(CVM_SIZESMISMATCH, bTrans1 ? m1.msize() : m1.nsize(), bTrans2 ? m2.nsize() : m2.msize());
14465  this->_gemm(bTrans1, m1, bTrans2, m2, dAlpha, dBeta);
14466  return *this;
14467  }
14468 
14469 // this = alpha*a*b + beta*this or this = alpha*b*a + beta*this where a is symmetric
14541  basic_rmatrix& symm(bool bLeft, const basic_srsmatrix<TR>& ms, const basic_rmatrix& m,
14542  TR dAlpha, TR dBeta) throw(cvmexception)
14543  {
14544  this->_check_symm();
14545  _check_ne(CVM_SIZESMISMATCH, this->msize(), m.msize());
14546  _check_ne(CVM_SIZESMISMATCH, this->nsize(), m.nsize());
14547  _check_ne(CVM_SIZESMISMATCH, bLeft ? this->msize() : this->nsize(), ms.msize());
14548  this->_symm(bLeft, ms, m, dAlpha, dBeta);
14549  return *this;
14550  }
14551 
14580  basic_rmatrix& vanish()
14581  {
14582  this->_vanish();
14583  return *this;
14584  }
14585 
14611  basic_rmatrix& randomize(TR dFrom, TR dTo)
14612  {
14613  this->_randomize(dFrom, dTo);
14614  return *this;
14615  }
14616 
14617  // 2-norm (maximum singular value)
14618  TR norm2() const override
14619  {
14620  RVector vS(_cvm_min<tint>(this->msize(), this->nsize()));
14621  this->_svd(vS, nullptr, nullptr);
14622  return vS[CVM0];
14623  }
14624 
14626  // ?gemm routines perform a matrix-matrix operation with general matrices. The operation is defined as
14627  // c := alpha*op(a)*op(b) + beta*c,
14628  // where: op(x) is one of op(x) = x or op(x) = x' or op(x) = conjg(x'),
14629  void _gemm(bool bTrans1, const basic_rmatrix& m1, bool bTrans2, const basic_rmatrix& m2, TR dAlpha, TR dBeta) throw(cvmexception) { // this = m1 * m2
14630  basic_rmatrix mTmp1, mTmp2;
14631  const TR* pD1 = m1.get();
14632  const TR* pD2 = m2.get();
14633  if (this->get() == pD1) mTmp1 << m1;
14634  if (this->get() == pD2) mTmp2 << m2;
14635  __gemm<TR,basic_rmatrix>(this->get() == pD1 ? mTmp1 : m1, bTrans1, this->get() == pD2 ? mTmp2 : m2,
14636  bTrans2, dAlpha, *this, dBeta);
14637  }
14638 
14639  // this = alpha*a*b + beta*this or this = alpha*b*a + beta*this where a is symmetric
14640  void _symm(bool bLeft, const basic_srsmatrix<TR>& ms, const basic_rmatrix& m, TR dAlpha, TR dBeta) throw(cvmexception)
14641  {
14642  basic_rmatrix mTmp;
14643  basic_srsmatrix<TR> msTmp;
14644  const TR* pD1 = ms.get();
14645  const TR* pD2 = m._pd();
14646  if (this->get() == pD1) msTmp << ms;
14647  if (this->get() == pD2) mTmp << m;
14648  __symm<TR,basic_srsmatrix<TR>, basic_rmatrix>(bLeft, this->get() == pD1 ? msTmp : ms,
14649  this->get() == pD2 ? mTmp : m, dAlpha, *this, dBeta);
14650  }
14651 
14652  // singular values in decreasing order
14653  virtual void _svd(RVector& vRes, basic_srmatrix<TR>* pmU, basic_srmatrix<TR>* pmVH) const throw(cvmexception)
14654  {
14655  if (pmU != nullptr && pmVH != nullptr) {
14656  _check_ne(CVM_SIZESMISMATCH, this->msize(), pmU->msize());
14657  _check_ne(CVM_SIZESMISMATCH, this->nsize(), pmVH->msize());
14658  }
14659  __svd<TR,basic_rmatrix, basic_srmatrix<TR> >(vRes, vRes.size(), vRes.incr(), *this, pmU, pmVH);
14660  }
14661 
14662  virtual void _pinv(basic_rmatrix& mX, TR threshold) const throw(cvmexception) {
14663  __pinv<TR,basic_rmatrix, basic_rmatrix>(mX, *this, threshold);
14664  }
14665 
14666  virtual void _check_submatrix() const {
14667  }
14668 
14669 protected:
14670  // protected constructors for inherited stuff
14671  basic_rmatrix(tint nM, tint nN, tint nLD, bool bZeroMemory)
14672  : BaseMatrix(nM, nN, nLD, bZeroMemory)
14673  {
14674  }
14675 
14676  basic_rmatrix(tint size, tint incr, tint msize, tint nsize, tint ld)
14677  : BaseMatrix(size, incr, msize, nsize, ld)
14678  {
14679  }
14680 
14681  // non-const version shares memory
14682  basic_rmatrix(TR* pd, tint nM, tint nN, tint nLD, tint nSize)
14683  : BaseMatrix(pd, nM, nN, nLD, nSize)
14684  {
14685  }
14686 
14687  // const version makes a copy
14688  basic_rmatrix(const TR* pd, tint nM, tint nN, tint nLD, tint nSize)
14689  : BaseMatrix(pd, nM, nN, nLD, nSize)
14690  {
14691  }
14692 
14693  // returns diagonal which IS l-value (shares memory)
14694  // 0 - main, negative - low, positive - up
14695  virtual RVector _diag(tint nDiag) throw(cvmexception)
14696  {
14697  tint nShift = 0;
14698  tint nSize = 0;
14699  this->_diag_helper(nDiag, nShift, nSize);
14700  return RVector(this->get() + nShift, nSize, this->ld() + 1);
14701  }
14702 
14703  // returns diagonal which IS NOT l-value (creates a copy)
14704  // 0 - main, negative - low, positive - up
14705  virtual const RVector _diag(tint nDiag) const throw(cvmexception)
14706  {
14707  tint nShift = 0;
14708  tint nSize = 0;
14709  this->_diag_helper(nDiag, nShift, nSize);
14710  return RVector(this->get() + nShift, nSize, this->ld() + 1);
14711  }
14712 
14713  // compares matrix elements (equal sizes assumed)
14714  bool _mequals(const basic_rmatrix& m) const {
14715  return ((*this) - m).norminf() <= basic_cvmMachMin<TR>();
14716  }
14717 
14718  // ?gemv routines perform a matrix-vector operation defined as
14719  // vRes = alpha*m*v + beta * vRes or vRes = alpha*v'*m + beta * vRes
14720  // not virtual since __gemv calls all virtual methods inside
14721  void _gemv(bool bLeft, TR dAlpha, const RVector& v, TR dBeta, RVector& vRes) const
14722 {
14723  RVector vTmp;
14724  basic_rmatrix mTmp;
14725  const TR* pDv = v;
14726  if (vRes.get() == pDv) vTmp << v;
14727  if (vRes.get() == this->get()) mTmp << *this;
14728  __gemv<TR,basic_rmatrix, RVector>(bLeft, vRes.get() == this->get() ? mTmp : *this, dAlpha,
14729  vRes.get() == pDv ? vTmp : v, dBeta, vRes);
14730  }
14731 
14732  // 0-based, returns l-value sharing memory
14733  virtual RVector _row(tint m) {
14734  return RVector(this->get() + m, this->nsize(), this->ld());
14735  }
14736 
14737  // 0-based, returns NOT l-value (copies memory)
14738  // looks like cut-n-paste, but it's not
14739  virtual const RVector _row(tint m) const {
14740  return RVector(this->get() + m, this->nsize(), this->ld());
14741  }
14742 
14743  // 0-based, returns l-value sharing memory
14744  virtual RVector _col(tint n) {
14745  return RVector(this->get() + this->ld() * n, this->msize());
14746  }
14747 
14748  // 0-based, returns NOT l-value (copies memory)
14749  // looks like cut-n-paste, but it's not
14750  virtual const RVector _col(tint n) const {
14751  return RVector(this->get() + this->ld() * n, this->msize());
14752  }
14753 
14754  virtual void _mult(const basic_rmatrix& m1, const basic_rmatrix& m2) throw(cvmexception)
14755  {
14756  _check_ne(CVM_SIZESMISMATCH, this->msize(), m1.msize());
14757  _check_ne(CVM_SIZESMISMATCH, this->nsize(), m2.nsize());
14758  _check_ne(CVM_SIZESMISMATCH, m1.nsize(), m2.msize());
14759  static const TR zero = TR(0.);
14760  static const TR one = TR(1.);
14761  this->_gemm(false, m1, false, m2, one, zero);
14762  }
14763 
14764  virtual void _multiply(RVector& vRes, const RVector& v, bool bLeft) const
14765  {
14766  static const TR zero = TR(0.);
14767  static const TR one = TR(1.);
14768  this->_gemv(bLeft, one, v, zero, vRes);
14769  }
14770 
14771  virtual void _randomize(TR dFrom, TR dTo)
14772  {
14773  if (this->_continuous()) {
14774  __randomize<TR>(this->get(), this->size(), this->incr(), dFrom, dTo);
14775  } else for (tint i = 0; i < this->nsize(); ++i) {
14776  __randomize<TR>(this->get() + this->ld() * i, this->msize(), this->incr(), dFrom, dTo);
14777  }
14778  }
14779 
14780  // QR factorization
14781  // Case 1: "economy" mode, A is (m x n) and Q is (m x n) and R is (n x n)
14782  void _qr_rs(basic_rmatrix<TR>& mQ, basic_srmatrix<TR>& mR) const throw(cvmexception)
14783  {
14784  _check_ne(CVM_SIZESMISMATCH, this->msize(), mQ.msize());
14785  _check_ne(CVM_SIZESMISMATCH, this->nsize(), mQ.nsize());
14786  _check_ne(CVM_SIZESMISMATCH, this->nsize(), mR.msize());
14787  __qre<basic_rmatrix, basic_srmatrix<TR> >(*this, mQ, mR);
14788  }
14789 
14790  // Case 2: full mode, A is (m x n) and Q is (m x m) and R is (m x n)
14791  void _qr_sr(basic_srmatrix<TR>& mQ, basic_rmatrix<TR>& mR) const throw(cvmexception)
14792  {
14793  _check_ne(CVM_SIZESMISMATCH, this->msize(), mQ.msize());
14794  _check_ne(CVM_SIZESMISMATCH, this->msize(), mR.msize());
14795  _check_ne(CVM_SIZESMISMATCH, this->nsize(), mR.nsize());
14796  __qrf<basic_rmatrix, basic_srmatrix<TR> >(*this, mQ, mR);
14797  }
14798 
14799  // RQ factorization
14800  // Case 1: "economy" mode, A is (m x n) and R is (m x m) and Q is (m x n)
14801  void _rq_sr(basic_srmatrix<TR>& mR, basic_rmatrix<TR>& mQ) const throw(cvmexception)
14802  {
14803  _check_gt(CVM_SIZESMISMATCH_GT, this->msize(), this->nsize());
14804  _check_ne(CVM_SIZESMISMATCH, this->msize(), mQ.msize());
14805  _check_ne(CVM_SIZESMISMATCH, this->nsize(), mQ.nsize());
14806  _check_ne(CVM_SIZESMISMATCH, this->msize(), mR.msize());
14807  __rqe<basic_rmatrix, basic_srmatrix<TR> >(*this, mR, mQ);
14808  }
14809 
14810  // Case 2: full mode, A is (m x n) and R is (m x n) and Q is (n x n)
14811  void _rq_rs(basic_rmatrix<TR>& mR, basic_srmatrix<TR>& mQ) const throw(cvmexception)
14812  {
14813  _check_gt(CVM_SIZESMISMATCH_GT, this->msize(), this->nsize());
14814  _check_ne(CVM_SIZESMISMATCH, this->nsize(), mQ.msize());
14815  _check_ne(CVM_SIZESMISMATCH, this->msize(), mR.msize());
14816  _check_ne(CVM_SIZESMISMATCH, this->nsize(), mR.nsize());
14817  __rqf<basic_rmatrix, basic_srmatrix<TR> >(*this, mR, mQ);
14818  }
14819 
14820  // LQ factorization
14821  // Case 1: "economy" mode, A is (m x n) and L is (m x m) and Q is (m x n)
14822  void _lq_sr(basic_srmatrix<TR>& mL, basic_rmatrix<TR>& mQ) const throw(cvmexception)
14823  {
14824  _check_ne(CVM_SIZESMISMATCH, this->msize(), mL.msize());
14825  _check_ne(CVM_SIZESMISMATCH, this->msize(), mQ.msize());
14826  _check_ne(CVM_SIZESMISMATCH, this->nsize(), mQ.nsize());
14827  __lqe<basic_rmatrix, basic_srmatrix<TR> >(*this, mL, mQ);
14828  }
14829 
14830  // Case 2: full mode, A is (m x n) and L is (m x n) and Q is (n x n)
14831  void _lq_rs(basic_rmatrix<TR>& mL, basic_srmatrix<TR>& mQ) const throw(cvmexception)
14832  {
14833  _check_ne(CVM_SIZESMISMATCH, this->msize(), mL.msize());
14834  _check_ne(CVM_SIZESMISMATCH, this->nsize(), mL.nsize());
14835  _check_ne(CVM_SIZESMISMATCH, this->nsize(), mQ.nsize());
14836  __lqf<basic_rmatrix, basic_srmatrix<TR> >(*this, mL, mQ);
14837  }
14838 
14839  // QL factorization
14840  // Case 1: "economy" mode, A is (m x n) and Q is (m x n) and L is (n x n)
14841  void _ql_rs(basic_rmatrix<TR>& mQ, basic_srmatrix<TR>& mL) const throw(cvmexception)
14842  {
14843  _check_lt(CVM_SIZESMISMATCH_LT, this->msize(), this->nsize());
14844  _check_ne(CVM_SIZESMISMATCH, this->msize(), mQ.msize());
14845  _check_ne(CVM_SIZESMISMATCH, this->nsize(), mQ.nsize());
14846  _check_ne(CVM_SIZESMISMATCH, this->nsize(), mL.msize());
14847  __qle<basic_rmatrix, basic_srmatrix<TR> >(*this, mQ, mL);
14848  }
14849 
14850  // Case 2: full mode, A is (m x n) and Q is (m x m) and L is (m x n)
14851  void _ql_sr(basic_srmatrix<TR>& mQ, basic_rmatrix<TR>& mL) const throw(cvmexception)
14852  {
14853  _check_lt(CVM_SIZESMISMATCH_LT, this->msize(), this->nsize());
14854  _check_ne(CVM_SIZESMISMATCH, this->msize(), mQ.msize());
14855  _check_ne(CVM_SIZESMISMATCH, this->msize(), mL.msize());
14856  _check_ne(CVM_SIZESMISMATCH, this->nsize(), mL.nsize());
14857  __qlf<basic_rmatrix, basic_srmatrix<TR> >(*this, mQ, mL);
14858  }
14859 
14860  virtual void _check_ger() { }
14861  virtual void _check_rank1update() { }
14862  virtual void _check_gemm() { }
14863  virtual void _check_symm() { }
14864 
14865 private:
14866  basic_rmatrix& _solve_helper(const basic_srmatrix<TR>& mA, const basic_rmatrix& mB,
14867  TR& dErr, int transp_mode) throw(cvmexception)
14868  {
14869  _check_ne(CVM_SIZESMISMATCH, this->msize(), mA.msize());
14870  _check_ne(CVM_SIZESMISMATCH, this->nsize(), mB.nsize());
14871  _check_ne(CVM_SIZESMISMATCH, mA.msize(), mB.msize());
14872  mA._solve(mB, *this, dErr, nullptr, nullptr, transp_mode);
14873  return *this;
14874  }
14875 
14876  // helper for svd and divide&conquer methods
14877  basic_rmatrix _gels_sd(bool svd, const basic_rmatrix& mB, basic_rvector<TR>& sv,
14878  tint& rank, TR tol) const throw(cvmexception)
14879  {
14880  _check_ne(CVM_SIZESMISMATCH, this->msize(), mB.msize());
14881  _check_ne(CVM_SIZESMISMATCH, sv.size(), _cvm_min<tint>(this->msize(), this->nsize()));
14882  basic_rmatrix mA(*this); // this algorithm overrides A
14883  basic_rmatrix mX;
14884  basic_rvector<TR> sv1(sv.size()); // to ensure that incr=1
14885  if (svd) {
14886  __gelss(mA, mB, mX, tol, sv1, rank);
14887  } else {
14888  __gelsd(mA, mB, mX, tol, sv1, rank);
14889  }
14890  sv = sv1;
14891  return mX;
14892  }
14893 
14894  void _gels_sd(bool svd, const basic_rmatrix& mA, const basic_rmatrix& mB, basic_rvector<TR>& sv,
14895  tint& rank, TR tol) throw(cvmexception)
14896  {
14897  _check_ne(CVM_SIZESMISMATCH, this->msize(), mA.nsize());
14898  _check_ne(CVM_SIZESMISMATCH, this->nsize(), mB.nsize());
14899  _check_ne(CVM_SIZESMISMATCH, mA.msize(), mB.msize());
14900  _check_ne(CVM_SIZESMISMATCH, sv.size(), _cvm_min<tint>(mA.msize(), mA.nsize()));
14901  basic_rmatrix mA2(mA); // this algorithm overrides A
14902  basic_rvector<TR> sv1(sv.size()); // to ensure that incr=1
14903  if (svd) {
14904  __gelss(mA2, mB, *this, tol, sv1, rank);
14905  } else {
14906  __gelsd(mA2, mB, *this, tol, sv1, rank);
14907  }
14908  sv = sv1;
14909  }
14910 
14911  basic_rvector<TR> _gels_sd(bool svd, const basic_rvector<TR>& vB, basic_rvector<TR>& sv,
14912  tint& rank, TR tol) const throw(cvmexception)
14913  {
14914  _check_ne(CVM_SIZESMISMATCH, this->msize(), vB.size());
14915  _check_ne(CVM_SIZESMISMATCH, _cvm_min<tint>(this->msize(), this->nsize()), sv.size());
14916  basic_rmatrix mA(*this); // this algorithm overrides A
14917  basic_rmatrix mB(vB, vB.size(), 1);
14918  basic_rmatrix mX;
14919  basic_rvector<TR> sv1(sv.size()); // to ensure that incr=1
14920  if (svd) {
14921  __gelss(mA, mB, mX, tol, sv1, rank);
14922  } else {
14923  __gelsd(mA, mB, mX, tol, sv1, rank);
14924  }
14925  sv = sv1;
14926  const TR* pResult = mX;
14927  return basic_rvector<TR>(pResult, mX.msize());
14928  }
14930 };
14931 
14932 
14939 template<typename TR>
14940 class basic_srmatrix : public basic_rmatrix<TR>, public SqMatrix<TR,TR>
14941 {
14942  typedef std::complex<TR> TC;
14943  typedef basic_rvector<TR> RVector;
14944  typedef basic_cvector<TR,TC> CVector;
14945  typedef basic_array<TR,TR> BaseArray;
14946  typedef Matrix<TR,TR> BaseMatrix;
14947  typedef SqMatrix<TR,TR> BaseSqMatrix;
14948  typedef basic_rmatrix<TR> BaseRMatrix;
14949 
14950 public:
14977  {
14978  }
14979 
15006  explicit basic_srmatrix(tint nDim)
15007  : BaseRMatrix(nDim, nDim)
15008  {
15009  }
15010 
15041  basic_srmatrix(TR* pd, tint nDim)
15042  : BaseRMatrix(pd, nDim, nDim)
15043  {
15044  }
15045 
15077  basic_srmatrix(const TR* pd, tint nDim)
15078  : BaseRMatrix(pd, nDim, nDim)
15079  {
15080  }
15081 
15111  : BaseRMatrix(m.msize(), m.nsize(), m.msize(), false), BaseSqMatrix()
15112  {
15113  this->_massign(m);
15114  }
15115 
15134  : BaseRMatrix(std::move(m))
15135  {
15136  }