aboutsummaryrefslogtreecommitdiffhomepage
path: root/unsupported
diff options
context:
space:
mode:
authorGravatar Jitse Niesen <jitse@maths.leeds.ac.uk>2013-08-07 16:34:34 +0100
committerGravatar Jitse Niesen <jitse@maths.leeds.ac.uk>2013-08-07 16:34:34 +0100
commit306ce33e1c08716005554228dd17138f11a778c0 (patch)
tree42e2e10e7ae8b2c6486736b76491d3e60fa48581 /unsupported
parent616f9cc593fb47e39e0a47b6a749f15d66a5c734 (diff)
BDCSVD: Streamline compute() and copyUV()
Diffstat (limited to 'unsupported')
-rw-r--r--unsupported/Eigen/src/SVD/BDCSVD.h61
1 files changed, 20 insertions, 41 deletions
diff --git a/unsupported/Eigen/src/SVD/BDCSVD.h b/unsupported/Eigen/src/SVD/BDCSVD.h
index 11d4882e4..dab419a47 100644
--- a/unsupported/Eigen/src/SVD/BDCSVD.h
+++ b/unsupported/Eigen/src/SVD/BDCSVD.h
@@ -139,7 +139,7 @@ public:
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;
}
@@ -201,7 +201,7 @@ private:
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 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:
MatrixXr m_naiveU, m_naiveV;
@@ -282,15 +282,8 @@ BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsigned int computationOp
internal::UpperBidiagonalization<MatrixX > bid(copy);
//**** step 2 Divide
- // this is ugly and has to be redone (care of complex cast)
- MatrixXr temp;
- 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);
+ m_computed.topRows(this->m_diagSize) = bid.bidiagonal().toDenseMatrix().transpose();
+ m_computed.template bottomRows<1>().setZero();
divide(0, this->m_diagSize - 1, 0, 0, 0);
//**** step 3 copy
@@ -307,45 +300,31 @@ BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsigned int computationOp
break;
}
}
- copyUV(m_naiveV, m_naiveU, bid.householderU(), bid.householderV());
+ copyUV(bid.householderU(), bid.householderV());
this->m_isInitialized = true;
return *this;
}// end compute
+// TODO: this function should accept householder sequences to save converting them to matrix
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()){
- MatrixX temp = MatrixX::Zero(naiveU.rows(), naiveU.cols());
- temp.real() = naiveU;
- if (this->m_computeThinU){
- this->m_matrixU = MatrixX::Identity(householderU.cols(), this->m_nonzeroSingularValues );
- this->m_matrixU.block(0, 0, this->m_diagSize, this->m_nonzeroSingularValues) =
- temp.block(0, 0, this->m_diagSize, this->m_nonzeroSingularValues);
- 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 ;
- }
+ Index Ucols = this->m_computeThinU ? this->m_nonzeroSingularValues : householderU.cols();
+ this->m_matrixU = MatrixX::Identity(householderU.cols(), Ucols);
+ Index blockCols = this->m_computeThinU ? this->m_nonzeroSingularValues : this->m_diagSize;
+ this->m_matrixU.block(0, 0, this->m_diagSize, blockCols) =
+ m_naiveV.template cast<Scalar>().block(0, 0, this->m_diagSize, blockCols);
+ this->m_matrixU = householderU * this->m_matrixU;
}
if (this->computeV()){
- MatrixX temp = MatrixX::Zero(naiveV.rows(), naiveV.cols());
- temp.real() = naiveV;
- if (this->m_computeThinV){
- this->m_matrixV = MatrixX::Identity(householderV.cols(),this->m_nonzeroSingularValues );
- this->m_matrixV.block(0, 0, this->m_nonzeroSingularValues, this->m_nonzeroSingularValues) =
- temp.block(0, 0, this->m_nonzeroSingularValues, this->m_nonzeroSingularValues);
- 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;
- }
+ Index Vcols = this->m_computeThinV ? this->m_nonzeroSingularValues : householderV.cols();
+ this->m_matrixV = MatrixX::Identity(householderV.cols(), Vcols);
+ Index blockCols = this->m_computeThinV ? this->m_nonzeroSingularValues : this->m_diagSize;
+ this->m_matrixV.block(0, 0, this->m_diagSize, blockCols) =
+ m_naiveU.template cast<Scalar>().block(0, 0, this->m_diagSize, blockCols);
+ this->m_matrixV = householderV * this->m_matrixV;
}
}