GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-12 11:49:02 +08:00
Code Issues Packages Projects Releases Wiki Activity

BDCSVD: Streamline compute() and copyUV()

Browse Source
This commit is contained in:
Jitse Niesen 2013-08-07 16:34:34 +01:00
parent 616f9cc593
commit 306ce33e1c
1 changed files with 20 additions and 41 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

57
unsupported/Eigen/src/SVD/BDCSVD.h
View File

@ -139,7 +139,7 @@ public:
void setSwitchSize(int s) void setSwitchSize(int s)
{ {
eigen_assert(s>3 && "BDCSVD the size of the algo switch has to be greater than 4"); eigen_assert(s>3 && "BDCSVD the size of the algo switch has to be greater than 3");
algoswap = s; algoswap = s;
} }
@ -201,7 +201,7 @@ private:
void deflation43(Index firstCol, Index shift, Index i, Index size); void deflation43(Index firstCol, Index shift, Index i, Index size);
void deflation44(Index firstColu , Index firstColm, Index firstRowW, Index firstColW, Index i, Index j, Index size); void deflation44(Index firstColu , Index firstColm, Index firstRowW, Index firstColW, Index i, Index j, Index size);
void deflation(Index firstCol, Index lastCol, Index k, Index firstRowW, Index firstColW, Index shift); void deflation(Index firstCol, Index lastCol, Index k, Index firstRowW, Index firstColW, Index shift);
void copyUV(MatrixXr naiveU, MatrixXr naiveV, MatrixX householderU, MatrixX houseHolderV); void copyUV(MatrixX householderU, MatrixX houseHolderV);
protected: protected:
MatrixXr m_naiveU, m_naiveV; MatrixXr m_naiveU, m_naiveV;
@ -282,15 +282,8 @@ BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsigned int computationOp
internal::UpperBidiagonalization<MatrixX > bid(copy); internal::UpperBidiagonalization<MatrixX > bid(copy);
//**** step 2 Divide //**** step 2 Divide
// this is ugly and has to be redone (care of complex cast) m_computed.topRows(this->m_diagSize) = bid.bidiagonal().toDenseMatrix().transpose();
MatrixXr temp; m_computed.template bottomRows<1>().setZero();
temp = bid.bidiagonal().toDenseMatrix().transpose();
m_computed.setZero();
for (int i=0; i<this->m_diagSize - 1; i++) {
m_computed(i, i) = temp(i, i);
m_computed(i + 1, i) = temp(i + 1, i);
}
m_computed(this->m_diagSize - 1, this->m_diagSize - 1) = temp(this->m_diagSize - 1, this->m_diagSize - 1);
divide(0, this->m_diagSize - 1, 0, 0, 0); divide(0, this->m_diagSize - 1, 0, 0, 0);
//**** step 3 copy //**** step 3 copy
@ -307,46 +300,32 @@ BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsigned int computationOp
break; break;
} }
} }
copyUV(m_naiveV, m_naiveU, bid.householderU(), bid.householderV()); copyUV(bid.householderU(), bid.householderV());
this->m_isInitialized = true; this->m_isInitialized = true;
return *this; return *this;
}// end compute }// end compute
// TODO: this function should accept householder sequences to save converting them to matrix
template<typename MatrixType> template<typename MatrixType>
void BDCSVD<MatrixType>::copyUV(MatrixXr naiveU, MatrixXr naiveV, MatrixX householderU, MatrixX householderV){ void BDCSVD<MatrixType>::copyUV(MatrixX householderU, MatrixX householderV){
// Note exchange of U and V: m_matrixU is set from m_naiveV and vice versa
if (this->computeU()){ if (this->computeU()){
MatrixX temp = MatrixX::Zero(naiveU.rows(), naiveU.cols()); Index Ucols = this->m_computeThinU ? this->m_nonzeroSingularValues : householderU.cols();
temp.real() = naiveU; this->m_matrixU = MatrixX::Identity(householderU.cols(), Ucols);
if (this->m_computeThinU){ Index blockCols = this->m_computeThinU ? this->m_nonzeroSingularValues : this->m_diagSize;
this->m_matrixU = MatrixX::Identity(householderU.cols(), this->m_nonzeroSingularValues ); this->m_matrixU.block(0, 0, this->m_diagSize, blockCols) =
this->m_matrixU.block(0, 0, this->m_diagSize, this->m_nonzeroSingularValues) = m_naiveV.template cast<Scalar>().block(0, 0, this->m_diagSize, blockCols);
temp.block(0, 0, this->m_diagSize, this->m_nonzeroSingularValues);
this->m_matrixU = householderU * this->m_matrixU; this->m_matrixU = householderU * this->m_matrixU;
} }
else
{
this->m_matrixU = MatrixX::Identity(householderU.cols(), householderU.cols());
this->m_matrixU.block(0, 0, this->m_diagSize, this->m_diagSize) = temp.block(0, 0, this->m_diagSize, this->m_diagSize);
this->m_matrixU = householderU * this->m_matrixU ;
}
}
if (this->computeV()){ if (this->computeV()){
MatrixX temp = MatrixX::Zero(naiveV.rows(), naiveV.cols()); Index Vcols = this->m_computeThinV ? this->m_nonzeroSingularValues : householderV.cols();
temp.real() = naiveV; this->m_matrixV = MatrixX::Identity(householderV.cols(), Vcols);
if (this->m_computeThinV){ Index blockCols = this->m_computeThinV ? this->m_nonzeroSingularValues : this->m_diagSize;
this->m_matrixV = MatrixX::Identity(householderV.cols(),this->m_nonzeroSingularValues ); this->m_matrixV.block(0, 0, this->m_diagSize, blockCols) =
this->m_matrixV.block(0, 0, this->m_nonzeroSingularValues, this->m_nonzeroSingularValues) = m_naiveU.template cast<Scalar>().block(0, 0, this->m_diagSize, blockCols);
temp.block(0, 0, this->m_nonzeroSingularValues, this->m_nonzeroSingularValues);
this->m_matrixV = householderV * this->m_matrixV; this->m_matrixV = householderV * this->m_matrixV;
} }
else
{
this->m_matrixV = MatrixX::Identity(householderV.cols(), householderV.cols());
this->m_matrixV.block(0, 0, this->m_diagSize, this->m_diagSize) = temp.block(0, 0, this->m_diagSize, this->m_diagSize);
this->m_matrixV = householderV * this->m_matrixV;
}
}
} }
// The divide algorithm is done "in place", we are always working on subsets of the same matrix. The divide methods takes as argument the // The divide algorithm is done "in place", we are always working on subsets of the same matrix. The divide methods takes as argument the