aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/Householder1
-rw-r--r--Eigen/src/Householder/BlockHouseholder.h76
-rw-r--r--Eigen/src/QR/HouseholderQR.h58
-rw-r--r--test/qr.cpp8
4 files changed, 137 insertions, 6 deletions
diff --git a/Eigen/Householder b/Eigen/Householder
index dcaeeeae8..11a6eefd1 100644
--- a/Eigen/Householder
+++ b/Eigen/Householder
@@ -17,6 +17,7 @@ namespace Eigen {
#include "src/Householder/Householder.h"
#include "src/Householder/HouseholderSequence.h"
+#include "src/Householder/BlockHouseholder.h"
} // namespace Eigen
diff --git a/Eigen/src/Householder/BlockHouseholder.h b/Eigen/src/Householder/BlockHouseholder.h
new file mode 100644
index 000000000..c17ca10e0
--- /dev/null
+++ b/Eigen/src/Householder/BlockHouseholder.h
@@ -0,0 +1,76 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Vincent Lejeune
+// Copyright (C) 2010 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_BLOCK_HOUSEHOLDER_H
+#define EIGEN_BLOCK_HOUSEHOLDER_H
+
+// This file contains some helper function to deal with block householder reflectors
+
+/** \internal */
+template<typename TriangularFactorType,typename VectorsType,typename CoeffsType>
+void ei_make_block_householder_triangular_factor(TriangularFactorType& triFactor, const VectorsType& vectors, const CoeffsType& hCoeffs)
+{
+ typedef typename TriangularFactorType::Index Index;
+ typedef typename VectorsType::Scalar Scalar;
+ const Index nbVecs = vectors.cols();
+ ei_assert(triFactor.rows() == nbVecs && triFactor.cols() == nbVecs && vectors.rows()>=nbVecs);
+
+ for(Index i = 0; i < nbVecs; i++)
+ {
+ Index rs = vectors.rows() - i;
+ Scalar Vii = vectors(i,i);
+ vectors.const_cast_derived().coeffRef(i,i) = Scalar(1);
+ triFactor.col(i).head(i).noalias() = vectors.block(i, 0, rs, i).adjoint()
+ * vectors.col(i).tail(rs);
+ triFactor.col(i).head(i) *= -hCoeffs(i);
+ vectors.const_cast_derived().coeffRef(i, i) = Vii;
+ // FIXME add .noalias() once the triangular product can work inplace
+ triFactor.col(i).head(i) = triFactor.block(0,0,i,i).template triangularView<Upper>()
+ * triFactor.col(i).head(i);
+ triFactor(i,i) = hCoeffs(i);
+ }
+}
+
+/** \internal */
+template<typename MatrixType,typename VectorsType,typename CoeffsType>
+void ei_apply_block_householder_on_the_left(MatrixType& mat, const VectorsType& vectors, const CoeffsType& hCoeffs)
+{
+ typedef typename MatrixType::Index Index;
+ enum { TFactorSize = MatrixType::ColsAtCompileTime };
+ Index nbVecs = vectors.cols();
+ Matrix<typename MatrixType::Scalar, TFactorSize, TFactorSize> T(nbVecs,nbVecs);
+ ei_make_block_householder_triangular_factor(T, vectors, hCoeffs);
+
+ const TriangularView<VectorsType, UnitLower>& V(vectors);
+
+ // A -= V T V^* A
+ Matrix<typename MatrixType::Scalar,Dynamic,Dynamic> tmp = V.adjoint() * mat;
+ // FIXME add .noalias() once the triangular product can work inplace
+ tmp = T.template triangularView<Upper>().adjoint() * tmp;
+ mat.noalias() -= V * tmp;
+}
+
+
+#endif // EIGEN_BLOCK_HOUSEHOLDER_H
diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h
index c977c2e33..29044e905 100644
--- a/Eigen/src/QR/HouseholderQR.h
+++ b/Eigen/src/QR/HouseholderQR.h
@@ -1,8 +1,9 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2008-2010 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2010 Vincent Lejeune
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@@ -229,6 +230,59 @@ void ei_householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typena
}
}
+/** \internal */
+template<typename MatrixQR, typename HCoeffs>
+void ei_householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs,
+ typename MatrixQR::Index maxBlockSize=32,
+ typename MatrixQR::Scalar* tempData = 0)
+{
+ typedef typename MatrixQR::Index Index;
+ typedef typename MatrixQR::Scalar Scalar;
+ typedef typename MatrixQR::RealScalar RealScalar;
+ typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
+
+ Index rows = mat.rows();
+ Index cols = mat.cols();
+ Index size = std::min(rows, cols);
+
+ typedef Matrix<Scalar,Dynamic,1,0,MatrixQR::MaxColsAtCompileTime,1> TempType;
+ TempType tempVector;
+ if(tempData==0)
+ {
+ tempVector.resize(cols);
+ tempData = tempVector.data();
+ }
+
+ Index blockSize = std::min(maxBlockSize,size);
+
+ int k = 0;
+ for (k = 0; k < size; k += blockSize)
+ {
+ Index bs = std::min(size-k,blockSize); // actual size of the block
+ Index tcols = cols - k - bs; // trailing columns
+ Index brows = rows-k; // rows of the block
+
+ // partition the matrix:
+ // A00 | A01 | A02
+ // mat = A10 | A11 | A12
+ // A20 | A21 | A22
+ // and performs the qr dec of [A11^T A12^T]^T
+ // and update [A21^T A22^T]^T using level 3 operations.
+ // Finally, the algorithm continue on A22
+
+ BlockType A11_21 = mat.block(k,k,brows,bs);
+ Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
+
+ ei_householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
+
+ if(tcols)
+ {
+ BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
+ ei_apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint());
+ }
+ }
+}
+
template<typename MatrixType>
HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
{
@@ -241,7 +295,7 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
m_temp.resize(cols);
- ei_householder_qr_inplace_unblocked(m_qr, m_hCoeffs, m_temp.data());
+ ei_householder_qr_inplace_blocked(m_qr, m_hCoeffs, 48, m_temp.data());
m_isInitialized = true;
return *this;
diff --git a/test/qr.cpp b/test/qr.cpp
index 90209639f..eae336e49 100644
--- a/test/qr.cpp
+++ b/test/qr.cpp
@@ -36,10 +36,10 @@ template<typename MatrixType> void qr(const MatrixType& m)
MatrixType a = MatrixType::Random(rows,cols);
HouseholderQR<MatrixType> qrOfA(a);
-
+
MatrixQType q = qrOfA.householderQ();
VERIFY_IS_UNITARY(q);
-
+
MatrixType r = qrOfA.matrixQR().template triangularView<Upper>();
VERIFY_IS_APPROX(a, qrOfA.householderQ() * r);
}
@@ -112,8 +112,8 @@ template<typename MatrixType> void qr_verify_assert()
void test_qr()
{
for(int i = 0; i < g_repeat; i++) {
- CALL_SUBTEST_1( qr(MatrixXf(47,40)) );
- CALL_SUBTEST_2( qr(MatrixXcd(17,7)) );
+ CALL_SUBTEST_1( qr(MatrixXf(ei_random<int>(1,200),ei_random<int>(1,200))) );
+ CALL_SUBTEST_2( qr(MatrixXcd(ei_random<int>(1,200),ei_random<int>(1,200))) );
CALL_SUBTEST_3(( qr_fixedsize<Matrix<float,3,4>, 2 >() ));
CALL_SUBTEST_4(( qr_fixedsize<Matrix<double,6,2>, 4 >() ));
CALL_SUBTEST_5(( qr_fixedsize<Matrix<double,2,5>, 7 >() ));
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-06 09:28:53 +0000