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"
43 # define TINT_ZERO (0)
45 # define CVM_TINT_FORMAT "%d"
49 #if defined (CVM_ZERO_BASED)
50 # define CVM0 TINT_ZERO
52 # define CVM0 TINT_ONE
55 #if !defined (CVM_NO_MT) // as of 5.7 it's by default
57 # if !defined (_PTHREADS)
63 #if defined (_MSC_VER)
65 # define WIN32_LEAN_AND_MEAN // Exclude rarely-used stuff from Windows headers
67 # define _WIN32_WINNT 0x500 // at least Win2000 is required for InitializeCriticalSectionAndSpinCount
74 # if (_MSC_VER < 1400)
75 # error "Please use stable version 5.2 for older MSVC compilers"
77 # if (_MSC_VER >= 1700) // since VC11 aka MS Visual Studio 2012
78 # define CVM_STD_MUTEX
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
87 # include <initializer_list>
89 # if (!defined(__INTEL_COMPILER) || !defined(_WIN64)) && !defined(CVM_ACML) && !(_MSC_VER >= 1500 && defined(_WIN64))
90 # define CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES
92 # if defined(CVM_ACML)
93 # define CVM_COMPLEX_NUMBER_RETURNED
96 # define noexcept throw()
98 # pragma warning(disable:4290)
99 # if defined (CVM_FLOAT)
100 # pragma warning(disable:4244)
102 # if defined (SRC_EXPORTS) && !defined (CVM_EXPORTS)
105 # include <unordered_map>
106 # define CVM_BLOCKS_MAP std::unordered_map
111 # define CVM_API __declspec(dllexport)
113 # define CVM_API __declspec(dllimport)
118 typedef __int64 CVM_LONGEST_INT;
121 typedef unsigned long long CVM_PTR_WRAPPER;
123 # if defined (CVM_ILP64)
124 # error "CVM_ILP64 is incompatible with 32 bit mode"
126 typedef unsigned long CVM_PTR_WRAPPER;
129 # define CVM_VSNPRINTF vsnprintf_s
130 # define CVM_VSNPRINTF_S_DEFINED
131 # define CVM_STRCPY_S_DEFINED
134 #elif defined (__GNUC__)
135 # if defined (__MINGW32__)
136 # define WIN32_LEAN_AND_MEAN
137 # include <windows.h>
138 # include <process.h>
140 # define CVM_USE_CRITICAL_SECTION_NOT_MUTEX
147 # include <semaphore.h>
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
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>
172 # define CVM_STDEXT stdext
173 # define CVM_BLOCKS_MAP std::map
175 typedef long long CVM_LONGEST_INT;
177 # if defined(__AMD64__) || defined(_WIN64) || defined(__x86_64__)
178 typedef unsigned long long CVM_PTR_WRAPPER;
180 # if defined (CVM_ILP64)
181 # error "CVM_ILP64 is incompatible with 32 bit mode"
183 typedef unsigned long CVM_PTR_WRAPPER;
186 # define CVM_VSNPRINTF vsnprintf
194 # if defined(__clang__) && __clang_major__ <= 5
195 # define CVM_NO_DEFAULT_RANDOM_ENGINE_SUPPORTED 1
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
205 # error "Unsupported compiler"
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
234 #if !defined(CVM_NO_DEFAULT_RANDOM_ENGINE_SUPPORTED)
238 #if defined (STLPORT)
239 # define CVM_USES_STLPORT
242 #if defined (_DEBUG) || defined (DEBUG)
245 # define CVM_ASSERT(p,n) _cvm_assert(p,n);
247 # define CVM_ASSERT(p,n)
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
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
286 #define CVM_THE_LAST_ERROR_CODE 37
287 #define CVM_MATRIX_ELEMENT_SEPARATOR " "
288 #define CVM_EOL std::endl
292 #ifdef CVM_NO_NAMESPACE
293 # define CVM_NAMESPACE_BEG
294 # define CVM_NAMESPACE_END
296 # define CVM_NAMESPACE_BEG namespace cvm {
297 # define CVM_NAMESPACE_END }
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;
328 template<
typename T,
typename TR>
class type_proxy;
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;
340 std::string msUnknown;
345 CVM_API
const std::string&
_get(
int nException);
346 CVM_API
bool _add(
int nNewCause,
const char* szNewMessage);
365 virtual const char* _get_message(
int nCause)
const {
420 const char*
what()
const throw()
override {
444 static bool add(
int nNewCause,
const char* szNewMessage) {
451 CVM_API
void __copy(
tint nSize,
const T* pFrom,
tint nFromIncr, T* pTo,
tint nToIncr);
453 CVM_API
void __swap(
tint nSize, T* p1,
tint n1Incr, T* p2,
tint n2Incr);
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);
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);
467 template<
typename TR,
typename TC>
468 CVM_API TR __norm(
const TC* pd,
tint nSize,
tint nIncr);
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);
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);
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);
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);
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);
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,
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>
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);
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);
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);
552 template<
typename TM>
553 CVM_API
tint __cholesky(TM& m);
554 template<
typename TM>
556 template<
typename TR,
typename TM,
typename TV>
557 CVM_API
void __poequ(
const TM& m, TV& vScalings, TR& dCond, TR& dMax);
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);
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);
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);
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);
594 inline const T& _cvm_min(
const T& x,
const T& y) {
595 return x < y ? x : y;
599 inline const T& _cvm_max(
const T& x,
const T& y) {
600 return x > y ? x : y;
604 inline const TR* __get_real_p(
const std::complex<TR>* c)
606 #if defined (CVM_USES_STLPORT)
609 return reinterpret_cast<const TR*
>(c);
614 inline const TR* __get_imag_p(
const std::complex<TR>* c)
616 #if defined (CVM_USES_STLPORT)
619 return reinterpret_cast<const TR*
>(c) + 1;
624 inline TR* __get_real_p(std::complex<TR>* c)
626 #if defined (CVM_USES_STLPORT)
629 return reinterpret_cast<TR*
>(c);
634 inline TR* __get_imag_p(std::complex<TR>* c)
636 #if defined (CVM_USES_STLPORT)
639 return reinterpret_cast<TR*
>(c) + 1;
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; }
652 static char mchars[15];
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;}
673 inline TR basic_cvmMachMin() {
674 static const TR _min = (std::numeric_limits<TR>::min)();
679 inline TR basic_cvmMachSp() {
680 static const TR _eps = (std::numeric_limits<TR>::epsilon)();
684 template<
typename TR>
685 inline TR basic_cvmMachMax() {
686 static const TR _max((std::numeric_limits<TR>::max)());
690 template<
typename TR>
691 inline TR basic_cvmLogMachMax() {
692 static const TR _log_max(::log(basic_cvmMachMax<TR>()));
698 void _check_negative(
int err_code, T val)
throw(
cvmexception)
700 static T zero = T(0);
707 void _check_positive(
int err_code, T val)
throw(
cvmexception)
709 static T zero = T(0);
716 void _check_positive(
int err_code, T val,
const char* func,
const char* file,
int line)
throw(
cvmexception)
718 static T zero = T(0);
725 void _check_lt(
int err_code, T idx, T limit)
throw(
cvmexception)
733 void _check_le(
int err_code, T idx, T limit)
throw(
cvmexception)
741 void _check_ge(
int err_code, T idx, T limit)
throw(
cvmexception)
749 void _check_gt(
int err_code, T idx, T limit)
throw(
cvmexception)
757 void _check_ne(
int err_code, T idx, T limit)
throw(
cvmexception)
765 void _check_lt_gt(
int err_code, T idx, T low_limit, T up_limit)
throw(
cvmexception)
767 if (idx < low_limit || idx > up_limit) {
773 void _check_lt_ge(
int err_code, T idx, T low_limit, T up_limit)
throw(
cvmexception)
775 if (idx < low_limit || idx >= up_limit) {
783 #ifdef CVM_USE_POOL_MANAGER
786 class CVM_API MemoryBlocks
792 BlockProperty(
size_t nSize,
tint nRefCount) : mnSize(nSize), mnRefCount(nRefCount) { }
795 typedef std::map<tbyte*, BlockProperty, std::less<tbyte*> > map_Blocks;
796 typedef map_Blocks::iterator itr_Blocks;
798 typedef std::multimap<size_t, tbyte*> map_FreeBs;
799 typedef map_FreeBs::iterator itr_FreeBs;
801 typedef std::map<tbyte*, itr_FreeBs, std::less<tbyte*> > map_FreeIt;
802 typedef map_FreeIt::iterator itr_FreeIt;
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);
815 void Assert(
const void* pvBlock,
size_t nBytes);
817 void AddNew(
tbyte* pBlock,
size_t nBytes);
824 class CVM_API MemoryPool
826 typedef std::list<tbyte*> list_blocks;
830 void operator () (T* p)
const {
835 list_blocks mOutBlocks;
836 MemoryBlocks mMemoryBlocks;
847 void Assert(
const void* pvBlock,
size_t nBytes) {
848 mMemoryBlocks.Assert(pvBlock, nBytes);
866 return reinterpret_cast<T*
>(_cvmMalloc(nEls *
sizeof(T)));
876 inline T* cvmAddRef(
const T* pd) {
877 return reinterpret_cast<T*
>(_cvmAddRef(reinterpret_cast<const tbyte*>(pd)));
889 return _cvmFree(reinterpret_cast<tbyte*&>(pd));
892 #else // CVM_USE_POOL_MANAGER
905 return nEls > 0 ? ::new T[nEls] :
nullptr;
924 #endif // !CVM_USE_POOL_MANAGER
944 memset(p, 0, nEls *
sizeof(T));
948 CVM_API
void _cvm_assert(
const void* pvBlock,
size_t nBytes);
950 inline float _abs(
const float& v) {
951 return static_cast<float>(fabs(v));
954 inline double _abs(
const double& v) {
958 inline long double _abs(
const long double& v) {
962 inline float _abs(
const std::complex<float>& v) {
963 return static_cast<float>(sqrt(v.real() * v.real() + v.imag() * v.imag()));
966 inline double _abs(
const std::complex<double>& v) {
967 return sqrt(v.real() * v.real() + v.imag() * v.imag());
970 inline long double _abs(
const std::complex<long double>& v) {
971 return sqrt(v.real() * v.real() + v.imag() * v.imag());
975 return v < 0 ? -v : v;
978 inline float _sqrt(
const float& v) {
979 return static_cast<float>(sqrt(v));
982 inline double _sqrt(
const double& v) {
986 inline std::complex<float> _conjugate(
const std::complex<float>& v) {
987 return std::complex<float>(v.real(), -v.imag());
990 inline std::complex<double> _conjugate(
const std::complex<double>& v) {
991 return std::complex<double>(v.real(), -v.imag());
994 template<
typename TR>
995 inline bool _conjugated(
const std::complex<TR>& v1,
const std::complex<TR>& v2, TR tol)
997 return _abs(v1.real() - v2.real()) <= tol &&
998 _abs(v1.imag() + v2.imag()) <= tol;
1004 template<
typename T,
typename TR>
1005 inline std::ostream& operator << (std::ostream& os, const type_proxy<T,TR>& mOut)
1007 #if defined (_MSC_VER) && !defined (CVM_USES_STLPORT)
1008 os.imbue(std::locale::empty());
1010 os << static_cast<T>(mOut);
1015 template<
typename T,
typename TR>
1018 #if defined (_MSC_VER) && !defined (CVM_USES_STLPORT)
1019 is.imbue(std::locale::empty());
1060 template<
typename TR,
typename TC>
1063 #if !defined (CVM_USES_STLPORT) && defined (_MSC_VER)
1064 is.imbue(std::locale::empty());
1068 for (
tint i = 0; i < nSize && is.good(); i += aIn.
incr()) {
1097 template<
typename TR,
typename TC>
1098 std::ostream& operator << (std::ostream& os, const basic_array<TR,TC>& aOut)
1100 #if !defined (CVM_USES_STLPORT) && defined (_MSC_VER)
1101 os.imbue(std::locale::empty());
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()) {
1150 template<
typename TR,
typename TC>
1153 #if defined (_MSC_VER) && !defined (CVM_USES_STLPORT)
1154 is.imbue(std::locale::empty());
1157 for (
tint i = 0; i < mIn.
msize(); ++i) {
1158 for (
tint j = 0; j < mIn.
nsize() && is.good(); ++j) {
1160 mIn._ij_proxy_val(i, j) = v;
1193 template<
typename TR,
typename TC>
1194 std::ostream& operator << (std::ostream& os, const Matrix<TR,TC>& mOut)
1196 #if defined (_MSC_VER) && !defined (CVM_USES_STLPORT)
1197 os.imbue(std::locale::empty());
1199 for (
tint i = 0; i < mOut.msize(); ++i) {
1200 for (
tint j = 0; j < mOut.nsize() && os.good(); ++j) {
1209 template<
typename TR,
typename TC,
typename RM,
typename RBM>
1210 inline void _copy_b_matrix(RM& m, RBM& mb,
bool bLeftToRight)
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;
1221 for (
tint i = 0; i < nN; ++i) {
1231 if (nN - i <= nKL) {
1232 nS -= nKL + 1 - (nN - i);
1235 pL = m.get() + i * nM + nShiftL;
1236 pR = mb._pb() + i * nCol + nShiftR;
1239 bLeftToRight ? pL : pR,
1241 bLeftToRight ? pR : pL,
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)
1253 static const TR two(2.);
1254 __scal<TR,TC>(pd, nSize, nIncr, two);
1256 __add<TC>(pd, nSize, nIncr, p2, nIncr2);
1260 __add<TC>(pd, nSize, nIncr, p1, nIncr1);
1262 __copy<TC>(nSize, p1, nIncr1, pd, nIncr);
1263 __add<TC>(pd, nSize, nIncr, p2, nIncr2);
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)
1275 static const TR zero(0.);
1276 __scal<TR,TC>(pd, nSize, nIncr, zero);
1278 __subtract<TC>(pd, nSize, nIncr, p2, nIncr2);
1282 static const TR mone(-1.);
1283 __subtract<TC>(pd, nSize, nIncr, p1, nIncr1);
1284 __scal<TR,TC>(pd, nSize, nIncr, mone);
1286 __copy<TC>(nSize, p1, nIncr1, pd, nIncr);
1287 __subtract<TC>(pd, nSize, nIncr, p2, nIncr2);
1293 template<
typename TR,
typename TC>
1294 inline void _incr(TC* pDm,
tint mnSize,
tint mnIncr,
const TC* pd,
tint nIncr)
1297 static const TR two(2.);
1298 __scal<TR,TC>(pDm, mnSize, mnIncr, two);
1300 __add<TC>(pDm, mnSize, mnIncr, pd, nIncr);
1305 template<
typename TR,
typename TC>
1306 inline void _decr(TC* pDm,
tint mnSize,
tint mnIncr,
const TC* pd,
tint nIncr)
1309 static const TR zero(0.);
1310 __scal<TR,TC>(pDm, mnSize, mnIncr, zero);
1312 __subtract<TC>(pDm, mnSize, mnIncr, pd, nIncr);
1317 template<
typename TR,
typename TC>
1318 void _set_real(TC* pDm,
tint mnSize,
tint mnIncr, TR d)
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) {
1330 template<typename TR, typename TC>
1331 void _set_imag(TC* pDm,
tint mnSize,
tint mnIncr, TR d)
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) {
1349 template<typename T, typename TR>
1364 type_proxy(T& ref,
bool read_only)
1366 mbReadOnly(read_only)
1377 type_proxy(
const T& ref,
bool read_only =
true)
1378 : mT(const_cast<T&>(ref)),
1379 mbReadOnly(read_only)
1389 type_proxy(
const type_proxy& p)
1391 mbReadOnly(p.mbReadOnly)
1401 operator T()
const {
1446 const T&
get()
const {
1457 type_proxy& operator = (
const type_proxy& p)
1496 const T* operator &()
const {
1507 template<
typename U>
1522 template<
typename U>
1537 template<
typename U>
1552 template<
typename U>
1567 template<
typename U>
1579 template<
typename U>
1591 template<
typename U>
1603 template<
typename U>
1670 return _real<T,TR>(mT);
1681 return _imag<T,TR>(mT);
1686 template<
typename T,
typename TR>
1692 template<
typename T,
typename TR>
1698 template<
typename T,
typename TR>
1704 template<
typename T,
typename TR>
1710 template<
typename U,
typename T,
typename TR>
1712 return T(u) + p.
val();
1716 template<
typename U,
typename T,
typename TR>
1718 return T(u) - p.
val();
1722 template<
typename U,
typename T,
typename TR>
1724 return T(u) * p.
val();
1728 template<
typename U,
typename T,
typename TR>
1730 return T(u) / p.
val();
1735 template<
typename T>
1739 #if !defined(CVM_NO_DEFAULT_RANDOM_ENGINE_SUPPORTED)
1740 std::random_device mre;
1741 std::uniform_int_distribution<int> mDist;
1749 #if defined(CVM_NO_DEFAULT_RANDOM_ENGINE_SUPPORTED)
1750 : mMax(static_cast<T>(RAND_MAX))
1752 srand(static_cast<unsigned int>(time(
nullptr)));
1755 : mMax(static_cast<T>((std::numeric_limits<int>::max)())),
1757 mDist(0, (std::numeric_limits<int>::max)())
1766 T _get(T dFrom, T dTo) noexcept
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;
1774 return dMin +
static_cast<T
>(mDist(mre)) * (dMax - dMin) / mMax;
1795 static T
get(T dFrom, T dTo) noexcept
1798 return r._get(dFrom, dTo);
1817 template<
typename TR,
typename TC>
1824 #ifdef CVM_USE_POOL_MANAGER
1866 #ifdef CVM_USE_POOL_MANAGER
1898 #ifdef CVM_USE_POOL_MANAGER
1906 if (bZeroMemory) cvmZeroMemory<TC>(this->
get(),
msz);
1938 #ifdef CVM_USE_POOL_MANAGER
1939 mpd(cvmAddRef<TC>(pd))
1975 #ifdef CVM_USE_POOL_MANAGER
1984 __copy<TC>(
msz, pd, nIncr, this->
get(), this->
incr());
2012 #ifdef CVM_USE_POOL_MANAGER
2044 #ifdef CVM_USE_POOL_MANAGER
2051 this->_assign(a.
get(), a.
incr());
2074 #ifdef CVM_USE_POOL_MANAGER
2081 _move(std::move(a));
2109 this->_assign(a.get(), a.incr());
2129 _move(std::move(a));
2141 #ifdef CVM_USE_POOL_MANAGER
2202 #ifdef CVM_USE_POOL_MANAGER
2206 return this->
mpf ==
nullptr ? this->
mp.get() : this->
mpf;
2242 const TC*
get()
const
2244 #ifdef CVM_USE_POOL_MANAGER
2248 return this->
mpf ==
nullptr ? this->
mp.get() : this->
mpf;
2322 operator const TC* ()
const {
2357 return this->
at(static_cast<size_type>(n -
CVM0));
2391 return this->
at(static_cast<size_type>(n -
CVM0));
2425 return this->
at(static_cast<size_type>(n -
CVM0));
2459 return this->
at(static_cast<size_type>(n -
CVM0));
2486 this->_assign(p, 1);
2547 this->_resize(nNewSize);
2606 return this->_indofmax();
2636 return this->_indofmin();
2675 return __norm<TR,TC>(this->
get(), this->
size(), this->
incr());
2714 return _abs(this->
get()[(this->_indofmax() -
CVM0) * this->
incr()]);
2756 for (
tint i = 0; i < nSize; i += this->
incr()) {
2757 dNorm += _abs(this->
get()[i]);
2803 return this->
norm();
2807 const tint* _pincr()
const {
2808 return &this->
mincr;
2811 const tint* _psize()
const {
2818 static const TR one(1.);
2819 if (_abs(d) <= basic_cvmMachMin<TR>()) {
2822 this->_scalr(one / d);
2831 virtual const TC* _pd()
const {
2882 for (
tint i = 0; i < (
tint)n; ++i) {
2883 this->
get()[i] = val;
2893 for (
tint i = 0; i < n; ++i) {
2894 this->
get()[i] = *(begin + i);
2907 __swap<TC>(this->
size(), this->
get(), this->
incr(), v.get(), v.incr());
2914 CVM_ASSERT(this->
get(), (n + 1) *
sizeof(TC))
2915 return this->
get()[n * this->
incr()];
2922 CVM_ASSERT(this->
get(), (n + 1) *
sizeof(TC))
2923 return this->
get()[n * this->
incr()];
2955 this->_resize(this->
size() + 1);
2956 this->
get()[this->
size() - 1] = x;
2986 if (this->
size() > 0) this->_resize(this->
size() - 1);
3028 this->_resize(this->
size() + 1);
3029 for (
tint i = this->
size() - 1; i > n; --i) {
3030 this->
get()[i] = this->
get()[i - 1];
3075 for (
tint i = n; i < this->
size() - 1; ++i) {
3076 this->
get()[i] = this->
get()[i + 1];
3078 if (this->
size() > 0) this->_resize(this->
size() - 1);
3087 #ifdef CVM_USE_POOL_MANAGER
3099 #ifdef CVM_USE_POOL_MANAGER
3103 if (a.mpf ==
nullptr) {
3104 if (
mpf ==
nullptr) {
3105 mp = std::move(a.mp);
3110 this->_assign(a.mp.get(), a.incr());
3114 this->_resize(a.size());
3115 this->_assign(a.mpf, a.incr());
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>()) {
3143 const TR dNorm = this->
norm();
3144 if (dNorm > basic_cvmMachMin<TR>()) {
3145 static const TR one(1.);
3146 this->_scalr(one / dNorm);
3153 #ifdef CVM_USE_POOL_MANAGER
3154 cvmFree<TC>(this->mpd);
3155 this->mpd = cvmMalloc<TC>(a.
size());
3158 this->
mpf =
nullptr;
3160 this->msz = a.
size();
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);
3176 #ifdef CVM_USE_POOL_MANAGER
3177 cvmFree<TC>(this->mpd);
3181 this->
mpf =
nullptr;
3183 this->msz = nNewSize;
3189 bool _is_empty()
const
3191 #ifdef CVM_USE_POOL_MANAGER
3192 return this->mpd ==
nullptr;
3194 return this->
mpf ==
nullptr && this->
mp.get() ==
nullptr;
3199 virtual void _scalr(TR d) {
3200 __scal<TR,TC>(this->
get(), this->
size(), this->
incr(), d);
3204 virtual tint _indofmax()
const {
3205 return __idamax<TC>(this->
get(), this->
size(), this->
incr());
3209 virtual tint _indofmin()
const {
3210 return __idamin<TC>(this->
get(), this->
size(), this->
incr());
3214 virtual void _set(TC d)
3218 for (
tint i = 0; i < nSize; i += this->incr()) {
3223 virtual void _assign(
const TC* pd,
tint nIncr)
3225 if (this->
get() != pd) {
3226 __copy<TC>(this->
size(), pd, nIncr, this->
get(), this->
incr());
3230 virtual void _assign_shifted(TC* pDshifted,
const TC* pd,
tint nSize,
tint nIncr,
tint)
3232 if (pDshifted != pd) {
3233 __copy<TC>(nSize, pd, nIncr, pDshifted, this->
incr());
3271 friend std::istream& operator >> <> (std::istream& is,
basic_array& aIn);
3296 friend std::ostream& operator << <> (std::ostream& os,
const basic_array& aOut);
3306 template<
typename TR>
3309 typedef std::complex<TR> TC;
3366 #if defined(CVM_USE_INITIALIZER_LISTS)
3392 TR* p = this->
mp.get();
3393 for (
auto it = list.begin(); it != list.end(); ++it) {
3573 __copy<TR>(this->
size(), v, v.
incr(), this->
get(), this->
incr());
3619 this->_assign(v, v.incr());
3666 this->_assign(pd, nIncr);
3705 this->_assign_shifted(this->
get() + this->
incr() * n, pd, this->
size() - n, nIncr, 0);
3746 this->_assign_shifted(this->
get() + this->
incr() * n, pd, _cvm_min<tint>(nSize, this->
size() - n), nIncr, 0);
3783 return assign(n, v, v.size(), v.incr());
3854 this->_resize(nNewSize);
3888 return this->_equals(v);
3964 __copy<TR>(this->
size(), v.get(), v.incr(), this->
get(), this->
incr());
4009 __add<TR>(vSum.get(), vSum.size(), vSum.incr(), v._pd(), v.incr());
4053 __subtract<TR>(vDiff.get(), vDiff.size(), vDiff.incr(), v._pd(), v.incr());
4098 _sum<TR,TR>(this->
get(), this->
size(), this->
incr(), v1, v1.incr(), v2, v2.incr());
4143 _diff<TR,TR>(this->
get(), this->
size(), this->
incr(), v1, v1.incr(), v2, v2.incr());
4189 _incr<TR,TR>(this->
get(), this->
size(), this->
incr(), v, v.incr());
4235 _decr<TR,TR>(this->
get(), this->
size(), this->
incr(), v, v.incr());
4264 static const TR mone(-1.);
4367 this->_scalr(dMult);
4476 return __dot<TR>(this->
get(), this->
size(), this->
incr(), v.get(), v.incr());
4516 m._multiply(vRes, *
this,
true);
4555 m._multiply(*
this, v,
true);
4594 m._multiply(*
this, v,
false);
4644 static const TR one(1.);
4645 __ger<TR,basic_rmatrix<TR>,
basic_rvector>(mRes, *
this, v, one);
4696 return _solve_helper(mA, vB, dErr, 0);
4747 return _solve_helper(mA, vB, dErr, 1);
4794 return this->
solve(mA, vB, dErr);
4886 return _operator_solve_helper(mA, 1);
4926 return _operator_solve_helper(mA, 0);
4997 mA._solve(vB, *
this, dErr, mLU, pPivots, 0);
5060 return this->
solve_lu(mA, mLU, pPivots, vB, dErr);
5126 __gels(transpose, mA2, mB, mX, vErr);
5186 __gelsy(mA2, mB, mX, tol, rank);
5246 _gels_sd(
true, mA, vB, sv, rank, tol);
5306 _gels_sd(
false, mA, vB, sv, rank, tol);
5384 mA._svd(*
this,
nullptr,
nullptr);
5462 mA._svd(*
this,
nullptr,
nullptr);
5546 mA._svd(*
this, &mU, &mVH);
5630 mA._svd(*
this, &mU, &mVH);
5709 __eig<basic_rvector, basic_srsmatrix<TR>,
basic_srmatrix<TR> >(*
this, mA,
nullptr,
true);
5790 __eig<basic_rvector, basic_srsmatrix<TR>,
basic_srmatrix<TR> >(*
this, mA, &mEigVect,
true);
6012 m._gemv(bLeft, dAlpha, v, dBeta, *
this);
6072 m._gbmv(bLeft, dAlpha, v, dBeta, *
this);
6101 __randomize<TR>(this->
get(), this->
size(), this->
incr(), dFrom, dTo);
6112 mA._solve(vB, *
this, dErr,
nullptr,
nullptr, transp_mode);
6121 mA._solve(*
this, vX, dErr,
nullptr,
nullptr, transp_mode);
6137 __gelss(mA2, mB, mX, tol, sv1, rank);
6139 __gelsd(mA2, mB, mX, tol, sv1, rank);
6154 template<
typename TR,
typename TC>
6215 #if defined(CVM_USE_INITIALIZER_LISTS)
6246 TC* p = this->
mp.get();
6247 for (
auto it = list.begin(); it != list.end(); ++it) {
6424 __copy<TC>(this->
size(), v, v.
incr(), this->
get(), this->
incr());
6489 :
BaseArray(nSize, pRe == nullptr || pIm == nullptr)
6491 __copy2<TR,TC>(this->
get(), this->
size(), this->
incr(), pRe, pIm, nIncrRe, nIncrIm);
6531 __copy2<TR,TC>(this->
get(), this->
size(), this->
incr(), vRe, vIm, vRe.
incr(), vIm.incr());
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);
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());
6703 this->_assign(v, v.incr());
6759 this->_assign(pd, nIncr);
6796 this->_assign_shifted(this->
get() + this->
incr() * n, pd, this->
size() - n, nIncr, 0);
6831 this->_assign_shifted(this->
get() + this->
incr() * n, pd, _cvm_min<tint>(nSize, this->
size() - n), nIncr, 0);
6865 return assign(n, v, v.size(), v.incr());
6926 __copy_real<TR,TC>(this->
get(), this->
size(), this->
incr(), vRe, vRe.incr());
6959 __copy_imag<TR,TC>(this->
get(), this->
size(), this->
incr(), vIm, vIm.incr());
6987 this->_set_real_number(d);
7015 this->_set_imag_number(d);
7056 this->_resize(nNewSize);
7090 return this->_equals(v);
7163 __copy<TC>(this->
size(), v.get(), v.incr(), this->
get(), this->
incr());
7207 __add<TC>(vSum.get(), vSum.size(), vSum.incr(), v._pd(), v.incr());
7250 __subtract<TC>(vDiff.get(), vDiff.size(), vDiff.incr(), v._pd(), v.incr());
7295 _sum<TR,TC>(this->
get(), this->
size(), this->
incr(), v1, v1.incr(), v2, v2.incr());
7340 _diff<TR,TC>(this->
get(), this->
size(), this->
incr(), v1, v1.incr(), v2, v2.incr());
7386 _incr<TR,TC>(this->
get(), this->
size(), this->
incr(), v, v.incr());
7432 _decr<TR,TC>(this->
get(), this->
size(), this->
incr(), v, v.incr());
7461 static const TR mone(-1.);
7629 this->_scalr(dMult);
7692 this->_scalc(cMult);
7827 this->_assign(v, v.incr());
7828 __conj<TC>(this->
get(), this->
size(), this->
incr());
7861 __conj<TC>(this->
get(), this->
size(), this->
incr());
7893 return __dotu<TC>(this->
get(), this->
size(), this->
incr(), v.get(), v.incr());
7926 return __dotc<TC>(this->
get(), this->
size(), this->
incr(), v.get(), v.incr());
7967 m._multiply(vRes, *
this,
true);
8016 m._multiply(*
this, v,
true);
8064 m._multiply(*
this, v,
false);
8114 static const TC one(1., 0.);
8115 __geru<TC, basic_cmatrix<TR,TC>,
basic_cvector>(mRes, *
this, v, one);
8171 static const TC one(1., 0.);
8172 __gerc<TC, basic_cmatrix<TR,TC>,
basic_cvector>(mRes, *
this, v, one);
8222 return _solve_helper(mA, vB, dErr, 0);
8274 return _solve_helper(mA, vB, dErr, 1);
8324 return _solve_helper(mA, vB, dErr, 2);
8369 return this->
solve(mA, vB, dErr);
8508 return _operator_solve_helper(mA, 1);
8548 return _operator_solve_helper(mA, 0);
8620 mA._solve(vB, *
this, dErr, mLU, pPivots, 0);
8683 return this->
solve_lu(mA, mLU, pPivots, vB, dErr);
8759 __gels(conjugate, mA2, mB, mX, vErr);
8821 __gelsy(mA2, mB, mX, tol, rank);
8883 _gels_sd(
true, mA, vB, sv, rank, tol);
8945 _gels_sd(
false, mA, vB, sv, rank, tol);
8997 mA._eig(*
this,
nullptr,
true);
9078 mA._eig(*
this, &mEigVect, bRightVect);
9123 mA._eig(*
this,
nullptr,
true);
9195 mA._eig(*
this, &mEigVect, bRightVect);
9338 __ggev(ma, mb, *
this, vBeta, &mEigVectLeft, &mEigVectRight);
9427 __ggev(ma, mb, *
this, vBeta, bRightVect ?
nullptr : &mEigVect, bRightVect ? &mEigVect :
nullptr);
9486 __ggev<basic_scmatrix<TR,TC>,
basic_scmatrix<TR,TC>,
basic_cvector, basic_cvector>(ma, mb, *
this, vBeta,
nullptr,
nullptr);
9576 __ggev(ma, mb, *
this, vBeta, &mEigVectLeft, &mEigVectRight);
9668 __ggev(ma, mb, *
this, vBeta, bRightVect ?
nullptr : &mEigVect, bRightVect ? &mEigVect :
nullptr);
9729 m._gemv(bLeft, dAlpha, v, dBeta, *
this);
9790 m._gbmv(bLeft, dAlpha, v, dBeta, *
this);
9819 __randomize_real<TC, TR>(this->
get(), this->
size(), this->
incr(), dFrom, dTo);
9848 __randomize_imag<TC, TR>(this->
get(), this->
size(), this->
incr(), dFrom, dTo);
9856 if (_abs(d) <= basic_cvmMachMin<TR>()) {
9859 static const TC one(1., 0.);
9860 this->_scalc(one / d);
9865 __scal<TC, TC>(this->
get(), this->
size(), this->
incr(), d);
9868 void _set_real_number(TR d)
9870 _set_real<TR,TC>(this->
get(), this->
size(), this->
incr(), d);
9873 void _set_imag_number(TR d)
9875 _set_imag<TR,TC>(this->
get(), this->
size(), this->
incr(), d);
9884 mA._solve(vB, *
this, dErr,
nullptr,
nullptr, transp_mode);
9893 mA._solve(*
this, vX, dErr,
nullptr,
nullptr, transp_mode);
9909 __gelss(mA2, mB, mX, tol, sv1, rank);
9911 __gelsd(mA2, mB, mX, tol, sv1, rank);
9926 template<
typename TR,
typename TC>
10016 mm(bBeColumn ? v.
size() : 1),
10017 mn(bBeColumn ? 1 : v.
size()),
10020 __copy<TC>(this->
size(), v, v.
incr(), this->
get(), this->
incr());
10041 #if defined(CVM_USE_DELEGATING_CONSTRUCTORS)
10042 :
Matrix(m.size(), m.incr(), m.msize(), m.nsize(), m.ld())
10050 _mmove(std::move(m));
10070 _mmove(std::move(m));
10166 return (this->_indofmax() -
CVM0) %
mm +
CVM0;
10193 return (this->_indofmin() -
CVM0) %
mm +
CVM0;
10220 return (this->_indofmax() -
CVM0) /
mm +
CVM0;
10247 return (this->_indofmin() -
CVM0) /
mm +
CVM0;
10253 TR rSum, rNorm(0.);
10255 for (j = 0; j < this->
nsize(); ++j) {
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]);
10264 if (rSum > rNorm) {
10274 TR rSum, rNorm(0.);
10276 for (i = 0; i < this->
msize(); ++i) {
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]);
10284 if (rSum > rNorm) {
10292 const tint* _pm()
const {
10296 const tint* _pn()
const {
10300 const tint* _pld()
const {
10304 TC* _sub_pointer_nocheck(
tint row,
tint col) {
10305 return this->_pd() + ((col -
CVM0) * this->
ld() + row -
CVM0);
10308 TC* _sub_pointer(
tint row,
tint col,
tint height,
tint width)
throw(cvmexception)
10316 return _sub_pointer_nocheck(row, col);
10319 virtual tint _ldm()
const {
10323 virtual const tint* _pldm()
const {
10324 return this->_pld();
10327 virtual bool _continuous()
const {
10328 return this->
msize() == this->
ld();
10331 void _check_ld()
const
10333 if (!this->_continuous()) {
10338 void _scalr(TR d)
override
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);
10357 void _mmove(
Matrix&& m) noexcept
10360 #ifdef CVM_USE_POOL_MANAGER
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;
10375 this->_resize2(m.msize(), m.nsize());
10381 void _diag_helper(
tint nDiag,
tint& nShift,
tint& nSize)
const throw(cvmexception)
10385 nShift = nDiag * this->
ld();
10387 this->
nsize() - nDiag;
10392 this->
msize() - nShift;
10396 tint _indofmax()
const override
10399 return __idamax<TC>(this->
get(), this->
size(), this->
incr());
10402 tint _indofmin()
const override
10405 return __idamin<TC>(this->
get(), this->
size(), this->
incr());
10408 void _assign(
const TC* pd,
tint nIncr)
override
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());
10419 void _assign_shifted(TC* pDshifted,
const TC* pd,
tint nRows,
tint nCols,
tint nLD)
override
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());
10428 void _set(TC d)
override
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;
10440 virtual void _massign(
const Matrix& m)
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());
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());
10455 virtual void _resize2(
tint nNewM,
tint nNewN)
throw(cvmexception)
10459 const bool is_empty = this->_is_empty();
10460 if (nNewM != this->
msize() || nNewN != this->
nsize() || is_empty) {
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);
10474 #ifdef CVM_USE_POOL_MANAGER
10475 cvmFree<TC>(this->mpd);
10479 this->
mpf =
nullptr;
10481 this->msz = nNewSize;
10490 virtual
tint _ld_for_replace()
const {
10494 virtual tint _size_for_replace()
const {
10495 return this->
size();
10498 void _replace(
const Matrix& m)
throw(cvmexception)
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());
10507 this->msz = m._size_for_replace();
10509 this->
mpf =
nullptr;
10515 this->
mld = m._ld_for_replace();
10518 void _transp_m(const
Matrix& m)
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);
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());
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);
10534 virtual TC _ij_val(
tint i,
tint j)
const {
10535 CVM_ASSERT(this->
get(), (this->
ld() * j + 1) *
sizeof(TC))
10536 return this->
get()[this->
ld() * j + i];
10539 virtual
void _swap_rows(
tint n1,
tint n2) throw(cvmexception)
10544 __swap<TC>(this->
nsize(), this->
get() + n1 -
CVM0, this->
ld(), this->
get() + n2 -
CVM0, this->
ld());
10548 virtual void _swap_cols(
tint n1,
tint n2)
throw(cvmexception)
10553 __swap<TC>(this->
msize(), this->
get() + (n1 -
CVM0) * this->
ld(), 1, this->
get() + (n2 -
CVM0) * this->
ld(), 1);
10557 virtual const TC* _pp(
const Matrix& m)
const {
10562 virtual void _vanish()
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));
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());
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());
10590 void _mincr(
const Matrix& m)
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());
10599 void _mdecr(
const Matrix& m)
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());
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);
10622 template<
typename TR,
typename TC>
10638 virtual tint _size()
const = 0;
10639 virtual tint _msize()
const = 0;
10640 virtual tint _nsize()
const = 0;
10641 virtual tint _ld()
const = 0;
10644 virtual const TC* _pv()
const = 0;
10645 virtual TC* _pv() = 0;
10650 const tint mm = this->_msize();
10651 const tint mld = this->_ld();
10652 TC* pd = this->_pv();
10654 const tint nM1 = mld + 1, nM1m = mld - 1, nM2m = mm - 1;
10655 tint i = 1, j = 1, m;
10658 __swap<TC>(m, pd + j, 1, pd + j + nM1m,
mld);
10669 void _sq_plus_plus()
10671 TC* pd = this->_pv();
10672 static const TC one(1.);
10673 const tint nSize = this->_size();
10674 const tint nNext = this->_msize() + 1;
10676 for (
tint i = 0; i < nSize; i += nNext) {
10682 void _sq_minus_minus()
10684 TC* pd = this->_pv();
10685 static const TC one(1.);
10686 const tint nSize = this->_size();
10687 const tint nNext = this->_msize() + 1;
10689 for (
tint i = 0; i < nSize; i += nNext) {
10695 void _clean_low_triangle()
10697 const tint mm = this->_msize();
10698 const tint mld = this->_ld();
10699 TC* pd = this->_pv();
10701 static const TR zero(0.);
10702 for (
tint i = 1; i <
mm; ++i) {
10703 __scal<TR,TC>(pd + n, mm - i, 1, zero);
10716 template<
typename TR>
10719 typedef std::complex<TR> TC;
10888 :
BaseMatrix(m.msize(), m.nsize(), m.msize(), false)
10985 :
BaseMatrix(m._sub_pointer(nRow, nCol, nHeight, nWidth), nHeight, nWidth, m.
ld(), nHeight * nWidth)
10987 m._check_submatrix();
11032 return this->_ij_proxy_val(nRow -
CVM0, nCol -
CVM0);
11069 return this->_ij_val(nRow -
CVM0, nCol -
CVM0);
11110 return this->_col(nCol -
CVM0);
11152 return this->_row(nRow -
CVM0);
11184 return this->_col(nCol -
CVM0);
11216 return this->_row(nRow -
CVM0);
11253 return this->_diag(nDiag);
11289 return this->_diag(nDiag);
11385 this->_assign(v, v.incr());
11418 this->_assign(pd, 1);
11458 this->_assign_shifted(this->_sub_pointer_nocheck(nRow, nCol), m._pd(), m.
msize(), m.
nsize(), m.
ld());
11535 this->_resize2(nNewM, nNewN);
11777 this->_msum(m1, m2);
11823 this->_mdiff(m1, m2);
11942 static const TR mone(-1.);
11972 mRes._scalr(dMult);
12039 this->_scalr(dMult);
12108 this->_normalize();
12142 mRes._transp_m(*
this);
12182 if (this->
get() == m.
get()) {
12184 this->_transp_m(mTmp);
12186 this->_transp_m(m);
12223 this->_resize2(this->
nsize(), this->
msize());
12224 this->_transp_m(mTmp);
12261 this->_multiply(vRes, v,
false);
12298 mRes.
mult(*
this, m);
12338 this->_mult(m1, m2);
12396 static const TR one(1.);
12397 this->_check_rank1update();
12401 __ger<TR,basic_rmatrix, RVector>(*
this, vCol, vRow, one);
12440 this->_swap_rows(n1, n2);
12477 this->_swap_cols(n1, n2);
12523 return _solve_helper(mA, mB, dErr, 0);
12575 return _solve_helper(mA, mB, dErr, 1);
12617 static TR dErr(0.);
12618 return this->
solve(mA, mB, dErr);
12668 static TR dErr(0.);
12748 mA._solve(mB, *
this, dErr, mLU, pPivots, 0);
12821 static TR dErr(0.);
12822 return this->
solve_lu(mA, mLU, pPivots, mB, dErr);
12870 this->_svd(vRes,
nullptr,
nullptr);
12952 this->_svd(vRes, &mU, &mVH);
13012 this->_pinv(mAx, threshold);
13077 mA._pinv(*
this, threshold);
13152 __gels(transpose, mA, mB, mX, vErr);
13230 __gels(transpose, mA2, mB, *
this, vErr);
13302 __gels(transpose, mA, mB, mX, vErr);
13304 const TR* pResult = mX;
13361 TR tol = basic_cvmMachSp<TR>())
const throw(cvmexception)
13366 __gelsy(mA, mB, mX, tol, rank);
13425 tint& rank, TR tol = basic_cvmMachSp<TR>()) throw(cvmexception)
13431 __gelsy(mA2, mB, *
this, tol, rank);
13488 TR tol = basic_cvmMachSp<TR>())
const throw(cvmexception)
13494 __gelsy(mA, mB, mX, tol, rank);
13495 const TR* pResult = mX;
13558 TR tol = basic_cvmMachSp<TR>())
const throw(cvmexception)
13560 return _gels_sd(
true, mB, sv, rank, tol);
13624 tint& rank, TR tol = basic_cvmMachSp<TR>()) throw(cvmexception)
13626 _gels_sd(
true, mA, mB, sv, rank, tol);
13689 tint& rank, TR tol = basic_cvmMachSp<TR>())
const throw(cvmexception)
13691 return _gels_sd(
true, vB, sv, rank, tol);
13753 TR tol = basic_cvmMachSp<TR>())
const throw(cvmexception)
13755 return _gels_sd(
false, mB, sv, rank, tol);
13819 tint& rank, TR tol = basic_cvmMachSp<TR>()) throw(cvmexception)
13821 _gels_sd(
false, mA, mB, sv, rank, tol);
13885 tint& rank, TR tol = basic_cvmMachSp<TR>())
const throw(cvmexception)
13887 return _gels_sd(
false, vB, sv, rank, tol);
13927 tint rank(TR tol = basic_cvmMachSp<TR>())
const throw(cvmexception)
13931 this->_svd(vS,
nullptr,
nullptr);
13933 for (; nRank < vS.
size(); ++nRank) {
13934 if (vS[nRank * this->
incr() +
CVM0] < tol)
break;
13990 this->_qr_rs(mQ, mR);
14042 this->_qr_sr(mQ, mR);
14094 this->_lq_sr(mL, mQ);
14145 this->_lq_rs(mL, mQ);
14192 this->_rq_sr(mR, mQ);
14238 this->_rq_rs(mR, mQ);
14284 this->_ql_rs(mQ, mL);
14329 this->_ql_sr(mQ, mL);
14394 this->_check_ger();
14397 __ger<TR,basic_rmatrix, RVector>(*
this, vCol, vRow, alpha);
14459 bool bTrans2, TR dAlpha, TR dBeta)
throw(cvmexception)
14461 this->_check_gemm();
14465 this->_gemm(bTrans1, m1, bTrans2, m2, dAlpha, dBeta);
14542 TR dAlpha, TR dBeta)
throw(cvmexception)
14544 this->_check_symm();
14548 this->_symm(bLeft, ms, m, dAlpha, dBeta);
14613 this->_randomize(dFrom, dTo);
14621 this->_svd(vS,
nullptr,
nullptr);
14629 void _gemm(
bool bTrans1,
const basic_rmatrix& m1,
bool bTrans2,
const basic_rmatrix& m2, TR dAlpha, TR dBeta)
throw(cvmexception) {
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);
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);
14655 if (pmU !=
nullptr && pmVH !=
nullptr) {
14659 __svd<TR,basic_rmatrix, basic_srmatrix<TR> >(vRes, vRes.size(), vRes.incr(), *
this, pmU, pmVH);
14662 virtual void _pinv(
basic_rmatrix& mX, TR threshold)
const throw(cvmexception) {
14663 __pinv<TR,basic_rmatrix, basic_rmatrix>(mX, *
this, threshold);
14666 virtual void _check_submatrix()
const {
14672 : BaseMatrix(nM, nN, nLD, bZeroMemory)
14677 : BaseMatrix(size, incr, msize, nsize, ld)
14683 : BaseMatrix(pd, nM, nN, nLD, nSize)
14689 : BaseMatrix(pd, nM, nN, nLD, nSize)
14695 virtual RVector _diag(
tint nDiag)
throw(cvmexception)
14699 this->_diag_helper(nDiag, nShift, nSize);
14700 return RVector(this->
get() + nShift, nSize, this->
ld() + 1);
14705 virtual const RVector _diag(
tint nDiag)
const throw(cvmexception)
14709 this->_diag_helper(nDiag, nShift, nSize);
14710 return RVector(this->
get() + nShift, nSize, this->
ld() + 1);
14715 return ((*
this) - m).norminf() <= basic_cvmMachMin<TR>();
14721 void _gemv(
bool bLeft, TR dAlpha,
const RVector& v, TR dBeta, RVector& vRes)
const
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);
14733 virtual RVector _row(
tint m) {
14734 return RVector(this->
get() + m, this->
nsize(), this->
ld());
14739 virtual const RVector _row(
tint m)
const {
14740 return RVector(this->
get() + m, this->
nsize(), this->
ld());
14744 virtual RVector _col(
tint n) {
14745 return RVector(this->
get() + this->
ld() * n, this->
msize());
14750 virtual const RVector _col(
tint n)
const {
14751 return RVector(this->
get() + this->
ld() * n, this->
msize());
14759 static const TR zero = TR(0.);
14760 static const TR one = TR(1.);
14761 this->_gemm(
false, m1,
false, m2, one, zero);
14764 virtual void _multiply(RVector& vRes,
const RVector& v,
bool bLeft)
const
14766 static const TR zero = TR(0.);
14767 static const TR one = TR(1.);
14768 this->_gemv(bLeft, one, v, zero, vRes);
14771 virtual void _randomize(TR dFrom, TR dTo)
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);
14787 __qre<basic_rmatrix, basic_srmatrix<TR> >(*
this, mQ, mR);
14796 __qrf<basic_rmatrix, basic_srmatrix<TR> >(*
this, mQ, mR);
14807 __rqe<basic_rmatrix, basic_srmatrix<TR> >(*
this, mR, mQ);
14817 __rqf<basic_rmatrix, basic_srmatrix<TR> >(*
this, mR, mQ);
14827 __lqe<basic_rmatrix, basic_srmatrix<TR> >(*
this, mL, mQ);
14836 __lqf<basic_rmatrix, basic_srmatrix<TR> >(*
this, mL, mQ);
14847 __qle<basic_rmatrix, basic_srmatrix<TR> >(*
this, mQ, mL);
14857 __qlf<basic_rmatrix, basic_srmatrix<TR> >(*
this, mQ, mL);
14860 virtual void _check_ger() { }
14861 virtual void _check_rank1update() { }
14862 virtual void _check_gemm() { }
14863 virtual void _check_symm() { }
14867 TR& dErr,
int transp_mode)
throw(cvmexception)
14872 mA._solve(mB, *
this, dErr,
nullptr,
nullptr, transp_mode);
14878 tint& rank, TR tol)
const throw(cvmexception)
14886 __gelss(mA, mB, mX, tol, sv1, rank);
14888 __gelsd(mA, mB, mX, tol, sv1, rank);
14895 tint& rank, TR tol)
throw(cvmexception)
14904 __gelss(mA2, mB, *
this, tol, sv1, rank);
14906 __gelsd(mA2, mB, *
this, tol, sv1, rank);
14912 tint& rank, TR tol)
const throw(cvmexception)
14921 __gelss(mA, mB, mX, tol, sv1, rank);
14923 __gelsd(mA, mB, mX, tol, sv1, rank);
14926 const TR* pResult = mX;
14939 template<
typename TR>
14942 typedef std::complex<TR> TC;
15170 :
BaseRMatrix(m.msize(), m.nsize(), m.msize(), false)
15204 __copy<TR>(this->
msize(), v, v.
incr(), this->
get(), this->
msize() + 1);
15238 m._check_submatrix();
15245 return this->_ij_proxy_val(nRow -
CVM0, nCol -
CVM0);
15252 return this->_ij_val(nRow -
CVM0, nCol -
CVM0);
15259 return this->_col(nCol -
CVM0);
15266 return this->_col(nCol -
CVM0);
15273 return this->_row(nRow -
CVM0);
15280 return this->_row(nRow -
CVM0);
15347 this->_assign(v, v.incr());
15354 this->_assign(pd, 1);
15395 this->_assign_shifted(this->_sub_pointer_nocheck(nRow, nCol), m._pd(), m.msize(), m.nsize(), m.ld());
15443 this->_resize2(nNewDim, nNewDim);
15622 this->_msum(m1, m2);
15667 this->_mdiff(m1, m2);
15788 static const TR mone(-1.);
15821 this->_plus_plus();
15852 this->_plus_plus();
15883 this->_minus_minus();
15914 this->_minus_minus();
15944 mRes._scalr(dMult);
15987 this->_scalr(dMult);
15999 this->_normalize();
16088 return this->transpose();
16170 mRes.mult(*
this, m);
16214 this->mult(mTmp, m);
16220 this->_swap_rows(n1, n2);
16226 this->_swap_cols(n1, n2);
16292 return _solve_helper(vB, dErr, 0);
16357 return _solve_helper(vB, dErr, 1);
16417 static TR dErr(0.);
16418 return this->solve(vB, dErr);
16479 static TR dErr(0.);
16480 return this->solve_tran(vB, dErr);
16543 return _solve_helper(mB, dErr, 0);
16606 return _solve_helper(mB, dErr, 1);
16666 static TR dErr(0.);
16667 return this->solve(mB, dErr);
16729 static TR dErr(0.);
16730 return this->solve_tran(mB, dErr);
16771 return vB / (*this);
16812 return vB % (*this);
16904 this->_solve(vB, vX, dErr, mLU, pPivots, 0);
16992 static TR dErr(0.);
16993 return this->solve_lu(mLU, pPivots, vB, dErr);
17084 this->_solve(mB, mX, dErr, mLU, pPivots, 0);
17171 static TR dErr(0.);
17172 return this->solve_lu(mLU, pPivots, mB, dErr);
17205 TR
det()
const throw(cvmexception) {
17206 return this->_det();
17287 this->_low_up(nPivots);
17364 mRes._low_up(nPivots);
17409 TR cond()
const throw(cvmexception)
17412 __cond_num<TR,basic_srmatrix>(*
this, dCondNum);
17459 __inv<basic_srmatrix>(*
this, m);
17506 __inv<basic_srmatrix>(mRes, *
this);
17573 __exp<basic_srmatrix, TR>(*
this, m, tol);
17638 __exp<basic_srmatrix, TR>(mRes, *
this, tol);
17710 if (v.incr() > 1) v1 << v;
17711 __polynom<TR,RVector>(this->
get(), this->
ld(), this->
msize(), m._pd(), m._ldm(), v.
incr() > 1 ? v1 : v);
17782 if (v.
incr() > 1) v1 << v;
17783 __polynom<TR,RVector>(mRes.
get(), mRes.
ld(), this->
msize(), this->
get(), this->
ld(), v.
incr() > 1 ? v1 : v);
17855 this->_eig(vEig, &mEigVect, bRightVect);
17901 this->_eig(vEig,
nullptr,
true);
17949 this->_check_cholesky();
17951 tint nOutInfo = __cholesky<basic_srmatrix>(*this);
17954 this->_clean_low_triangle();
17989 this->_check_bunch_kaufman();
17991 __bunch_kaufman<basic_srmatrix>(*
this, nPivots);
18023 this->_qr_ss(mQ, mR);
18054 this->_lq_ss(mL, mQ);
18085 this->_ql_ss(mQ, mL);
18116 this->_rq_ss(mR, mQ);
18148 this->_plus_plus();
18188 this->_randomize(dFrom, dTo);
18193 virtual void _eig(CVector& vEig,
basic_scmatrix<TR,TC>* mEigVect,
bool bRightVect)
const throw(cvmexception)
18195 __eig<CVector, basic_srmatrix, basic_scmatrix<TR,TC> >(vEig, *
this, mEigVect, bRightVect);
18198 virtual void _solve(
const RVector& vB, RVector& vX,
18199 TR& dErr,
const TR* pLU,
const tint* pPivots,
int transp_mode)
const throw(cvmexception)
18204 if (vB.
incr() > 1) vB1 << vB;
18205 if (vX.incr() > 1) vX1 << vX;
18206 __solve<TR,TR, basic_srmatrix>(*
this, 1, vB.
incr() > 1 ? vB1 : vB, vB.
size(),
18207 vX.incr() > 1 ? vX1 : vX, vX.size(), dErr, pLU, pPivots, transp_mode);
18208 if (vX.incr() > 1) vX = vX1;
18211 virtual void _solve(
const BaseRMatrix& mB, BaseRMatrix& mX,
18212 TR& dErr,
const TR* pLU,
const tint* pPivots,
int transp_mode)
const throw(cvmexception)
18215 __solve<TR,TR, basic_srmatrix>(*
this, mB.
nsize(), mB, mB.
ld(), mX, mX.ld(), dErr, pLU, pPivots, transp_mode);
18218 const TR* _pv()
const override {
return this->
get();}
18219 TR* _pv()
override {
return this->
get();}
18224 : BaseRMatrix(nDim, nDim, nLD, bZeroMemory)
18229 : BaseRMatrix(size, incr, msize, nsize, ld)
18235 : BaseRMatrix(pd, nDim, nDim, nLD, nSize)
18241 : BaseRMatrix(pd, nDim, nDim, nLD, nSize)
18245 tint _size()
const override {
return this->
size();}
18246 tint _msize()
const override {
return this->
msize();}
18247 tint _nsize()
const override {
return this->
nsize();}
18248 tint _ld()
const override {
return this->
ld();}
18252 RVector _diag(
tint nDiag)
throw(cvmexception)
override
18254 const tint nD = _abs(nDiag);
18256 return RVector(this->
get() + (nDiag > 0 ? nDiag * this->
ld() : nD), this->
msize() - nD, this->
ld() + 1);
18261 const RVector _diag(
tint nDiag)
const throw(cvmexception)
override
18263 const tint nD = _abs(nDiag);
18265 return RVector(this->
get() + (nDiag > 0 ? nDiag * this->
ld() : nD), this->
msize() - nD, this->
ld() + 1);
18276 virtual void _transp() {
18277 this->_sq_transp();
18280 virtual void _plus_plus() {
18281 this->_sq_plus_plus();
18284 virtual void _minus_minus() {
18285 this->_sq_minus_minus();
18288 virtual TR _det()
const throw(cvmexception)
18291 switch (this->
msize()) {
18295 dDet = this->_ij_val(0, 0);
18298 dDet = this->_ij_val(0, 0) * this->_ij_val(1, 1) -
18299 this->_ij_val(1, 0) * this->_ij_val(0, 1);
18303 static const TR one(1.);
18305 RVector vUpDiag = this->_low_up_diag(naPivots);
18309 dDet *= vUpDiag[i];
18310 if (i + (1 -
CVM0) != naPivots[i]) {
18315 catch (
const cvmexception& e) {
18323 virtual void _low_up(
tint* nPivots)
throw(cvmexception) {
18324 __low_up<basic_srmatrix>(*
this, nPivots);
18359 virtual void _check_cholesky() { }
18360 virtual void _check_bunch_kaufman() { }
18363 RVector _solve_helper(
const RVector& vB, TR& dErr,
int transp_mode)
const throw(cvmexception)
18366 RVector vX(this->
msize());
18367 this->_solve(vB, vX, dErr,
nullptr,
nullptr, transp_mode);
18371 BaseRMatrix _solve_helper(
const BaseRMatrix& mB, TR& dErr,
int transp_mode)
const throw(cvmexception)
18375 this->_solve(mB, mX, dErr,
nullptr,
nullptr, transp_mode);
18388 template<
typename TR,
typename TC>
18558 :
BaseMatrix(m.msize(), m.nsize(), m.msize(), false)
18646 :
BaseMatrix(m.msize(), m.nsize(), m.msize(), true)
18649 __copy2<TR,TC>(this->
get(), this->
size(), this->
incr(), m._pd(),
nullptr);
18652 __copy2<TR,TC>(this->
get(), this->
size(), this->
incr(),
nullptr, m._pd());
18698 __copy2<TR,TC>(this->
get(), this->
size(), this->
incr(), pRe, pIm, 1, 1);
18732 :
BaseMatrix(mRe.msize(), mRe.nsize(), mRe.msize(), false)
18737 __copy2<TR,TC>(this->
get(), this->
size(), this->
incr(), mRe, mIm, mRe.
incr(), mIm.incr());
18770 :
BaseMatrix(m._sub_pointer(nRow, nCol, nHeight, nWidth), nHeight, nWidth, m.
ld(), nHeight * nWidth)
18772 m._check_submatrix();
18815 return this->_ij_proxy_val(nRow -
CVM0, nCol -
CVM0);
18850 return this->_ij_val(nRow -
CVM0, nCol -
CVM0);
18890 return this->_col(nCol -
CVM0);
18930 return this->_row(nRow -
CVM0);
18961 return this->_col(nCol -
CVM0);
18992 return this->_row(nRow -
CVM0);
19040 return this->_diag(nDiag);
19077 return this->_diag(nDiag);
19107 __copy<TR>(this->
size(), __get_real_p<TR>(this->
get()), this->
incr() * 2, mRet, mRet.
incr());
19138 __copy<TR>(this->
size(), __get_imag_p<TR>(this->
get()), this->
incr() * 2, mRet, mRet.
incr());
19227 this->_assign(v, v.incr());
19258 this->_assign(pd, 1);
19298 this->_assign_shifted(this->_sub_pointer_nocheck(nRow, nCol), m._pd(), m.
msize(), m.
nsize(), m.
ld());
19356 __copy_real<TR,TC>(this->
get(), this->
size(), this->
incr(), mRe._pd(), mRe.incr());
19387 __copy_imag<TR,TC>(this->
get(), this->
size(), this->
incr(), mIm._pd(), mIm.incr());
19413 this->_set_real_number(d);
19439 this->_set_imag_number(d);
19487 this->_resize2(nNewM, nNewN);
19723 this->_msum(m1, m2);
19770 this->_mdiff(m1, m2);
19887 static const TR mone(-1.);
19919 mRes._scalr(dMult);
19986 mRes._scalc(cMult);
20052 this->_scalr(dMult);
20117 this->_scalc(cMult);
20183 this->_normalize();
20225 mRes._transp_m(*
this);
20267 mRes._transp_m(*
this);
20268 __conj<TC>(mRes.
get(), mRes.
size(), mRes.
incr());
20313 if (this->
get() == m.
get()) {
20315 this->_transp_m(mTmp);
20317 this->_transp_m(m);
20361 this->transpose(m);
20362 __conj<TC>(this->
get(), this->
size(), this->
incr());
20405 this->_resize2(this->
nsize(), this->
msize());
20406 this->_transp_m(mTmp);
20449 __conj<TC>(this->
get(), this->
size(), this->
incr());
20482 this->_multiply(vRes, v,
false);
20514 mRes.
mult(*
this, m);
20546 this->_mult(m1, m2);
20595 static const TC one(1., 0.);
20596 this->_check_rank1update_u();
20600 __geru<TC, basic_cmatrix, CVector>(*
this, vCol, vRow, one);
20650 static const TC one(1., 0.);
20651 this->_check_rank1update_c();
20655 __gerc<TC, basic_cmatrix, CVector>(*
this, vCol, vRow, one);
20690 this->_swap_rows(n1, n2);
20723 this->_swap_cols(n1, n2);
20771 return _solve_helper(mA, mB, dErr, 0);
20819 return _solve_helper(mA, mB, dErr, 1);
20867 return _solve_helper(mA, mB, dErr, 2);
20910 static TR dErr(0.);
20911 return this->
solve(mA, mB, dErr);
20957 static TR dErr(0.);
21004 static TR dErr(0.);
21084 mA._solve(mB, *
this, dErr, mLU, pPivots, 0);
21157 static TR dErr(0.);
21158 return this->
solve_lu(mA, mLU, pPivots, mB, dErr);
21230 this->_svd(vRes,
nullptr,
nullptr);
21308 this->_svd(vRes, &mU, &mVH);
21371 this->_pinv(mAx, threshold);
21436 mA._pinv(*
this, threshold);
21514 __gels(conjugate, mA, mB, mX, vErr);
21595 __gels(conjugate, mA2, mB, *
this, vErr);
21671 __gels(conjugate, mA, mB, mX, vErr);
21673 const TC* pResult = mX;
21733 TR tol = basic_cvmMachSp<TR>())
const throw(cvmexception)
21738 __gelsy(mA, mB, mX, tol, rank);
21800 tint& rank, TR tol = basic_cvmMachSp<TR>()) throw(cvmexception)
21806 __gelsy(mA2, mB, *
this, tol, rank);
21866 TR tol = basic_cvmMachSp<TR>())
const throw(cvmexception)
21872 __gelsy(mA, mB, mX, tol, rank);
21873 const TC* pResult = mX;
21939 TR tol = basic_cvmMachSp<TR>())
const throw(cvmexception)
21941 return _gels_sd(
true, mB, sv, rank, tol);
22008 tint& rank, TR tol = basic_cvmMachSp<TR>()) throw(cvmexception)
22010 _gels_sd(
true, mA, mB, sv, rank, tol);
22076 tint& rank, TR tol = basic_cvmMachSp<TR>())
const throw(cvmexception)
22078 return _gels_sd(
true, vB, sv, rank, tol);
22143 TR tol = basic_cvmMachSp<TR>())
const throw(cvmexception)
22145 return _gels_sd(
false, mB, sv, rank, tol);
22212 tint& rank, TR tol = basic_cvmMachSp<TR>()) throw(cvmexception)
22214 _gels_sd(
false, mA, mB, sv, rank, tol);
22280 tint& rank, TR tol = basic_cvmMachSp<TR>())
const throw(cvmexception)
22282 return _gels_sd(
false, vB, sv, rank, tol);
22316 tint rank(TR tol = basic_cvmMachSp<TR>())
const throw(cvmexception)
22320 this->_svd(vS,
nullptr,
nullptr);
22322 for (; nRank < vS.
size(); ++nRank) {
22323 if (vS[nRank * this->
incr() +
CVM0] < tol)
break;
22378 this->_qr_rs(mQ, mR);
22430 this->_qr_sr(mQ, mR);
22482 this->_lq_sr(mL, mQ);
22533 this->_lq_rs(mL, mQ);
22579 this->_rq_sr(mR, mQ);
22624 this->_rq_rs(mR, mQ);
22669 this->_ql_rs(mQ, mL);
22713 this->_ql_sr(mQ, mL);
22775 this->_check_geru();
22778 __geru<TC, basic_cmatrix, CVector>(*
this, vCol, vRow, alpha);
22841 this->_check_gerc();
22844 __gerc<TC, basic_cmatrix, CVector>(*
this, vCol, vRow, alpha);
22909 bool bConj2, TC cAlpha, TC cBeta)
throw(cvmexception)
22911 this->_check_gemm();
22915 this->_gemm(bConj1, m1, bConj2, m2, cAlpha, cBeta);
22991 TC cAlpha, TC cBeta)
throw(cvmexception)
22993 this->_check_hemm();
22997 this->_hemm(bLeft, ms, m, cAlpha, cBeta);
23061 this->_randomize_real(dFrom, dTo);
23090 this->_randomize_imag(dFrom, dTo);
23098 this->_svd(vS,
nullptr,
nullptr);
23103 void _div(TC d)
throw(cvmexception)
23105 if (_abs(d) <= basic_cvmMachMin<TR>()) {
23108 static const TC one(1., 0.);
23109 this->_scalc(one / d);
23115 void _gemm(
bool bTrans1,
const basic_cmatrix& m1,
bool bTrans2,
const basic_cmatrix& m2, TC dAlpha, TC dBeta)
throw(cvmexception)
23118 const TC* pD1 = m1.
get();
23119 const TC* pD2 = m2.
get();
23120 if (this->
get() == pD1) mTmp1 << m1;
23121 if (this->
get() == pD2) mTmp2 << m2;
23122 __gemm<TC, basic_cmatrix>(this->
get() == pD1 ? mTmp1 : m1, bTrans1,
23123 this->
get() == pD2 ? mTmp2 : m2, bTrans2, dAlpha, *
this, dBeta);
23131 const TC* pD1 = ms.
get();
23132 const TC* pD2 = m._pd();
23133 if (this->
get() == pD1) msTmp << ms;
23134 if (this->
get() == pD2) mTmp << m;
23135 __hemm<TC, basic_schmatrix<TR,TC>,
basic_cmatrix>(bLeft, this->
get() == pD1 ? msTmp : ms,
23136 this->
get() == pD2 ? mTmp : m, dAlpha, *
this, dBeta);
23142 if (pmU !=
nullptr && pmVH !=
nullptr) {
23146 __svd<TR,basic_cmatrix, basic_scmatrix<TR,TC> >(vRes, vRes.size(), vRes.incr(), *
this, pmU, pmVH);
23149 virtual void _pinv(
basic_cmatrix& mX, TR threshold)
const throw(cvmexception) {
23150 __pinv<TR,basic_cmatrix, basic_cmatrix>(mX, *
this, threshold);
23153 virtual void _check_submatrix()
const {
23156 virtual void _scalc(TC d)
23158 if (this->_continuous()) {
23159 __scal<TC, TC>(this->
get(), this->
size(), this->
incr(), d);
23160 }
else for (
tint i = 0; i < this->
nsize(); ++i) {
23161 __scal<TC, TC>(this->
get() + this->
ld() * i, this->
msize(), this->
incr(), d);
23168 : BaseMatrix(nM, nN, nLD, bZeroMemory)
23173 : BaseMatrix(size, incr, msize, nsize, ld)
23179 : BaseMatrix(pd, nM, nN, nLD, nSize)
23185 : BaseMatrix(pd, nM, nN, nLD, nSize)
23189 virtual CVector _diag(
tint nDiag)
throw(cvmexception)
23193 this->_diag_helper(nDiag, nShift, nSize);
23194 return CVector(this->
get() + nShift, nSize, this->
ld() + 1);
23197 virtual const CVector _diag(
tint nDiag)
const throw(cvmexception)
23201 this->_diag_helper(nDiag, nShift, nSize);
23202 return CVector(this->
get() + nShift, nSize, this->
ld() + 1);
23208 return ((*
this) - m).norminf() <= basic_cvmMachMin<TR>();
23214 void _gemv(
bool bLeft, TC dAlpha,
const CVector& v, TC dBeta, CVector& vRes)
const
23219 if (vRes.get() == pDv) vTmp << v;
23220 if (vRes.get() == this->
get()) mTmp << *
this;
23221 __gemv<TC, basic_cmatrix, CVector>(bLeft, vRes.get() == this->
get() ? mTmp : *
this, dAlpha,
23222 vRes.
get() == pDv ? vTmp : v, dBeta, vRes);
23226 virtual CVector _row(
tint m) {
23227 return CVector(this->
get() + m, this->
nsize(), this->
ld());
23231 virtual const CVector _row(
tint m)
const {
23232 return CVector(this->
get() + m, this->
nsize(), this->
ld());
23236 virtual CVector _col(
tint n) {
23237 return CVector(this->
get() + this->
ld() * n, this->
msize());
23241 virtual const CVector _col(
tint n)
const {
23242 return CVector(this->
get() + this->
ld() * n, this->
msize());
23250 static const TR zero = TR(0.);
23251 static const TR one = TR(1.);
23252 this->_gemm(
false, m1,
false, m2, one, zero);
23255 virtual void _multiply(CVector& vRes,
const CVector& v,
bool bLeft)
const
23257 static const TR zero = TR(0.);
23258 static const TR one = TR(1.);
23259 this->_gemv(bLeft, one, v, zero, vRes);
23262 void _set_real_number(TR d)
23264 if (this->_continuous()) {
23265 _set_real<TR,TC>(this->
get(), this->
size(), this->
incr(), d);
23266 }
else for (
tint i = 0; i < this->
nsize(); ++i) {
23267 _set_real<TR,TC>(this->
get() + this->
ld() * i, this->
msize(), this->
incr(), d);
23271 virtual void _set_imag_number(TR d)
23273 if (this->_continuous()) {
23274 _set_imag<TR,TC>(this->
get(), this->
size(), this->
incr(), d);
23275 }
else for (
tint i = 0; i < this->
nsize(); ++i) {
23276 _set_imag<TR,TC>(this->
get() + this->
ld() * i, this->
msize(), this->
incr(), d);
23280 virtual void _randomize_real(TR dFrom, TR dTo)
23282 if (this->_continuous()) {
23283 __randomize_real<TC, TR>(this->
get(), this->
size(), this->
incr(), dFrom, dTo);
23284 }
else for (
tint i = 0; i < this->
nsize(); ++i) {
23285 __randomize_real<TC, TR>(this->
get() + this->
ld() * i, this->
msize(), this->
incr(), dFrom, dTo);
23289 virtual void _randomize_imag(TR dFrom, TR dTo)
23291 if (this->_continuous()) {
23292 __randomize_imag<TC, TR>(this->
get(), this->
size(), this->
incr(), dFrom, dTo);
23293 }
else for (
tint i = 0; i < this->
nsize(); ++i) {
23294 __randomize_imag<TC, TR>(this->
get() + this->
ld() * i, this->
msize(), this->
incr(), dFrom, dTo);
23305 __qre<basic_cmatrix, basic_scmatrix<TR,TC> >(*
this, mQ, mR);
23314 __qrf<basic_cmatrix, basic_scmatrix<TR,TC> >(*
this, mQ, mR);
23325 __rqe<basic_cmatrix, basic_scmatrix<TR,TC> >(*
this, mR, mQ);
23335 __rqf<basic_cmatrix, basic_scmatrix<TR,TC> >(*
this, mR, mQ);
23345 __lqe<basic_cmatrix, basic_scmatrix<TR,TC> >(*
this, mL, mQ);
23354 __lqf<basic_cmatrix, basic_scmatrix<TR,TC> >(*
this, mL, mQ);
23365 __qle<basic_cmatrix, basic_scmatrix<TR,TC> >(*
this, mQ, mL);
23375 __qlf<basic_cmatrix, basic_scmatrix<TR,TC> >(*
this, mQ, mL);
23378 virtual void _check_geru() { }
23379 virtual void _check_gerc() { }
23380 virtual void _check_rank1update_u() { }
23381 virtual void _check_rank1update_c() { }
23382 virtual void _check_gemm() { }
23383 virtual void _check_hemm() { }
23387 TR& dErr,
int transp_mode)
throw(cvmexception)
23392 mA._solve(mB, *
this, dErr,
nullptr,
nullptr, transp_mode);
23398 tint& rank, TR tol)
const throw(cvmexception)
23406 __gelss(mA, mB, mX, tol, sv1, rank);
23408 __gelsd(mA, mB, mX, tol, sv1, rank);
23415 tint& rank, TR tol)
throw(cvmexception)
23424 __gelss(mA2, mB, *
this, tol, sv1, rank);
23426 __gelsd(mA2, mB, *
this, tol, sv1, rank);
23432 tint& rank, TR tol)
const throw(cvmexception)
23441 __gelss(mA, mB, mX, tol, sv1, rank);
23443 __gelsd(mA, mB, mX, tol, sv1, rank);
23446 const TC* pResult = mX;
23459 template<
typename TR,
typename TC>
23684 :
BaseCMatrix(m.msize(), m.nsize(), m.msize(), false)
23716 __copy<TC>(this->
msize(), v, v.
incr(), this->
get(), this->
msize() + 1);
23747 :
BaseCMatrix(m.msize(), m.msize(), m.msize(), true)
23750 __copy2<TR,TC>(this->
get(), this->
size(), this->
incr(), m._pd(),
nullptr);
23753 __copy2<TR,TC>(this->
get(), this->
size(), this->
incr(),
nullptr, m._pd());
23824 :
BaseCMatrix(mRe.msize(), mRe.nsize(), mRe.msize(), false)
23828 __copy2<TR,TC>(this->
get(), this->
size(), this->
incr(), mRe, mIm, mRe.
incr(), mIm.incr());
23862 m._check_submatrix();
23869 return this->_ij_proxy_val(nRow -
CVM0, nCol -
CVM0);
23876 return this->_ij_val(nRow -
CVM0, nCol -
CVM0);
23883 return this->_col(nCol -
CVM0);
23890 return this->_col(nCol -
CVM0);
23897 return this->_row(nRow -
CVM0);
23904 return this->_row(nRow -
CVM0);
23911 __copy<TR>(this->
size(), __get_real_p<TR>(this->
get()), this->
incr() * 2, mRet, mRet.
incr());
23919 __copy<TR>(this->
size(), __get_imag_p<TR>(this->
get()), this->
incr() * 2, mRet, mRet.
incr());
23983 this->_assign(v, v.incr());
23990 this->_assign(pd, 1);
24031 this->_assign_shifted(this->_sub_pointer_nocheck(nRow, nCol), m._pd(), m.msize(), m.nsize(), m.ld());
24070 __copy_real<TR,TC>(this->
get(), this->
size(), this->
incr(), mRe._pd(), mRe.incr());
24102 __copy_imag<TR,TC>(this->
get(), this->
size(), this->
incr(), mIm._pd(), mIm.incr());
24109 this->_set_real_number(d);
24116 this->_set_imag_number(d);
24157 this->_resize2(nNewDim, nNewDim);
24331 this->_msum(m1, m2);
24377 this->_mdiff(m1, m2);
24499 static const TR mone(-1.);
24532 this->_plus_plus();
24563 this->_plus_plus();
24594 this->_minus_minus();
24625 this->_minus_minus();
24656 mRes._scalr(dMult);
24720 mRes._scalc(cMult);
24783 this->_scalr(dMult);
24845 this->_scalc(cMult);
24883 this->_normalize();
24969 return mRes.
conj();
25013 return this->transpose();
25058 return this->conj();
25220 mRes.mult(*
this, m);
25262 this->mult(mTmp, m);
25268 this->_swap_rows(n1, n2);
25274 this->_swap_cols(n1, n2);
25340 return _solve_helper(vB, dErr, 0);
25409 return _solve_helper(vB, dErr, 1);
25478 return _solve_helper(vB, dErr, 2);
25538 static TR dErr(0.);
25539 return this->solve(vB, dErr);
25605 static TR dErr(0.);
25606 return this->solve_tran(vB, dErr);
25672 static TR dErr(0.);
25673 return this->solve_conj(vB, dErr);
25736 return _solve_helper(mB, dErr, 0);
25805 return _solve_helper(mB, dErr, 1);
25874 return _solve_helper(mB, dErr, 2);
25934 static TR dErr(0.);
25935 return this->solve(mB, dErr);
26002 static TR dErr(0.);
26003 return this->solve_tran(mB, dErr);
26070 static TR dErr(0.);
26071 return this->solve_conj(mB, dErr);
26109 return vB / (*this);
26147 return vB % (*this);
26235 const CVector& vB, TR& dErr)
const throw(cvmexception)
26240 this->_solve(vB, vX, dErr, mLU, pPivots, 0);
26328 static TR dErr(0.);
26329 return this->solve_lu(mLU, pPivots, vB, dErr);
26416 const BaseCMatrix& mB, TR& dErr)
const throw(cvmexception)
26421 this->_solve(mB, mX, dErr, mLU, pPivots, 0);
26508 static TR dErr(0.);
26509 return this->solve_lu(mLU, pPivots, mB, dErr);
26541 TC
det()
const throw(cvmexception) {
26542 return this->_det();
26625 this->_low_up(nPivots);
26707 mRes._low_up(nPivots);
26753 TR cond()
const throw(cvmexception)
26756 __cond_num<TR,basic_scmatrix>(*
this, dCondNum);
26804 __inv<basic_scmatrix>(*
this, m);
26852 __inv<basic_scmatrix>(mRes, *
this);
26929 __exp<basic_scmatrix, TR>(*
this, m, tol);
27005 __exp<basic_scmatrix, TR>(mRes, *
this, tol);
27092 if (v.incr() > 1) v1 << v;
27093 __polynom<TC, CVector>(this->
get(), this->
ld(), this->
msize(), m._pd(), m._ldm(), v.
incr() > 1 ? v1 : v);
27179 if (v.
incr() > 1) v1 << v;
27180 __polynom<TC, CVector>(mRes.
get(), mRes.
ld(), this->
msize(), this->
get(), this->
ld(), v.
incr() > 1 ? v1 : v);
27245 this->_eig(vEig, &mEigVect, bRightVect);
27289 this->_eig(vEig,
nullptr,
false);
27340 this->_check_cholesky();
27342 tint nOutInfo = __cholesky<basic_scmatrix>(*this);
27345 this->_clean_low_triangle();
27380 this->_check_bunch_kaufman();
27382 __bunch_kaufman<basic_scmatrix>(*
this, nPivots);
27415 this->_qr_ss(mQ, mR);
27447 this->_lq_ss(mL, mQ);
27479 this->_ql_ss(mQ, mL);
27511 this->_rq_ss(mR, mQ);
27542 this->_plus_plus();
27581 this->_randomize_real(dFrom, dTo);
27587 this->_randomize_imag(dFrom, dTo);
27592 virtual void _eig(CVector& vEig,
basic_scmatrix<TR,TC>* mEigVect,
bool bRightVect)
const throw(cvmexception)
27594 __eig<CVector, basic_scmatrix, basic_scmatrix>(vEig, *
this, mEigVect, bRightVect);
27597 virtual void _solve(
const CVector& vB, CVector& vX,
27598 TR& dErr,
const TC* pLU,
const tint* pPivots,
int transp_mode)
const throw(cvmexception)
27603 if (vB.
incr() > 1) vB1 << vB;
27604 if (vX.incr() > 1) vX1 << vX;
27605 __solve<TR,TC, basic_scmatrix>(*
this, 1, vB.
incr() > 1 ? vB1 : vB, vB.
size(), vX.incr() > 1 ? vX1 : vX, vX.size(),
27606 dErr, pLU, pPivots, transp_mode);
27607 if (vX.incr() > 1) vX = vX1;
27610 virtual void _solve(
const BaseCMatrix& mB, BaseCMatrix& mX,
27611 TR& dErr,
const TC* pLU,
const tint* pPivots,
int transp_mode)
const throw(cvmexception)
27614 __solve<TR,TC, basic_scmatrix>(*
this, mB.
nsize(), mB, mB.
ld(), mX, mX.ld(), dErr, pLU, pPivots, transp_mode);
27617 virtual const TC* _pv()
const override {
return this->
get(); }
27618 virtual TC* _pv()
override {
return this->
get(); }
27623 : BaseCMatrix(nDim, nDim, nLD, bZeroMemory)
27628 : BaseCMatrix(size, incr, msize, nsize, ld)
27634 : BaseCMatrix(pd, nDim, nDim, nLD, nSize)
27640 : BaseCMatrix(pd, nDim, nDim, nLD, nSize)
27644 tint _size()
const override {
return this->
size();}
27645 tint _msize()
const override {
return this->
msize();}
27646 tint _nsize()
const override {
return this->
nsize();}
27647 tint _ld()
const override {
return this->
ld();}
27651 CVector _diag(
tint nDiag)
throw(cvmexception)
override
27653 const tint nD = _abs(nDiag);
27655 return CVector(this->
get() + (nDiag > 0 ? nDiag * this->
ld() : nD), this->
msize() - nD, this->
ld() + 1);
27660 const CVector _diag(
tint nDiag)
const throw(cvmexception)
override
27662 const tint nD = _abs(nDiag);
27664 return CVector(this->
get() + (nDiag > 0 ? nDiag * this->
ld() : nD), this->
msize() - nD, this->
ld() + 1);
27675 virtual void _transp() {
27676 this->_sq_transp();
27680 virtual void _conj()
27683 __conj<TC>(this->
get(), this->
size(), this->
incr());
27686 virtual void _plus_plus() {
27687 this->_sq_plus_plus();
27690 virtual void _minus_minus() {
27691 this->_sq_minus_minus();
27694 virtual TC _det()
const throw(cvmexception)
27697 switch (this->
msize()) {
27701 cDet = this->_ij_val(0, 0);
27704 cDet = this->_ij_val(0, 0) * this->_ij_val(1, 1) -
27705 this->_ij_val(1, 0) * this->_ij_val(0, 1);
27709 static const TC one(1., 0.);
27711 CVector vUpDiag = this->_low_up_diag(naPivots);
27715 cDet *= vUpDiag[i];
27716 if (i + (1 -
CVM0) != naPivots[i]) {
27721 catch (
const cvmexception& e) {
27729 virtual void _low_up(
tint* nPivots)
throw(cvmexception) {
27730 __low_up<basic_scmatrix>(*
this, nPivots);
27765 virtual void _check_cholesky() { }
27766 virtual void _check_bunch_kaufman() { }
27769 CVector _solve_helper(
const CVector& vB, TR& dErr,
int transp_mode)
const throw(cvmexception)
27772 CVector vX(this->
msize());
27773 this->_solve(vB, vX, dErr,
nullptr,
nullptr, transp_mode);
27777 BaseCMatrix _solve_helper(
const BaseCMatrix& mB, TR& dErr,
int transp_mode)
const throw(cvmexception)
27781 this->_solve(mB, mX, dErr,
nullptr,
nullptr, transp_mode);
27810 template<
typename TR,
typename TC>
27837 virtual tint _size()
const = 0;
27838 virtual tint _msize()
const = 0;
27839 virtual tint _nsize()
const = 0;
27840 virtual tint _ld()
const = 0;
27841 virtual void _set_p(TC* pf) = 0;
27845 virtual const TC* _pb()
const = 0;
27846 virtual TC* _pb() = 0;
27848 #ifdef CVM_USE_POOL_MANAGER
27849 void _bresize(
tint nNewM)
throw(cvmexception)
27851 void _bresize(std::shared_ptr<TC>&
mp,
tint nNewM)
throw (cvmexception)
27854 const tint mm = this->_msize();
27855 const tint nSize = this->_size();
27858 const tint nNewSize = nNewM * (1 + this->lsize() + this->usize());
27859 TC* pd = cvmMalloc<TC>(nNewSize);
27860 if (nNewSize > nSize) cvmZeroMemory<TC>(pd, nNewSize);
27861 const tint nMinSize = _cvm_min<tint>(nSize, nNewSize);
27862 if (nMinSize > 0) {
27863 __copy<TC>(nMinSize, this->_pb(), 1, pd, 1);
27866 #ifdef CVM_USE_POOL_MANAGER
27867 TC* pB = this->_pb();
27869 this->_set(pd, nNewSize, nNewM, nNewM, 1, nNewSize > 0 ? 1 + this->lsize() + this->usize() : 0);
27872 this->_set(
nullptr, nNewSize, nNewM, nNewM, 1, nNewSize > 0 ? 1 + this->lsize() + this->usize() : 0);
27877 void _check_dimensions()
const throw(cvmexception)
27883 void _check_dimensions(
const BandMatrix& m)
const throw(cvmexception) {
27889 bool _bcontinuous()
const {
27890 return this->_msize() == 0 || this->_ld() == 1 + this->lsize() + this->usize();
27897 if (pd != m.
get()) {
27898 __copy<TC>(this->_size(), m, m.
incr(), pd, 1);
27903 void _b_for_each(TC d,
bool bRandom,
bool bReal, TR dFrom, TR dTo)
27905 tint i, jc, ju, jl;
27906 const tint nCol = 1 + this->lsize() + this->usize();
27908 const tint mn = this->_nsize();
27909 for (
tint j = 0; j <
mn; ++j) {
27911 ju = this->usize() - j;
27913 for (i = 0; i < nCol; ++i) {
27914 if (i >= ju && i < jl) {
27917 TR* pd = bReal ? __get_real_p<TR>(pB + (jc + i)) : __get_imag_p<TR>(pB + (jc + i));
27930 static const TR zero = TR(0.);
27931 this->_b_for_each(d,
false,
false, zero, zero);
27935 void _b_randomize(
bool bReal, TR dFrom, TR dTo)
27937 static const TC zero(0.);
27938 this->_b_for_each(zero,
true, bReal, dFrom, dTo);
27949 const tint mn = this->_nsize();
27952 for (j = 0; j <
mn; ++j) {
27954 this->_get_col(j, rv, 1, &nLen, &nShift);
27957 for (i = nShift + 1; i <= nLen; ++i) {
27958 rSum += _abs(rv[i - (1 -
CVM0)]);
27961 if (rSum > rNorm) {
27969 TR _bnorminf()
const
27976 const tint mm = this->_msize();
27979 for (i = 0; i <
mm; ++i) {
27981 this->_get_row(i, rv, 1, &nLen, &nShift);
27984 for (j = nShift + 1; j <= nLen; ++j) {
27985 rSum += _abs(rv[j - (1 -
CVM0)]);
27988 if (rSum > rNorm) {
27996 this->_bset(TC(0.));
28002 static const TC zero = TC(0.);
28004 const tint nA = j - this->usize();
28005 CVM_ASSERT(pd, (i + j * (1 + this->lsize() + this->usize()) - nA + 1) *
sizeof(TC))
28006 return(nA < 0 || i >= nA) && i <= this->lsize() + j ? type_proxy<TC, TR>(pd[i + j * (1 + this->lsize() + this->usize()) - nA], false) :
28007 type_proxy<TC, TR>(zero, true);
28011 TC _b_ij_val(
tint i,
tint j)
const
28013 static const TC zero = TC(0.);
28014 const TC* pd = _pb();
28015 const tint nA = j - this->usize();
28016 CVM_ASSERT(pd, (i + j * (1 + this->lsize() + this->usize()) - nA + 1) *
sizeof(TC))
28017 return(nA < 0 || i >= nA) && i <= this->lsize() + j ? pd[i + j * (1 + this->lsize() + this->usize()) - nA] : zero;
28020 void _get_col(
tint i, TC* pCol,
tint nIncr,
tint* pnLen = nullptr,
tint* pnShift = nullptr)
const
28022 const TC* pd = _pb();
28023 const tint mm = this->_msize();
28024 const tint mn = this->_nsize();
28025 const tint nCol = 1 + this->lsize() + this->usize();
28027 tint nShiftSrc = 0;
28028 tint nShiftDest = 0;
28030 if (i > this->usize()) {
28031 nShiftDest = i - this->usize();
28033 nShiftSrc = this->usize() - i;
28037 if (mm - i <= this->lsize()) {
28038 nS -= this->lsize() + 1 - (mn - i);
28042 pd + i * nCol + nShiftSrc,
28047 if (pnLen) *pnLen = nS;
28048 if (pnShift) *pnShift = nShiftDest;
28051 void _get_row(
tint i, TC* pCol,
tint nIncr,
tint* pnLen =
nullptr,
tint* pnShift =
nullptr)
const
28053 const TC* pd = _pb();
28054 const tint mm = this->_msize();
28055 const tint mn = this->_nsize();
28056 const tint nCol = this->lsize() + this->usize();
28058 tint nShiftSrc = i + this->usize();
28059 tint nShiftDest = 0;
28061 if (i > this->lsize()) {
28062 nShiftDest = i - this->lsize();
28063 nShiftSrc += nShiftDest * nCol;
28066 if (mn - i > this->usize()) {
28067 nS -= (mm - i) - this->usize() - 1;
28076 if (pnLen) *pnLen = nS;
28077 if (pnShift) *pnShift = nShiftDest;
28081 #ifdef CVM_USE_POOL_MANAGER
28082 void _btransp() throw(cvmexception)
28084 void _btransp(std::shared_ptr<TC>&
mp)
throw (cvmexception)
28088 const tint mn = this->_nsize();
28089 if (this->lsize() > 0 || this->usize() > 0) {
28090 const tint nLU = this->lsize() + this->usize();
28091 const tint nCol = 1 + nLU;
28092 const tint nEls = nCol *
mn;
28093 tint nS, nShiftSrc;
28096 TC* pd = cvmMalloc<TC>(nEls);
28099 for (
tint i = 0; i <
mn; ++i) {
28103 if (i < this->usize()) {
28104 nShiftSrc = this->usize() - i;
28107 if (mn - i <= this->lsize()) {
28108 nS -= this->lsize() + 1 - (mn - i);
28111 pL = mpd + i * nCol + nShiftSrc;
28112 pR = pd + (i > this->usize() ? (i - this->usize() + 1) * nCol - 1 : this->lsize() + i);
28114 __copy<TC>(nS, pL, 1, pR, nLU);
28117 #ifdef CVM_USE_POOL_MANAGER
28122 this->_set_p(
nullptr);
28127 void _b_plus_plus()
28130 static const TC one(1.);
28131 const tint mn = this->_nsize();
28132 const tint nNext = 1 + this->lsize() + this->usize();
28133 const tint nSize = nNext *
mn;
28135 for (
tint i = this->usize(); i < nSize; i += nNext) {
28140 void _b_minus_minus()
28143 static const TC one(1.);
28144 const tint mn = this->_nsize();
28145 const tint nNext = 1 + this->lsize() + this->usize();
28146 const tint nSize = nNext *
mn;
28148 for (
tint i = this->usize(); i < nSize; i += nNext) {
28154 #ifdef CVM_USE_POOL_MANAGER
28155 void _b_replace(
const BandMatrix& m)
throw(cvmexception)
28157 void _b_replace(std::shared_ptr<TC>&
mp,
const BandMatrix& m)
throw(cvmexception)
28160 #ifdef CVM_USE_POOL_MANAGER
28164 tint nSize = m._size();
28165 TC* pd = cvmMalloc<TC>(nSize);
28169 #ifdef CVM_USE_POOL_MANAGER
28170 this->_set(pd, nSize, m._msize(), m._nsize(), 1, m._ld());
28173 this->_set(
nullptr, nSize, m._msize(), m._nsize(), 1, m._ld());
28177 #ifdef CVM_USE_POOL_MANAGER
28178 void _resize_lu(
tint nNewKL,
tint nNewKU)
throw(cvmexception)
28180 void _resize_lu(std::shared_ptr<TC>&
mp,
tint nNewKL,
tint nNewKU)
throw(cvmexception)
28183 if (nNewKL != this->lsize() || nNewKU != this->usize()) {
28184 if (this->lsize() < 0)
throw cvmexception(
CVM_WRONGSIZE, this->lsize());
28185 if (this->usize() < 0)
throw cvmexception(
CVM_WRONGSIZE, this->usize());
28187 const tint mm = this->_msize();
28188 const tint mn = this->_nsize();
28189 const tint nOldLD = 1 + this->lsize() + this->usize();
28190 const tint nNewLD = 1 + nNewKL + nNewKU;
28191 const tint nMinKL = _cvm_min<tint>(this->lsize(), nNewKL);
28192 const tint nMinKU = _cvm_min<tint>(this->usize(), nNewKU);
28193 const tint nNewSize = mn * (1 + nNewKL + nNewKU);
28194 TC* pd = cvmMalloc<TC>(nNewSize);
28197 for (
tint i = -nMinKU; i <= nMinKL; ++i) {
28198 __copy<TC>(
mn, pB + (this->usize() + i), nOldLD, pd + (nNewKU + i), nNewLD);
28202 #ifdef CVM_USE_POOL_MANAGER
28204 this->_set(pd, nNewSize, mm, mn, 1, nNewLD);
28207 this->_set(
nullptr, nNewSize, mm, mn, 1, nNewLD);
28269 const tint* _pl()
const {
28273 const tint* _pu()
const {
28308 template<
typename TR>
28311 typedef std::complex<TR> TC;
28419 this->_check_dimensions();
28465 :
BaseSRMatrix(pd, nDim, 1 + nKL + nKU, nDim * (1 + nKL + nKU)),
28468 this->_check_dimensions();
28514 :
BaseSRMatrix(pd, nDim, 1 + nKL + nKU, nDim * (1 + nKL + nKU)),
28517 this->_check_dimensions();
28549 :
BaseSRMatrix(m.msize(), 1 + m.lsize() + m.usize(), false),
28576 this->_mmove(std::move(m));
28615 this->_check_dimensions();
28616 _copy_b_matrix<TR,TR, BaseRMatrix, basic_srbmatrix>(
const_cast<BaseRMatrix&
>(m), *
this,
true);
28647 __copy<TR>(this->
msize(), v, v.
incr(), this->
get(), 1);
28654 return this->_ij_proxy_val(nRow -
CVM0, nCol -
CVM0);
28661 return this->_ij_val(nRow -
CVM0, nCol -
CVM0);
28668 return this->_col(nCol -
CVM0);
28675 return this->_row(nRow -
CVM0);
28712 this->_check_dimensions(m);
28732 this->_check_dimensions(m);
28741 this->_assign(v, v.incr());
28748 this->_assign(pd, 1);
28761 this->_resize2(nNewDim, nNewDim);
28805 #ifdef CVM_USE_POOL_MANAGER
28806 this->_resize_lu(nNewKL, nNewKU);
28808 this->_resize_lu(this->
mp, nNewKL, nNewKU);
28925 #ifdef CVM_USE_POOL_MANAGER
28926 this->_b_replace(m);
28928 this->_b_replace(this->
mp, m);
28975 this->_check_dimensions(m);
29022 this->_check_dimensions(m);
29076 this->_check_dimensions(m1);
29077 this->_check_dimensions(m2);
29078 this->_msum(m1, m2);
29130 this->_check_dimensions(m1);
29131 this->_check_dimensions(m2);
29132 this->_mdiff(m1, m2);
29181 this->_check_dimensions(m);
29231 this->_check_dimensions(m);
29264 static const TR mone(-1.);
29299 this->_plus_plus();
29332 this->_plus_plus();
29365 this->_minus_minus();
29398 this->_minus_minus();
29405 mRes._scalr(dMult);
29418 this->_scalr(dMult);
29430 this->_normalize();
29497 this->_multiply(vRes, v,
false);
29584 _cvm_min<tint>(this->msize() - 1, this->
lsize() + m.lsize()),
29585 _cvm_min<tint>(this->
msize() - 1, this->
usize() + m.usize()));
29591 return v % (*this);
29689 this->_low_up(nPivots);
29696 mRes._low_up(nPivots);
29703 this->_plus_plus();
29715 this->_b_randomize(
true, dFrom, dTo);
29722 TR dNorm = TR(0.), d;
29725 for (i = 0; i <= this->
lsize(); ++i) {
29730 for (i = 1; i <= this->
usize(); ++i) {
29736 return _sqrt(dNorm);
29741 return this->_bnorm1();
29746 return this->_bnorminf();
29754 this->_mult(m1, m2);
29772 const TR* _pd()
const override
29779 void _svd(RVector& vRes, BaseSRMatrix* pmU, BaseSRMatrix* pmVH)
const throw(cvmexception)
override
29781 if (pmU !=
nullptr && pmVH !=
nullptr) {
29785 __svd<TR,basic_srbmatrix, BaseSRMatrix>(vRes, vRes.size(), vRes.incr(), *
this, pmU, pmVH);
29788 void _pinv(BaseRMatrix& mX, TR threshold)
const throw(cvmexception)
override {
29789 __pinv<TR,basic_srbmatrix, BaseRMatrix>(mX, *
this, threshold);
29792 void _eig(CVector& vEig,
basic_scmatrix<TR,TC>* mEigVect,
bool bRightVect)
const throw(cvmexception)
override
29795 mSM._eig(vEig, mEigVect, bRightVect);
29798 void _solve(
const RVector& vB, RVector& vX,
29799 TR& dErr,
const TR* pLU,
const tint* pPivots,
int transp_mode)
const throw(cvmexception)
override
29804 if (vB.
incr() > 1) vB1 << vB;
29805 if (vX.incr() > 1) vX1 << vX;
29806 __solve<TR,TR, basic_srbmatrix>(*
this, 1, vB.
incr() > 1 ? vB1 : vB, vB.
size(),
29807 vX.incr() > 1 ? vX1 : vX, vX.size(), dErr, pLU, pPivots, transp_mode);
29808 if (vX.incr() > 1) vX = vX1;
29811 void _solve(
const BaseRMatrix& mB, BaseRMatrix& mX, TR& dErr,
29812 const TR* pLU,
const tint* pPivots,
int transp_mode)
const throw(cvmexception)
override
29815 __solve<TR,TR, basic_srbmatrix>(*
this, mB.
nsize(), mB, mB.
ld(), mX, mX.ld(), dErr, pLU, pPivots, transp_mode);
29821 void _gbmv(
bool bLeft, TR dAlpha,
const RVector& v, TR dBeta, RVector& vRes)
const
29826 if (vRes.get() == pDv) vTmp << v;
29827 if (vRes.get() == this->
get()) mTmp << *
this;
29828 __gbmv<TR,basic_srbmatrix, RVector>(bLeft, vRes.get() == this->
get() ? mTmp : *
this, dAlpha,
29829 vRes.
get() == pDv ? vTmp : v, dBeta, vRes);
29832 void _check_submatrix()
const throw(cvmexception)
override {
29836 const TR* _pb()
const override {
29837 return this->
get();
29840 TR* _pb()
override {
29841 return this->
get();
29845 void _bake_SM()
const
29849 _copy_b_matrix<TR,TR, BaseSRMatrix, basic_srbmatrix>(
mSM, *
const_cast<basic_srbmatrix*
>(
this),
false);
29852 tint _size()
const override {
return this->
size(); }
29853 tint _msize()
const override {
return this->
msize(); }
29854 tint _nsize()
const override {
return this->
nsize(); }
29855 tint _ld()
const override {
return this->
ld(); }
29857 void _set_p(TR* pf)
override
29859 #ifdef CVM_USE_POOL_MANAGER
29868 #ifdef CVM_USE_POOL_MANAGER
29876 this->mincr = nIncr;
29881 const TR* _pp(
const BaseMatrix& m)
const override {
29885 tint _ldm()
const override {
29889 const tint* _pldm()
const override {
29894 RVector _row(
tint m)
override
29896 RVector vRet(this->
msize());
29897 this->_get_row(m, vRet, vRet.incr());
29903 const RVector _row(
tint m)
const override
29905 RVector vRet(this->
msize());
29906 this->_get_row(m, vRet, vRet.incr());
29911 RVector _col(
tint n)
override
29913 RVector vRet(this->
msize());
29914 this->_get_col(n, vRet, vRet.incr());
29920 const RVector _col(
tint n)
const override
29922 RVector vRet(this->
msize());
29923 this->_get_col(n, vRet, vRet.incr());
29927 RVector _diag(
tint nDiag)
throw(cvmexception)
override
29929 const tint nD = _abs(nDiag);
29937 return RVector(this->
get() + this->
usize() + (nDiag < 0 ? nD : nD * nLU), this->
msize() - nD, nLU + 1);
29940 const RVector _diag(
tint nDiag)
const throw(cvmexception)
override
29942 const tint nD = _abs(nDiag);
29950 return RVector(this->
get() + this->
usize() + (nDiag < 0 ? nD : nD * nLU), this->
msize() - nD, nLU + 1);
29955 void _check_ger() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"ger",
"srbmatrix");}
29956 void _check_rank1update() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"rank1update",
"srbmatrix");}
29957 void _check_gemm() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"gemm",
"srbmatrix");}
29958 void _check_symm() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"symm",
"srbmatrix");}
29959 void _check_cholesky() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"cholesky",
"srbmatrix");}
29960 void _check_bunch_kaufman() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"bunch_kaufman",
"srbmatrix");}
29962 void _resize2(
tint nNewM,
tint) throw(cvmexception)
override
29964 #ifdef CVM_USE_POOL_MANAGER
29965 this->_bresize(nNewM);
29967 this->_bresize(this->
mp, nNewM);
29971 bool _continuous()
const override {
29972 return this->_bcontinuous();
29975 void _massign(
const BaseMatrix& m)
override {
29976 this->_mbassign(m);
29979 void _set(TR d)
override {
29983 void _vanish()
override {
29989 return this->_b_ij_proxy_val(i, j);
29993 TR _ij_val(
tint i,
tint j)
const override {
29994 return this->_b_ij_val(i, j);
29997 void _transp() throw(cvmexception)
override
29999 #ifdef CVM_USE_POOL_MANAGER
30002 this->_btransp(this->
mp);
30006 void _plus_plus()
override {
30007 this->_b_plus_plus();
30010 void _minus_minus()
override {
30011 this->_b_minus_minus();
30014 tint _indofmax()
const override
30021 tint _indofmin()
const override
30036 void _scalr(TR d)
override {
30037 __scal<TR,TR>(this->
get(), this->
size(), this->
incr(), d);
30040 void _mult(
const BaseRMatrix& m1,
const BaseRMatrix& m2)
throw(cvmexception)
override
30045 BaseSRMatrix mR(this->
msize());
30047 #ifdef CVM_USE_POOL_MANAGER
30048 this->_resize_lu(this->
msize() - 1, this->
msize() - 1);
30050 this->_resize_lu(this->
mp, this->
msize() - 1, this->
msize() - 1);
30052 _copy_b_matrix<TR,TR, BaseSRMatrix, basic_srbmatrix>(
const_cast<BaseSRMatrix&
>(mR), *
this,
true);
30055 void _multiply(RVector& vRes,
const RVector& v,
bool bLeft)
const override
30057 static const TR zero(0.);
30058 static const TR one(1.);
30059 this->_gbmv(bLeft, one, v, zero, vRes);
30062 void _randomize(TR dFrom, TR dTo)
override {
30063 __randomize<TR>(this->
get(), this->
size(), this->
incr(), dFrom, dTo);
30066 void _low_up(
tint* nPivots)
throw(cvmexception)
override {
30067 __low_up<basic_srbmatrix>(*
this, nPivots);
30070 tint _ld_for_replace()
const override {
30074 tint _size_for_replace()
const override {
30075 return this->mm * this->
mn;
30109 template<
typename TR,
typename TC>
30213 this->_check_dimensions();
30260 :
BaseSCMatrix(pd, nDim, 1 + nKL + nKU, nDim * (1 + nKL + nKU)),
30263 this->_check_dimensions();
30310 :
BaseSCMatrix(pd, nDim, 1 + nKL + nKU, nDim * (1 + nKL + nKU)),
30313 this->_check_dimensions();
30346 :
BaseSCMatrix(m.msize(), 1 + m.lsize() + m.usize(), false),
30373 this->_mmove(std::move(m));
30411 this->_check_dimensions();
30412 _copy_b_matrix<TR,TC, BaseCMatrix, basic_scbmatrix>(
const_cast<BaseCMatrix&
>(m), *
this,
true);
30441 __copy<TC>(this->
msize(), v, v.
incr(), this->
get(), 1);
30479 __copy2<TR,TC>(this->
get(), this->
size(), this->
incr(), m.
get(),
nullptr);
30482 __copy2<TR,TC>(this->
get(), this->
size(), this->
incr(),
nullptr, m.
get());
30522 __copy2<TR,TC>(this->
get(), this->
size(), this->
incr(), mRe, mIm, mRe.
incr(), mIm.incr());
30529 return this->_ij_proxy_val(nRow -
CVM0, nCol -
CVM0);
30536 return this->_ij_val(nRow -
CVM0, nCol -
CVM0);
30543 return this->_col(nCol -
CVM0);
30550 return this->_row(nRow -
CVM0);
30557 __copy<TR>(this->
size(), __get_real_p<TR>(this->
get()), this->
incr() * 2, mRet, mRet.
incr());
30565 __copy<TR>(this->
size(), __get_imag_p<TR>(this->
get()), this->
incr() * 2, mRet, mRet.
incr());
30603 this->_check_dimensions(m);
30623 this->_check_dimensions(m);
30632 this->_assign(v, v.incr());
30639 this->_assign(pd, 1);
30680 __copy_real<TR,TC>(this->
get(), this->
size(), this->
incr(), mRe, mRe.incr());
30715 __copy_imag<TR,TC>(this->
get(), this->
size(), this->
incr(), mIm, mIm.incr());
30722 this->_set_real_number(d);
30729 this->_set_imag_number(d);
30735 this->_resize2(nNewDim, nNewDim);
30777 #ifdef CVM_USE_POOL_MANAGER
30778 this->_resize_lu(nNewKL, nNewKU);
30780 this->_resize_lu(this->
mp, nNewKL, nNewKU);
30819 this->lsize() == m.
lsize() && this->usize() == m.
usize() && this->_mequals(m);
30854 return !(this->operator == (m));
30899 #ifdef CVM_USE_POOL_MANAGER
30900 this->_b_replace(m);
30902 this->_b_replace(this->
mp, m);
30956 this->_check_dimensions(m);
31010 this->_check_dimensions(m);
31060 this->_check_dimensions(m1);
31061 this->_check_dimensions(m2);
31062 this->_msum(m1, m2);
31110 this->_check_dimensions(m1);
31111 this->_check_dimensions(m2);
31112 this->_mdiff(m1, m2);
31161 this->_check_dimensions(m);
31211 this->_check_dimensions(m);
31218 static const TR mone(-1.);
31227 this->_plus_plus();
31234 this->_plus_plus();
31241 this->_minus_minus();
31248 this->_minus_minus();
31255 mRes._scalr(dMult);
31269 mRes._scalc(cMult);
31282 this->_scalr(dMult);
31294 this->_scalc(cMult);
31306 this->_normalize();
31321 return mRes.
conj();
31446 this->_multiply(vRes, v,
false);
31532 _cvm_min<tint>(this->msize() - 1, this->lsize() + m.lsize()),
31533 _cvm_min<tint>(this->
msize() - 1, this->usize() + m.usize()));
31539 return v % (*this);
31636 this->_low_up(nPivots);
31643 mRes._low_up(nPivots);
31650 this->_plus_plus();
31662 this->_b_randomize(
true, dFrom, dTo);
31668 this->_b_randomize(
false, dFrom, dTo);
31675 TR dNorm = TR(0.), d;
31677 for (i = 0; i <= this->lsize(); ++i) {
31682 for (i = 1; i <= this->usize(); ++i) {
31687 return _sqrt(dNorm);
31691 return this->_bnorm1();
31695 return this->_bnorminf();
31702 this->_mult(m1, m2);
31719 const TC* _pd()
const override
31726 void _svd(
basic_rvector<TR>& vRes, BaseSCMatrix* pmU, BaseSCMatrix* pmVH)
const throw(cvmexception)
override
31728 if (pmU !=
nullptr && pmVH !=
nullptr) {
31732 __svd<TR,basic_scbmatrix, BaseSCMatrix>(vRes, vRes.size(), vRes.incr(), *
this, pmU, pmVH);
31735 void _pinv(BaseCMatrix& mX, TR threshold)
const throw(cvmexception)
override {
31736 __pinv<TR,basic_scbmatrix, BaseCMatrix>(mX, *
this, threshold);
31739 void _eig(CVector& vEig,
basic_scmatrix<TR,TC>* mEigVect,
bool bRightVect)
const throw(cvmexception)
override
31742 mSM._eig(vEig, mEigVect, bRightVect);
31745 void _solve(
const CVector& vB, CVector& vX,
31746 TR& dErr,
const TC* pLU,
const tint* pPivots,
int transp_mode)
const throw(cvmexception)
override
31751 if (vB.
incr() > 1) vB1 << vB;
31752 if (vX.incr() > 1) vX1 << vX;
31753 __solve<TR,TC, basic_scbmatrix>(*
this, 1, vB.
incr() > 1 ? vB1 : vB, vB.
size(), vX.incr() > 1 ? vX1 : vX, vX.size(),
31754 dErr, pLU, pPivots, transp_mode);
31755 if (vX.incr() > 1) vX = vX1;
31758 void _solve(
const BaseCMatrix& mB, BaseCMatrix& mX,
31759 TR& dErr,
const TC* pLU,
const tint* pPivots,
int transp_mode)
const throw(cvmexception)
override
31762 __solve<TR,TC, basic_scbmatrix>(*
this, mB.
nsize(), mB, mB.
ld(), mX, mX.ld(), dErr, pLU, pPivots, transp_mode);
31768 void _gbmv(
bool bLeft, TC dAlpha,
const CVector& v, TC dBeta, CVector& vRes)
const
31773 if (vRes.get() == pDv) vTmp << v;
31774 if (vRes.get() == this->
get()) mTmp << *
this;
31775 __gbmv<TC, basic_scbmatrix, CVector>(bLeft, vRes.get() == this->
get() ? mTmp : *
this, dAlpha,
31776 vRes.
get() == pDv ? vTmp : v, dBeta, vRes);
31779 void _check_submatrix()
const throw(cvmexception)
override {
31783 const TC* _pb()
const override {
31784 return this->
get();
31787 TC* _pb()
override {
31788 return this->
get();
31792 void _bake_SM()
const
31794 mSM.resize(this->
msize());
31796 _copy_b_matrix<TR,TC, BaseSCMatrix, basic_scbmatrix>(mSM, *
const_cast<basic_scbmatrix*
>(
this),
false);
31799 tint _size()
const override {
return this->
size(); }
31800 tint _msize()
const override {
return this->
msize(); }
31801 tint _nsize()
const override {
return this->
nsize(); }
31802 tint _ld()
const override {
return this->
ld(); }
31804 void _set_p(TC* pf)
override
31806 #ifdef CVM_USE_POOL_MANAGER
31815 #ifdef CVM_USE_POOL_MANAGER
31823 this->mincr = nIncr;
31827 const TC* _pp(
const BaseMatrix& m)
const override {
31831 tint _ldm()
const override {
31835 const tint* _pldm()
const override {
31840 CVector _row(
tint m)
override
31842 CVector vRet(this->
msize());
31843 this->_get_row(m, vRet, vRet.incr());
31849 const CVector _row(
tint m)
const override
31851 CVector vRet(this->
msize());
31852 this->_get_row(m, vRet, vRet.incr());
31857 CVector _col(
tint n)
override
31859 CVector vRet(this->
msize());
31860 this->_get_col(n, vRet, vRet.incr());
31866 const CVector _col(
tint n)
const override
31868 CVector vRet(this->
msize());
31869 this->_get_col(n, vRet, vRet.incr());
31873 CVector _diag(
tint nDiag)
throw(cvmexception)
override
31875 const tint nD = _abs(nDiag);
31882 const tint nLU = this->lsize() + this->usize();
31883 return CVector(this->
get() + this->usize() + (nDiag < 0 ? nD : nD * nLU), this->
msize() - nD, nLU + 1);
31886 const CVector _diag(
tint nDiag)
const throw(cvmexception)
override
31888 const tint nD = _abs(nDiag);
31895 const tint nLU = this->lsize() + this->usize();
31896 return CVector(this->
get() + this->usize() + (nDiag < 0 ? nD : nD * nLU), this->
msize() - nD, nLU + 1);
31899 void _resize2(
tint nNewM,
tint) throw(cvmexception)
override
31901 #ifdef CVM_USE_POOL_MANAGER
31902 this->_bresize(nNewM);
31904 this->_bresize(this->
mp, nNewM);
31908 bool _continuous()
const override
31910 return this->_bcontinuous();
31913 void _massign(
const BaseMatrix& m)
override {
31914 this->_mbassign(m);
31917 void _set(TC d)
override {
31921 void _vanish()
override {
31927 return this->_b_ij_proxy_val(i, j);
31931 TC _ij_val(
tint i,
tint j)
const override {
31932 return this->_b_ij_val(i, j);
31935 void _transp() throw(cvmexception)
override
31937 #ifdef CVM_USE_POOL_MANAGER
31940 this->_btransp(this->
mp);
31944 void _plus_plus()
override {
31945 this->_b_plus_plus();
31948 void _minus_minus()
override {
31949 this->_b_minus_minus();
31952 tint _indofmax()
const override
31956 return mSM.indofmax();
31959 tint _indofmin()
const override
31963 return mSM.indofmin();
31973 void _scalr(TR d)
override {
31974 __scal<TR,TC>(this->
get(), this->
size(), this->
incr(), d);
31977 void _scalc(TC d)
override {
31978 __scal<TC, TC>(this->
get(), this->
size(), this->
incr(), d);
31981 void _mult(
const BaseCMatrix& m1,
const BaseCMatrix& m2)
throw(cvmexception)
override
31986 BaseSCMatrix mR(this->
msize());
31988 #ifdef CVM_USE_POOL_MANAGER
31989 this->_resize_lu(this->
msize() - 1, this->
msize() - 1);
31991 this->_resize_lu(this->
mp, this->
msize() - 1, this->
msize() - 1);
31993 _copy_b_matrix<TR,TC, BaseSCMatrix, basic_scbmatrix>(
const_cast<BaseSCMatrix&
>(mR), *
this,
true);
31996 void _multiply(CVector& vRes,
const CVector& v,
bool bLeft)
const override
31998 static const TC zero(0., 0.);
31999 static const TC one(1., 0.);
32000 this->_gbmv(bLeft, one, v, zero, vRes);
32003 void _low_up(
tint* nPivots)
throw(cvmexception)
override {
32004 __low_up<basic_scbmatrix>(*
this, nPivots);
32007 tint _ld_for_replace()
const override {
32011 tint _size_for_replace()
const override {
32012 return this->mm * this->
mn;
32017 void _check_geru() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"geru",
"scbmatrix");}
32018 void _check_gerc() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"gerc",
"scbmatrix");}
32019 void _check_rank1update_u() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"rank1update_u",
"scbmatrix");}
32020 void _check_rank1update_c() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"rank1update_c",
"scbmatrix");}
32021 void _check_gemm() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"gemm",
"scbmatrix");}
32022 void _check_hemm() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"hemm",
"scbmatrix");}
32023 void _check_cholesky() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"cholesky",
"scbmatrix");}
32024 void _check_bunch_kaufman() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"bunch_kaufman",
"scbmatrix");}
32035 template<
typename TR>
32141 this->_check_symmetric(tol);
32180 this->_check_symmetric(tol);
32267 this->_check_symmetric(tol);
32326 :
BaseSRMatrix(m._sub_pointer(nRowCol, nRowCol, nDim, nDim), nDim, m.
ld(), nDim * nDim)
32334 return this->_ij_val(nRow -
CVM0, nCol -
CVM0);
32341 return this->_col(nCol -
CVM0);
32348 return this->_row(nRow -
CVM0);
32354 return this->_diag(nDiag);
32439 this->_assign(v, v.
incr());
32440 this->_check_symmetric(tol);
32475 this->_assign(pd, 1);
32476 this->_check_symmetric(tol);
32514 this->_assign_shifted(this->_sub_pointer_nocheck(nRowCol, nRowCol), m._pd(), m.
msize(), m.
nsize(), m.
ld());
32552 this->_set_at(nRow -
CVM0, nCol -
CVM0, d);
32588 this->_diag(nDiag) = vDiag;
32590 this->_diag(-nDiag) = vDiag;
32597 this->_resize2(nNewDim, nNewDim);
32660 return !(this->operator == (m));
32815 this->_msum(m1, m2);
32855 this->_mdiff(m1, m2);
32968 static const TR mone(-1.);
32977 this->_plus_plus();
32984 this->_plus_plus();
32991 this->_minus_minus();
32998 this->_minus_minus();
33005 mRes._scalr(dMult);
33018 this->_scalr(dMult);
33030 this->_normalize();
33061 this->_multiply(vRes, v,
false);
33069 static const TR zero = TR(0.);
33070 static const TR one = TR(1.);
33072 mRes._symm(
true, *
this, m, one, zero);
33116 static const TR zero = TR(0.);
33117 static const TR one = TR(1.);
33119 mRes._symm(
true, *
this, m, one, zero);
33126 return v % (*this);
33180 __syrk<TR,basic_srsmatrix>(
false, alpha, 1, vc, vc.size(), beta, *
this);
33249 __syrk<TR,basic_srsmatrix>(bTransp, alpha, bTransp ? m.
msize() : m.
nsize(), m, m.
ld(), beta, *
this);
33305 __syr2k<TR,basic_srsmatrix>(
false, alpha, 1, v1c, v1c.size(), v2c, v2c.
size(), beta, *
this);
33387 __syr2k<TR,basic_srsmatrix>(bTransp, alpha, bTransp ? m1.msize() : m1.nsize(), m1, m1.ld(), m2, m2.
ld(), beta, *
this);
33425 __inv<basic_srsmatrix>(*
this, m);
33461 __inv<basic_srsmatrix>(mRes, *
this);
33528 __exp_symm<basic_srsmatrix, TR>(*
this, m, tol);
33592 __exp_symm<basic_srsmatrix, TR>(msRes, *
this, tol);
33667 if (v.
incr() > 1) v1 << v;
33668 __polynom<TR,RVector>(this->
get(), this->
ld(), this->
msize(), m._pd(), m._ldm(), v.
incr() > 1 ? v1 : v);
33740 if (v.
incr() > 1) v1 << v;
33741 __polynom<TR,RVector>(msRes, msRes.
ld(), this->
msize(), this->
get(), this->
ld(), v.
incr() > 1 ? v1 : v);
33790 __eig<RVector, basic_srsmatrix, BaseSRMatrix>(vEig, *
this, &mEigVect,
true);
33827 __eig<RVector, basic_srsmatrix, BaseSRMatrix>(vEig, *
this,
nullptr,
true);
33873 tint nOutInfo = __cholesky<BaseSRMatrix>(mRes);
33876 mRes._clean_low_triangle();
33911 __bunch_kaufman<BaseSRMatrix>(mRes, nPivots);
33918 this->_plus_plus();
33930 this->_randomize(dFrom, dTo);
33936 return this->
norm1();
33944 bool is_positive_definite()
const
33946 static const TR zero = TR(0.);
33947 const TR* pd = this->_pv();
33948 const tint nSize = this->_size();
33949 const tint nNext = this->_msize() + 1;
33951 for (
tint i = 0; i < nSize; i += nNext) {
33952 if (pd[i] <= zero) {
33969 bool bRes = this->_equilibrate(vScalings);
33972 vB[i] *= vScalings[i];
33988 bool bRes = this->_equilibrate(vScalings);
33992 mB(i, j) *= vScalings[i];
34004 tint nOutInfo = __cholesky<BaseSRMatrix>(*this);
34006 if (nOutInfo > 0) {
34008 __bunch_kaufman<BaseSRMatrix>(*
this, nPivots);
34009 bPositiveDefinite =
false;
34011 bPositiveDefinite =
true;
34019 if (this->
msize() > 1) {
34020 const tint nM1 = this->
ld() + 1, nM1m = this->
ld() - 1, nM2m = this->
msize() - 1;
34021 tint i = 1, j = 1, m;
34023 m = this->
msize() - i;
34024 __copy<TR>(m, this->
get() + j + nM1m, this->
ld(), this->
get() + j, 1);
34037 return this->
get();
34041 const TR* _pd()
const override {
34042 return this->
get();
34045 void _check_submatrix()
const throw(cvmexception)
override {
34050 const TR* _pp(
const BaseMatrix& m)
const override {
34054 void _check_ger() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"ger",
"srsmatrix");}
34055 void _check_rank1update() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"rank1update",
"srsmatrix");}
34056 void _check_gemm() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"gemm",
"srsmatrix");}
34061 void _transp()
override {
34070 void _scalr(TR d)
override {
34071 __scal<TR,TR>(this->
get(), this->
size(), this->
incr(), d);
34074 void _multiply(RVector& vRes,
const RVector& v,
bool)
const override
34078 static const TR zero = TR(0.);
34079 static const TR one = TR(1.);
34080 const TR* pDm = this->
get();
34082 if (vRes.get() == pDv) vTmp << v;
34083 if (vRes.get() == pDm) mTmp << *
this;
34084 __symv<TR,basic_srsmatrix, RVector>(vRes.get() == pDm ? mTmp : *
this, one,
34085 vRes.
get() == pDv ? vTmp : v, zero, vRes);
34088 void _randomize(TR dFrom, TR dTo)
override
34090 this->BaseSRMatrix::_randomize(dFrom, dTo);
34094 void _solve(
const RVector& vB, RVector& vX,
34095 TR& dErr,
const TR* pLU,
const tint* pPivots,
int transp_mode)
const throw(cvmexception)
override
34098 if (!this->is_positive_definite()) {
34099 BaseSRMatrix::_solve(vB, vX, dErr, pLU, pPivots, transp_mode);
34103 RVector vScalings(this->
msize());
34105 const bool bEquilibrated = m.equilibrate(vScalings, vB1);
34107 __solve<TR,TR, basic_srsmatrix>(m, 1, vB1, vB1.size(), vX1, vX1.size(), dErr, pLU, pPivots, 0);
34108 if (bEquilibrated) {
34109 for (
tint i = 1; i <= this->
msize(); ++i) {
34110 vX[i] = vX1[i] * vScalings[i];
34117 void _solve(
const BaseRMatrix& mB, BaseRMatrix& mX,
34118 TR& dErr,
const TR* pLU,
const tint* pPivots,
int transp_mode)
const throw(cvmexception)
override
34121 if (!this->is_positive_definite()) {
34122 BaseSRMatrix::_solve(mB, mX, dErr, pLU, pPivots, transp_mode);
34125 BaseRMatrix mB1(mB);
34126 RVector vScalings(this->
msize());
34128 const bool bEquilibrated = m.equilibrate(vScalings, mB1);
34130 __solve<TR,TR, basic_srsmatrix>(m, mB.
nsize(), mB, mB.
ld(), mX, mX.ld(), dErr, pLU, pPivots, 0);
34131 if (bEquilibrated) {
34132 for (
tint j = 1; j <= mX.nsize(); ++j) {
34133 for (
tint i = 1; i <= this->
msize(); ++i) {
34134 mX(i, j) *= vScalings[i];
34141 TR _det()
const throw(cvmexception)
override
34144 switch (this->
msize()) {
34148 dDet = this->_ij_val(0, 0);
34151 dDet = this->_ij_val(0, 0) * this->_ij_val(1, 1) -
34152 this->_ij_val(1, 0) * this->_ij_val(0, 1);
34156 static const TR one(1.);
34157 bool bPositiveDefinite =
false;
34160 m._factorize(*
this, nPivots, bPositiveDefinite);
34164 if (bPositiveDefinite) {
34165 const RVector vUpDiag = m.diag(0);
34166 for (i = 1; i <= this->
msize(); ++i) {
34167 dDet *= vUpDiag[i] * vUpDiag[i];
34170 const RVector vEig = this->eig();
34171 for (i = 1; i <= this->
msize(); ++i) {
34176 catch (
const cvmexception& e) {
34185 bool _equilibrate(RVector& vScalings)
throw(cvmexception)
34191 static const TR sp = basic_cvmMachSp<TR>();
34192 static const TR sp_inv = TR(1.) / sp;
34194 __poequ<TR,basic_srsmatrix, RVector>(*
this, vScalings, dCond, dMax);
34196 if (dCond < TR(0.1) || _abs(dMax) <= sp || _abs(dMax) >= sp_inv) {
34198 for (
tint i = 0; i < this->
msize(); ++i) {
34199 for (
tint j = i; j < this->
msize(); ++j) {
34200 this->
get()[this->
ld() * j + i] *= vScalings[i +
CVM0] * vScalings[j +
CVM0];
34209 void _set_at(
tint nRow,
tint nCol, TR val)
throw(cvmexception)
34213 this->
get()[this->
ld() * nCol + nRow] = val;
34214 if (nRow != nCol) {
34215 this->
get()[this->
ld() * nRow + nCol] = val;
34219 void _check_symmetric(TR tol)
const throw(cvmexception)
34221 for (
tint j = 0; j < this->
nsize(); ++j) {
34222 for (
tint i = 0; i < this->
msize(); ++i) {
34223 if (i != j && _abs(this->
get()[this->
ld() * j + i] - this->
get()[this->
ld() * i + j]) > tol) {
34240 template<
typename TR,
typename TC>
34345 this->_check_hermitian(tol);
34386 this->_check_hermitian(tol);
34478 this->_check_hermitian(tol);
34510 __copy2<TR,TC>(this->
get(), v.
size(), this->
ld() + 1, v._pd(),
nullptr, v.
incr());
34537 __copy2<TR,TC>(this->
get(), this->
size(), this->
incr(), m._pd(),
nullptr);
34572 this->_check_hermitian(tol);
34607 this->_check_hermitian(tol);
34637 :
BaseSCMatrix(m._sub_pointer(nRowCol, nRowCol, nDim, nDim), nDim, m.
ld(), nDim * nDim)
34645 return this->_ij_val(nRow -
CVM0, nCol -
CVM0);
34652 return this->_col(nCol -
CVM0);
34659 return this->_row(nRow -
CVM0);
34665 return this->_diag(nDiag);
34672 __copy<TR>(this->
size(), __get_real_p<TR>(this->
get()), this->
incr() * 2, mRet, mRet.
incr());
34680 __copy<TR>(this->
size(), __get_imag_p<TR>(this->
get()), this->
incr() * 2, mRet, mRet.
incr());
34774 this->_assign(v, v.
incr());
34775 this->_check_hermitian(tol);
34811 this->_assign(pd, 1);
34812 this->_check_hermitian(tol);
34850 this->_assign_shifted(this->_sub_pointer_nocheck(nRowCol, nRowCol), m._pd(), m.
msize(), m.
nsize(), m.
ld());
34885 this->_set_at(nRow -
CVM0, nCol -
CVM0, c);
34925 const tint nD = _abs(nDiag);
34926 const tint nSize = vDiag.size();
34928 const tint nShift = nDiag * this->
ld();
34929 const tint nIncr = this->
ld() + 1;
34930 TC* pD1 = this->
get() + (nDiag > 0 ? nShift : nD);
34931 TC* pD2 = this->
get() + (nDiag > 0 ? nD : nShift);
34932 __copy<TC>(nSize, vDiag, vDiag.incr(), pD1, nIncr);
34933 __copy<TC>(nSize, vDiag, vDiag.incr(), pD2, nIncr);
34934 __conj<TC>(pD2, nSize, nIncr);
34968 __copy_real<TR,TC>(this->
get(), this->
msize(), this->
ld() + 1, vDiag, vDiag.incr());
35002 __copy_real<TR,TC>(this->
get(), this->
size(), this->
incr(), mRe._pd(), mRe.
incr());
35009 this->_set_real_number(d);
35015 this->_resize2(nNewDim, nNewDim);
35077 return !(this->operator == (m));
35227 this->_msum(m1, m2);
35264 this->_mdiff(m1, m2);
35354 static const TR mone(-1.);
35363 this->_plus_plus();
35370 this->_plus_plus();
35377 this->_minus_minus();
35384 this->_minus_minus();
35391 mRes._scalr(dMult);
35405 mRes._scalc(cMult);
35418 this->_scalr(dMult);
35431 this->_normalize();
35528 return this->transpose();
35596 this->_multiply(vRes, v,
false);
35604 static const TR zero = TR(0.);
35605 static const TR one = TR(1.);
35607 mRes._hemm(
true, *
this, m, one, zero);
35646 static const TR zero = TR(0.);
35647 static const TR one = TR(1.);
35649 mRes._hemm(
true, *
this, m, one, zero);
35656 return v % (*this);
35711 __herk<TR,TC, basic_schmatrix>(
false, alpha, 1, vc, vc.size(), beta, *
this);
35784 __herk<TR,TC, basic_schmatrix>(bTransp, alpha, bTransp ? m.
msize() : m.
nsize(), m, m.
ld(), beta, *
this);
35842 __her2k<TR,TC, basic_schmatrix>(
false, alpha, 1, v1c, v1c.size(), v2c, v2c.
size(), beta, *
this);
35936 __her2k<TR,TC, basic_schmatrix>(bTransp, alpha, bTransp ? m1.msize() : m1.nsize(), m1, m1.ld(), m2, m2.
ld(), beta, *
this);
35975 __inv<basic_schmatrix>(*
this, m);
36012 __inv<basic_schmatrix>(mRes, *
this);
36102 __exp_symm<basic_schmatrix, TR>(*
this, m, tol);
36190 __exp_symm<basic_schmatrix, TR>(msRes, *
this, tol);
36288 __polynom<TC, CVector>(this->
get(), this->
ld(), this->
msize(), m._pd(), m._ldm(), vc);
36384 __polynom<TC, CVector>(msRes, msRes.
ld(), this->
msize(), this->
get(), this->
ld(), vc);
36445 __eig<RVector, basic_schmatrix, BaseSCMatrix>(vEig, *
this, &mEigVect,
true);
36484 __eig<RVector, basic_schmatrix, BaseSCMatrix>(vEig, *
this,
nullptr,
true);
36531 tint nOutInfo = __cholesky<BaseSCMatrix>(mRes);
36534 mRes._clean_low_triangle();
36569 __bunch_kaufman<BaseSCMatrix>(mRes, nPivots);
36576 this->_plus_plus();
36588 this->_randomize_real(dFrom, dTo);
36594 this->_randomize_imag(dFrom, dTo);
36601 return this->
norm1();
36609 bool is_positive_definite()
const
36611 static const TR zero = TR(0.);
36612 const TC* pd = this->_pv();
36613 const tint nSize = this->_size();
36614 const tint nNext = this->_msize() + 1;
36616 for (
tint i = 0; i < nSize; i += nNext) {
36617 if (pd[i].real() <= zero) {
36634 bool bRes = this->_equilibrate(vScalings);
36637 vB[i] *= vScalings[i];
36653 bool bRes = this->_equilibrate(vScalings);
36657 mB(i, j) *= vScalings[i];
36669 tint nOutInfo = __cholesky<BaseSCMatrix>(*this);
36671 if (nOutInfo > 0) {
36673 __bunch_kaufman<BaseSCMatrix>(*
this, nPivots);
36674 bPositiveDefinite =
false;
36676 bPositiveDefinite =
true;
36684 if (this->
msize() > 1) {
36685 const tint nM1 = this->
ld() + 1, nM1m = this->
ld() - 1, nM2m = this->
msize() - 1;
36686 tint i = 1, j = 1, m;
36688 m = this->
msize() - i;
36689 __copy<TC>(m, this->
get() + j + nM1m, this->
ld(), this->
get() + j, 1);
36690 __conj<TC>(this->
get() + j, m, 1);
36691 if (i >= nM2m)
break;
36699 TC* _pd()
override {
36700 return this->
get();
36704 const TC* _pd()
const override {
36705 return this->
get();
36708 void _check_submatrix()
const throw(cvmexception)
override {
36713 const TC* _pp(
const BaseMatrix& m)
const override {
36724 void _set(TC) throw(cvmexception)
override {
36728 void _randomize_real(TR dFrom, TR dTo)
override
36730 this->BaseSCMatrix::_randomize_real(dFrom, dTo);
36732 this->_make_main_diag_real();
36735 void _randomize_imag(TR dFrom, TR dTo)
override
36737 this->BaseSCMatrix::_randomize_imag(dFrom, dTo);
36739 this->_make_main_diag_real();
36742 void _transp()
override {
36743 this->_sq_transp();
36747 void _conj()
override {
36757 void _scalr(TR d)
override {
36758 __scal<TR,TC>(this->
get(), this->
size(), this->
incr(), d);
36761 void _scalc(TC d)
override {
36762 __scal<TC, TC>(this->
get(), this->
size(), this->
incr(), d);
36765 void _multiply(CVector& vRes,
const CVector& v,
bool bLeft)
const override
36769 BaseCMatrix::_multiply(vRes, v, bLeft);
36773 static const TR zero = TR(0.);
36774 static const TR one = TR(1.);
36775 const TC* pDm = this->
get();
36777 if (vRes.get() == pDv) vTmp << v;
36778 if (vRes.get() == pDm) mTmp << *
this;
36779 __shmv<TC, basic_schmatrix, CVector>(vRes.get() == pDm ? mTmp : *
this, one,
36780 vRes.
get() == pDv ? vTmp : v, zero, vRes);
36786 void _set_at(
tint nRow,
tint nCol, TC val)
throw(cvmexception)
36790 if (nRow == nCol && _abs(val.imag()) > basic_cvmMachMin<TR>()) {
36793 this->
get()[this->
ld() * nCol + nRow] = val;
36794 if (nRow != nCol) {
36795 this->
get()[this->
ld() * nRow + nCol] = _conjugate(val);
36799 void _check_hermitian(TR tol)
const throw(cvmexception)
36803 for (
tint j = 0; j < this->
nsize(); ++j) {
36804 for (
tint i = 0; i < this->
msize(); ++i) {
36806 discr = _abs(this->
get()[this->
ld() * j + i].imag());
36812 c1 = this->
get()[this->
ld() * j + i];
36813 c2 = this->
get()[this->
ld() * i + j];
36814 if (!_conjugated<TR>(c1, c2, tol)) {
36815 throw cvmexception(
CVM_NOT_CONJUGATED, c1.real(), c1.imag(), c2.real(), c2.imag(), tol);
36821 void _solve(
const CVector& vB, CVector& vX, TR& dErr,
36822 const TC* pLU,
const tint* pPivots,
int transp_mode)
const throw(cvmexception)
override
36825 if (transp_mode == 1 || !this->is_positive_definite()) {
36826 BaseSCMatrix::_solve(vB, vX, dErr, pLU, pPivots, transp_mode);
36832 const bool bEquilibrated = m.equilibrate(vScalings, vB1);
36834 __solve<TR,TC, basic_schmatrix>(m, 1, vB1, vB1.size(), vX1, vX1.size(), dErr, pLU, pPivots, 0);
36836 if (bEquilibrated) {
36838 vX[i] = vX1[i] * vScalings[i];
36845 void _solve(
const BaseCMatrix& mB, BaseCMatrix& mX, TR& dErr,
36846 const TC* pLU,
const tint* pPivots,
int transp_mode)
const throw(cvmexception)
override
36849 if (transp_mode == 1 || !this->is_positive_definite()) {
36850 BaseSCMatrix::_solve(mB, mX, dErr, pLU, pPivots, transp_mode);
36853 BaseCMatrix mB1(mB);
36856 const bool bEquilibrated = m.equilibrate(vScalings, mB1);
36858 __solve<TR,TC, basic_schmatrix>(m, mB.
nsize(), mB, mB.
ld(), mX, mX.ld(), dErr, pLU, pPivots, 0);
36859 if (bEquilibrated) {
36862 mX(i, j) *= vScalings[i];
36869 TC _det()
const throw(cvmexception)
override
36872 switch (this->
msize()) {
36876 dDet = this->_ij_val(0, 0);
36879 dDet = this->_ij_val(0, 0) * this->_ij_val(1, 1) -
36880 this->_ij_val(1, 0) * this->_ij_val(0, 1);
36884 static const TC one(1., 0.);
36885 bool bPositiveDefinite =
false;
36888 m._factorize(*
this, nPivots, bPositiveDefinite);
36892 if (bPositiveDefinite) {
36893 const CVector vUpDiag = m.diag(0);
36895 dDet *= vUpDiag[i] * vUpDiag[i];
36904 catch (
const cvmexception& e) {
36919 static const TR sp = basic_cvmMachSp<TR>();
36920 static const TR sp_inv = TR(1.) / sp;
36922 __poequ<TR,basic_schmatrix, basic_rvector<TR> >(*
this, vScalings, dCond, dMax);
36924 if (dCond < TR(0.1) || _abs(dMax) <= sp || _abs(dMax) >= sp_inv) {
36926 for (
tint i = 0; i < this->
msize(); ++i) {
36927 for (
tint j = i; j < this->
msize(); ++j) {
36928 this->
get()[this->
ld() * j + i] *= vScalings[i +
CVM0] * vScalings[j +
CVM0];
36935 void _make_main_diag_real()
36937 static const TR zero = TR(0.);
36938 __scal<TR,TR>(__get_imag_p<TR>(this->
get()), this->
msize(), (this->
ld() + 1) * 2, zero);
36941 void _check_gerc() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"gerc",
"schmatrix");}
36942 void _check_geru() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"geru",
"schmatrix");}
36943 void _check_rank1update_c() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"rank1update_c",
"schmatrix");}
36944 void _check_rank1update_u() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"rank1update_u",
"schmatrix");}
36945 void _check_gemm() throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"gemm",
"schmatrix");}
36948 void _set_imag_number(TR) throw(cvmexception)
override {
throw cvmexception(
CVM_METHODNOTAVAILABLE,
"set_imag",
"schmatrix");}
36954 template<
typename TR>
36959 template<
typename TR>
36964 template<
typename TR>
36969 template<
typename TR>
36974 template<
typename TR>
36980 template<
typename TR,
typename TC>
36985 template<
typename TR,
typename TC>
36990 template<
typename TR,
typename TC>
36995 template<
typename TR,
typename TC>
37000 template<
typename TR,
typename TC>
37006 template<
typename TR,
typename TC>
37011 template<
typename TR,
typename TC>
37016 template<
typename TR,
typename TC>
37021 template<
typename TR,
typename TC>
37026 template<
typename TR,
typename TC>
37032 template<
typename TR>
37034 return v *
static_cast<TR
>(d);
37037 template<
typename TR>
37039 return m *
static_cast<TR
>(d);
37042 template<
typename TR>
37044 return m *
static_cast<TR
>(d);
37047 template<
typename TR>
37049 return m *
static_cast<TR
>(d);
37052 template<
typename TR>
37054 return m *
static_cast<TR
>(d);
37058 template<
typename TR,
typename TC>
37060 return v *
static_cast<TR
>(d);
37063 template<
typename TR,
typename TC>
37065 return m *
static_cast<TR
>(d);
37068 template<
typename TR,
typename TC>
37070 return m *
static_cast<TR
>(d);
37073 template<
typename TR,
typename TC>
37075 return m *
static_cast<TR
>(d);
37078 template<
typename TR,
typename TC>
37080 return m *
static_cast<TR
>(d);
37084 #if defined (CVM_FLOAT)
37085 typedef float treal;
37106 template<
typename TR>
37114 template<
typename TR,
typename TC>
37124 return basic_eye_real<treal>(nM);
37128 return basic_eye_complex<treal, tcomplex>(nM);
37151 return basic_cvmMachMin<treal>();
37176 return basic_cvmMachSp<treal>();
37181 #if !defined (CVM_STD_MUTEX)
37183 class CriticalSection
37185 #if defined (CVM_MT)
37189 #if defined (WIN32) || defined (_WIN32)
37190 #if defined (CVM_USE_CRITICAL_SECTION_NOT_MUTEX)
37191 ::CRITICAL_SECTION mCriticalSection;
37195 #else // POSIX Threads library assumed
37196 pthread_mutex_t mMutex;
37197 pthread_mutexattr_t mMutexAttr;
37203 ~CriticalSection();
37206 #if defined (CVM_MT)
37207 throw(cvmexception)
37211 #if defined (CVM_MT)
37212 throw(cvmexception)
37219 CriticalSection& mcs;
37221 Lock(CriticalSection& cs) : mcs(cs) {
37236 #if defined(CVM_USE_USER_LITERALS)
37238 constexpr std::complex<treal>
operator "" _i(
long double im)
37240 return { 0.,
static_cast<treal
>(im) };
37243 constexpr std::complex<treal>
operator "" _i(
unsigned long long int im)
37245 return { 0.,
static_cast<treal
>(im) };
37248 template<
typename T,
typename TR>
37249 inline std::complex<TR>
operator + (T re,
const std::complex<TR>& c) {
37250 return {
static_cast<TR
>(re) + c.real(), c.imag() };
37253 template<
typename T,
typename TR>
37254 inline std::complex<TR>
operator - (T re,
const std::complex<TR>& c) {
37255 return {
static_cast<TR
>(re) - c.real(), - c.imag() };
37258 template<
typename T,
typename TR>
37259 inline std::complex<TR>
operator + (
const std::complex<TR>& c,
const T& re) {
37260 return { c.real() +
static_cast<TR
>(re), c.imag() };
37263 template<
typename T,
typename TR>
37264 inline std::complex<TR>
operator - (
const std::complex<TR>& c,
const T& re) {
37265 return { c.real() -
static_cast<TR
>(re), c.imag() };
37275 #if !defined (_MSC_VER)
37276 # define XERBLA xerbla_
37281 void __stdcall
XERBLA(
const char* szSubName,
37282 #
if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
37285 const tint* pnParam)
throw(cvm::cvmexception);