CFD/eigen
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Make UpperBidiagonalization accept row-major matrices (bug #769)

Browse Source 
* Give temporary workspace the same storage order as original matrix
* Take storage order into account when determining inner stride
  of rows and columns
* Change one test to use a row-major matrix.
This commit is contained in:
Jitse Niesen 2014-09-12 14:52:35 +01:00
parent e6d55c081b
commit 9452eb38f8
2 changed files with 23 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
Eigen/src/SVD/UpperBidiagonalization.h
View File

@ -154,14 +154,19 @@ void upperbidiagonalization_blocked_helper(MatrixType& A,
typename MatrixType::RealScalar *diagonal,
typename MatrixType::RealScalar *upper_diagonal,
typename MatrixType::Index bs,
Ref<Matrix<typename MatrixType::Scalar, Dynamic, Dynamic> > X,
Ref<Matrix<typename MatrixType::Scalar, Dynamic, Dynamic> > Y)
Ref<Matrix<typename MatrixType::Scalar, Dynamic, Dynamic,
traits<MatrixType>::Flags & RowMajorBit> > X,
Ref<Matrix<typename MatrixType::Scalar, Dynamic, Dynamic,
traits<MatrixType>::Flags & RowMajorBit> > Y)
{
typedef typename MatrixType::Index Index;
typedef typename MatrixType::Scalar Scalar;
typedef Ref<Matrix<Scalar, Dynamic, 1> > SubColumnType;
typedef Ref<Matrix<Scalar, 1, Dynamic>, 0, InnerStride<> > SubRowType;
typedef Ref<Matrix<Scalar, Dynamic, Dynamic> > SubMatType;
enum { StorageOrder = traits<MatrixType>::Flags & RowMajorBit };
typedef InnerStride<int(StorageOrder) == int(ColMajor) ? 1 : Dynamic> ColInnerStride;
typedef InnerStride<int(StorageOrder) == int(ColMajor) ? Dynamic : 1> RowInnerStride;
typedef Ref<Matrix<Scalar, Dynamic, 1>, 0, ColInnerStride> SubColumnType;
typedef Ref<Matrix<Scalar, 1, Dynamic>, 0, RowInnerStride> SubRowType;
typedef Ref<Matrix<Scalar, Dynamic, Dynamic, StorageOrder > > SubMatType;
Index brows = A.rows();
Index bcols = A.cols();
@ -288,8 +293,18 @@ void upperbidiagonalization_inplace_blocked(MatrixType& A, BidiagType& bidiagona
Index cols = A.cols();
Index size = (std::min)(rows, cols);
Matrix<Scalar,MatrixType::RowsAtCompileTime,Dynamic,ColMajor,MatrixType::MaxRowsAtCompileTime> X(rows,maxBlockSize);
Matrix<Scalar,MatrixType::ColsAtCompileTime,Dynamic,ColMajor,MatrixType::MaxColsAtCompileTime> Y(cols,maxBlockSize);
// X and Y are work space
enum { StorageOrder = traits<MatrixType>::Flags & RowMajorBit };
Matrix<Scalar,
MatrixType::RowsAtCompileTime,
Dynamic,
StorageOrder,
MatrixType::MaxRowsAtCompileTime> X(rows,maxBlockSize);
Matrix<Scalar,
MatrixType::ColsAtCompileTime,
Dynamic,
StorageOrder,
MatrixType::MaxColsAtCompileTime> Y(cols,maxBlockSize);
Index blockSize = (std::min)(maxBlockSize,size);
Index k = 0;

2
test/upperbidiagonalization.cpp
View File

@ -35,7 +35,7 @@ void test_upperbidiagonalization()
CALL_SUBTEST_1( upperbidiag(MatrixXf(3,3)) );
CALL_SUBTEST_2( upperbidiag(MatrixXd(17,12)) );
CALL_SUBTEST_3( upperbidiag(MatrixXcf(20,20)) );
CALL_SUBTEST_4( upperbidiag(MatrixXcd(16,15)) );
CALL_SUBTEST_4( upperbidiag(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(16,15)) );
CALL_SUBTEST_5( upperbidiag(Matrix<float,6,4>()) );
CALL_SUBTEST_6( upperbidiag(Matrix<float,5,5>()) );
CALL_SUBTEST_7( upperbidiag(Matrix<double,4,3>()) );