diff options
author | Christoph Hertzberg <chtz@informatik.uni-bremen.de> | 2015-08-07 23:10:56 +0200 |
---|---|---|
committer | Christoph Hertzberg <chtz@informatik.uni-bremen.de> | 2015-08-07 23:10:56 +0200 |
commit | cac6b23033807187382a5ace7d49d39a3436fe1b (patch) | |
tree | 8e00fae3f2f60998fd88a862cd20bc5eff3f4b14 /Eigen/src/SparseLU | |
parent | 7e0d7a76b868b702b20b896f69cea8c6a8f74730 (diff) |
bug #1053: SparseLU failed with EIGEN_DEFAULT_TO_ROW_MAJOR
Diffstat (limited to 'Eigen/src/SparseLU')
-rw-r--r-- | Eigen/src/SparseLU/SparseLU.h | 4 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLUImpl.h | 2 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h | 6 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_column_bmod.h | 4 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_kernel_bmod.h | 4 | ||||
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_panel_bmod.h | 8 |
6 files changed, 15 insertions, 13 deletions
diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index 71a699e51..8cdd29c7b 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -740,8 +740,8 @@ struct SparseLUMatrixUReturnType : internal::no_assignment_operator } else { - Map<const Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) ); - Map< Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); + Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) ); + Map< Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); U = A.template triangularView<Upper>().solve(U); } diff --git a/Eigen/src/SparseLU/SparseLUImpl.h b/Eigen/src/SparseLU/SparseLUImpl.h index 731d1652c..fc0cfc4de 100644 --- a/Eigen/src/SparseLU/SparseLUImpl.h +++ b/Eigen/src/SparseLU/SparseLUImpl.h @@ -22,6 +22,8 @@ class SparseLUImpl public: typedef Matrix<Scalar,Dynamic,1> ScalarVector; typedef Matrix<StorageIndex,Dynamic,1> IndexVector; + typedef Matrix<Scalar,Dynamic,Dynamic,ColMajor> ScalarMatrix; + typedef Map<ScalarMatrix, 0, OuterStride<> > MappedMatrixBlock; typedef typename ScalarVector::RealScalar RealScalar; typedef Ref<Matrix<Scalar,Dynamic,1> > BlockScalarVector; typedef Ref<Matrix<StorageIndex,Dynamic,1> > BlockIndexVector; diff --git a/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h b/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h index b37b93cf1..96e88fda1 100644 --- a/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h +++ b/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h @@ -239,7 +239,7 @@ void MappedSuperNodalMatrix<Scalar,Index_>::solveInPlace( MatrixBase<Dest>&X) co Index n = int(X.rows()); Index nrhs = Index(X.cols()); const Scalar * Lval = valuePtr(); // Nonzero values - Matrix<Scalar,Dynamic,Dynamic> work(n, nrhs); // working vector + Matrix<Scalar,Dynamic,Dynamic, ColMajor> work(n, nrhs); // working vector work.setZero(); for (Index k = 0; k <= nsuper(); k ++) { @@ -270,12 +270,12 @@ void MappedSuperNodalMatrix<Scalar,Index_>::solveInPlace( MatrixBase<Dest>&X) co Index lda = colIndexPtr()[fsupc+1] - luptr; // Triangular solve - Map<const Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) ); + Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) ); Map< Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); U = A.template triangularView<UnitLower>().solve(U); // Matrix-vector product - new (&A) Map<const Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > ( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) ); + new (&A) Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > ( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) ); work.block(0, 0, nrow, nrhs) = A * U; //Begin Scatter diff --git a/Eigen/src/SparseLU/SparseLU_column_bmod.h b/Eigen/src/SparseLU/SparseLU_column_bmod.h index be190997d..b57f06802 100644 --- a/Eigen/src/SparseLU/SparseLU_column_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_column_bmod.h @@ -163,11 +163,11 @@ Index SparseLUImpl<Scalar,StorageIndex>::column_bmod(const Index jcol, const Ind // points to the beginning of jcol in snode L\U(jsupno) ufirst = glu.xlusup(jcol) + d_fsupc; Index lda = glu.xlusup(jcol+1) - glu.xlusup(jcol); - Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(glu.lusup.data()[luptr]), nsupc, nsupc, OuterStride<>(lda) ); + MappedMatrixBlock A( &(glu.lusup.data()[luptr]), nsupc, nsupc, OuterStride<>(lda) ); VectorBlock<ScalarVector> u(glu.lusup, ufirst, nsupc); u = A.template triangularView<UnitLower>().solve(u); - new (&A) Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > ( &(glu.lusup.data()[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) ); + new (&A) MappedMatrixBlock ( &(glu.lusup.data()[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) ); VectorBlock<ScalarVector> l(glu.lusup, ufirst+nsupc, nrow); l.noalias() -= A * u; diff --git a/Eigen/src/SparseLU/SparseLU_kernel_bmod.h b/Eigen/src/SparseLU/SparseLU_kernel_bmod.h index 9513f8369..b2619c925 100644 --- a/Eigen/src/SparseLU/SparseLU_kernel_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_kernel_bmod.h @@ -56,7 +56,7 @@ EIGEN_DONT_INLINE void LU_kernel_bmod<SegSizeAtCompileTime>::run(const Index seg // Dense triangular solve -- start effective triangle luptr += lda * no_zeros + no_zeros; // Form Eigen matrix and vector - Map<Matrix<Scalar,SegSizeAtCompileTime,SegSizeAtCompileTime>, 0, OuterStride<> > A( &(lusup.data()[luptr]), segsize, segsize, OuterStride<>(lda) ); + Map<Matrix<Scalar,SegSizeAtCompileTime,SegSizeAtCompileTime, ColMajor>, 0, OuterStride<> > A( &(lusup.data()[luptr]), segsize, segsize, OuterStride<>(lda) ); Map<Matrix<Scalar,SegSizeAtCompileTime,1> > u(tempv.data(), segsize); u = A.template triangularView<UnitLower>().solve(u); @@ -65,7 +65,7 @@ EIGEN_DONT_INLINE void LU_kernel_bmod<SegSizeAtCompileTime>::run(const Index seg luptr += segsize; const Index PacketSize = internal::packet_traits<Scalar>::size; Index ldl = internal::first_multiple(nrow, PacketSize); - Map<Matrix<Scalar,Dynamic,SegSizeAtCompileTime>, 0, OuterStride<> > B( &(lusup.data()[luptr]), nrow, segsize, OuterStride<>(lda) ); + Map<Matrix<Scalar,Dynamic,SegSizeAtCompileTime, ColMajor>, 0, OuterStride<> > B( &(lusup.data()[luptr]), nrow, segsize, OuterStride<>(lda) ); Index aligned_offset = internal::first_aligned(tempv.data()+segsize, PacketSize); Index aligned_with_B_offset = (PacketSize-internal::first_aligned(B.data(), PacketSize))%PacketSize; Map<Matrix<Scalar,Dynamic,1>, 0, OuterStride<> > l(tempv.data()+segsize+aligned_offset+aligned_with_B_offset, nrow, OuterStride<>(ldl) ); diff --git a/Eigen/src/SparseLU/SparseLU_panel_bmod.h b/Eigen/src/SparseLU/SparseLU_panel_bmod.h index bd3cf87b9..2a5b15b09 100644 --- a/Eigen/src/SparseLU/SparseLU_panel_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_panel_bmod.h @@ -102,7 +102,7 @@ void SparseLUImpl<Scalar,StorageIndex>::panel_bmod(const Index m, const Index w, if(nsupc >= 2) { Index ldu = internal::first_multiple<Index>(u_rows, PacketSize); - Map<Matrix<Scalar,Dynamic,Dynamic>, Aligned, OuterStride<> > U(tempv.data(), u_rows, u_cols, OuterStride<>(ldu)); + Map<ScalarMatrix, Aligned, OuterStride<> > U(tempv.data(), u_rows, u_cols, OuterStride<>(ldu)); // gather U Index u_col = 0; @@ -136,17 +136,17 @@ void SparseLUImpl<Scalar,StorageIndex>::panel_bmod(const Index m, const Index w, Index lda = glu.xlusup(fsupc+1) - glu.xlusup(fsupc); no_zeros = (krep - u_rows + 1) - fsupc; luptr += lda * no_zeros + no_zeros; - Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A(glu.lusup.data()+luptr, u_rows, u_rows, OuterStride<>(lda) ); + MappedMatrixBlock A(glu.lusup.data()+luptr, u_rows, u_rows, OuterStride<>(lda) ); U = A.template triangularView<UnitLower>().solve(U); // update luptr += u_rows; - Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > B(glu.lusup.data()+luptr, nrow, u_rows, OuterStride<>(lda) ); + MappedMatrixBlock B(glu.lusup.data()+luptr, nrow, u_rows, OuterStride<>(lda) ); eigen_assert(tempv.size()>w*ldu + nrow*w + 1); Index ldl = internal::first_multiple<Index>(nrow, PacketSize); Index offset = (PacketSize-internal::first_aligned(B.data(), PacketSize)) % PacketSize; - Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > L(tempv.data()+w*ldu+offset, nrow, u_cols, OuterStride<>(ldl)); + MappedMatrixBlock L(tempv.data()+w*ldu+offset, nrow, u_cols, OuterStride<>(ldl)); L.setZero(); internal::sparselu_gemm<Scalar>(L.rows(), L.cols(), B.cols(), B.data(), B.outerStride(), U.data(), U.outerStride(), L.data(), L.outerStride()); |