30 mRowPivots(src.mRowPivots),
31 mPivotInverse(src.mPivotInverse),
32 mSwapVector(src.mSwapVector),
33 mIndependent(src.mIndependent)
52 C_INT LDA = std::max<C_INT>(1, NumCols);
68 for (i = 0; i < NumRows; i++)
161 JPVT.array(), TAU.
array(), WORK.
array(), &LWORK, &INFO);
165 LWORK = (
C_INT) WORK[0];
169 JPVT.array(), TAU.
array(), WORK.
array(), &LWORK, &INFO);
176 for (i = 0; i < NumRows; i++)
184 C_INT independent = 0;
186 if (fabs(M(0, 0)) != 0.0)
202 C_FLOAT64 RCOND =
std::max(NumRows, NumCols) * smax * std::numeric_limits< C_FLOAT64 >::epsilon();
204 C_FLOAT64 sminpr, s1, c1, smaxpr, s2, c2;
206 while (independent < mn && independent < maxRank)
208 dlaic1_(&imin, &independent, pIsmin, &smin, &M(independent, 0), &M(independent, independent), &sminpr, &s1, &c1);
209 dlaic1_(&imax, &independent, pIsmax, &smax, &M(independent, 0), &M(independent, independent), &smaxpr, &s2, &c2);
211 if (smaxpr * RCOND > sminpr)
216 for (
size_t i = 0; i < independent - 1; ++i)
222 *(pIsmin + independent) = c1;
223 *(pIsmax + independent) = c2;
235 resize(NumRows - independent, independent);
237 if (NumRows != independent && independent != 0)
301 dtrtri_(&cL, &cU, &independent, M.
array(), &LDA, &INFO);
313 for (j = 0; j < NumRows - independent; j++)
314 for (i = 0; i < independent; i++, pTmp1++)
316 pTmp2 = &M(j + independent, i);
322 for (k = i; k < independent; k++, pTmp2++, pTmp3 += NumCols)
327 *pTmp1 += *pTmp3 * *pTmp2;
330 if (fabs(*pTmp1) < 100.0 * std::numeric_limits< C_FLOAT64 >::epsilon()) *pTmp1 = 0.0;
348 for (
size_t i = 0; pCurrentIndex != pEnd; ++pCurrentIndex, ++i)
354 size_t * pCurrentColumn = CurrentColumn.
array();
361 for (; pCurrentIndex != pEnd; ++pPivot, ++pJPVT, ++pCurrentIndex, ++pCurrentColumn)
368 size_t * pToColumn = & CurrentColumn[*pToIndex];
369 size_t * pFromColumn = pCurrentColumn;
372 *pJPVT = *pToIndex + 1;
374 size_t tmp = *pFromIndex;
375 *pFromIndex = *pToIndex;
379 *pFromColumn = *pToColumn;
390 for (; pPivot != pPivotEnd; ++pPivot, ++Index)
428 for (; pRowP < pEnd; pRowP += LDC, pRowM += LDB)
430 memcpy(pRowP, pRowM,
sizeof(
C_FLOAT64) * LDC);
436 dgemm_(&T, &T, &M, &N, &K,
437 const_cast< C_FLOAT64 * >(&alpha),
438 const_cast< C_FLOAT64 * >(
array()), &LDA,
439 const_cast< C_FLOAT64 * >(m.
array()) + LDA, &LDB,
440 const_cast< C_FLOAT64 * >(&alpha), p.
array(), &LDC);
478 dgemm_(&T, &T, &M, &N, &K, &Alpha,
479 const_cast< C_FLOAT64 * >(m.
array()), &LDA,
480 const_cast< C_FLOAT64 * >(
array()), &LDB,
519 size_t i, imax = pivots.
size();
524 for (i = 0; i < imax; i++)
532 memcpy(Tmp.array(), matrix[to],
sizeof(
C_FLOAT64) * numCols);
536 memcpy(matrix[to], matrix[from],
sizeof(
C_FLOAT64) * numCols);
563 const C_INT & incr)
const
577 const_cast< C_INT * >(&incr));
591 mpNumIndependent(&A.getNumIndependent())
621 os <<
"Matrix(" << imax <<
"x" << jmax <<
")" << std::endl;
623 for (i = 0; i < imax; i++)
625 for (j = 0; j < jmax; j++)
626 os <<
"\t" << A(i, j);
bool build(const CMatrix< C_FLOAT64 > &matrix, size_t maxRank=C_INVALID_INDEX)
size_t getNumDependent() const
CVector< size_t > mRowPivots
virtual size_t numRows() const
bool applyColumnPivot(CMatrix< C_FLOAT64 > &matrix, const C_INT &incr) const
bool leftMultiply(const CMatrix< C_FLOAT64 > &M, CMatrix< C_FLOAT64 > &P) const
void resize(size_t size, const bool ©=false)
const size_t & getNumIndependent() const
bool rightMultiply(const C_FLOAT64 &alpha, const CMatrix< C_FLOAT64 > &M, CMatrix< C_FLOAT64 > &P) const
std::ostream & operator<<(std::ostream &os, const CLinkMatrixView &A)
CVector< C_INT > mSwapVector
int dtrtri_(char *uplo, char *diag, integer *n, doublereal *a, integer *lda, integer *info)
CLinkMatrixView(const CLinkMatrix &A)
static const C_FLOAT64 mUnit
const CVector< size_t > & getRowPivots() const
bool doRowPivot(CMatrix< C_FLOAT64 > &matrix) const
bool applyRowPivot(CVectorCore< CType > &vector) const
int dlaswp_(integer *n, doublereal *a, integer *lda, integer *k1, integer *k2, integer *ipiv, integer *incx)
int dgeqp3_(integer *m, integer *n, doublereal *a, integer *lda, integer *jpvt, doublereal *tau, doublereal *work, integer *lwork, integer *info)
int dlaic1_(integer *job, integer *j, doublereal *x, doublereal *sest, doublereal *w, doublereal *gamma, doublereal *sestpr, doublereal *s, doublereal *c__)
int dgemm_(char *transa, char *transb, integer *m, integer *n, integer *k, doublereal *alpha, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *beta, doublereal *c, integer *ldc)
virtual void resize(size_t rows, size_t cols, const bool ©=false)
bool doColumnPivot(CMatrix< C_FLOAT64 > &matrix) const
virtual size_t size() const
bool undoColumnPivot(CMatrix< C_FLOAT64 > &matrix) const
void completePivotInformation()
static const C_FLOAT64 mZero
const size_t * mpNumIndependent
virtual size_t numCols() const
CVector< size_t > mPivotInverse
CLinkMatrixView & operator=(const CLinkMatrixView &rhs)
bool undoRowPivot(CMatrix< C_FLOAT64 > &matrix) const
#define CONSTRUCTOR_TRACE