GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-14 04:35:57 +08:00
Code Issues Packages Projects Releases Wiki Activity

merge

Browse Source
This commit is contained in:
Gael Guennebaud 2015-08-09 21:24:20 +02:00
commit 23aab82c0c
6 changed files with 15 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Eigen/src/SparseLU/SparseLU.h
View File

@ -740,8 +740,8 @@ struct SparseLUMatrixUReturnType : internal::no_assignment_operator
} }
else else
{ {
Map<const Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) ); Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 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< Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) );
U = A.template triangularView<Upper>().solve(U); U = A.template triangularView<Upper>().solve(U);
} }

2
Eigen/src/SparseLU/SparseLUImpl.h
View File

@ -22,6 +22,8 @@ class SparseLUImpl
public: public:
typedef Matrix<Scalar,Dynamic,1> ScalarVector; typedef Matrix<Scalar,Dynamic,1> ScalarVector;
typedef Matrix<StorageIndex,Dynamic,1> IndexVector; 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 typename ScalarVector::RealScalar RealScalar;
typedef Ref<Matrix<Scalar,Dynamic,1> > BlockScalarVector; typedef Ref<Matrix<Scalar,Dynamic,1> > BlockScalarVector;
typedef Ref<Matrix<StorageIndex,Dynamic,1> > BlockIndexVector; typedef Ref<Matrix<StorageIndex,Dynamic,1> > BlockIndexVector;

6
Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h
View File

@ -239,7 +239,7 @@ void MappedSuperNodalMatrix<Scalar,Index_>::solveInPlace( MatrixBase<Dest>&X) co
Index n = int(X.rows()); Index n = int(X.rows());
Index nrhs = Index(X.cols()); Index nrhs = Index(X.cols());
const Scalar * Lval = valuePtr(); // Nonzero values 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(); work.setZero();
for (Index k = 0; k <= nsuper(); k ++) 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; Index lda = colIndexPtr()[fsupc+1] - luptr;
// Triangular solve // 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) ); Map< Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) );
U = A.template triangularView<UnitLower>().solve(U); U = A.template triangularView<UnitLower>().solve(U);
// Matrix-vector product // 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; work.block(0, 0, nrow, nrhs) = A * U;
//Begin Scatter //Begin Scatter

4
Eigen/src/SparseLU/SparseLU_column_bmod.h
View File

@ -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) // points to the beginning of jcol in snode L\U(jsupno)
ufirst = glu.xlusup(jcol) + d_fsupc; ufirst = glu.xlusup(jcol) + d_fsupc;
Index lda = glu.xlusup(jcol+1) - glu.xlusup(jcol); 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); VectorBlock<ScalarVector> u(glu.lusup, ufirst, nsupc);
u = A.template triangularView<UnitLower>().solve(u); 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); VectorBlock<ScalarVector> l(glu.lusup, ufirst+nsupc, nrow);
l.noalias() -= A * u; l.noalias() -= A * u;

4
Eigen/src/SparseLU/SparseLU_kernel_bmod.h
View File

@ -56,7 +56,7 @@ EIGEN_DONT_INLINE void LU_kernel_bmod<SegSizeAtCompileTime>::run(const Index seg
// Dense triangular solve -- start effective triangle // Dense triangular solve -- start effective triangle
luptr += lda * no_zeros + no_zeros; luptr += lda * no_zeros + no_zeros;
// Form Eigen matrix and vector // 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); Map<Matrix<Scalar,SegSizeAtCompileTime,1> > u(tempv.data(), segsize);
u = A.template triangularView<UnitLower>().solve(u); 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; luptr += segsize;
const Index PacketSize = internal::packet_traits<Scalar>::size; const Index PacketSize = internal::packet_traits<Scalar>::size;
Index ldl = internal::first_multiple(nrow, PacketSize); 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_default_aligned(tempv.data()+segsize, PacketSize); Index aligned_offset = internal::first_default_aligned(tempv.data()+segsize, PacketSize);
Index aligned_with_B_offset = (PacketSize-internal::first_default_aligned(B.data(), PacketSize))%PacketSize; Index aligned_with_B_offset = (PacketSize-internal::first_default_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) ); Map<Matrix<Scalar,Dynamic,1>, 0, OuterStride<> > l(tempv.data()+segsize+aligned_offset+aligned_with_B_offset, nrow, OuterStride<>(ldl) );

8
Eigen/src/SparseLU/SparseLU_panel_bmod.h
View File

@ -102,7 +102,7 @@ void SparseLUImpl<Scalar,StorageIndex>::panel_bmod(const Index m, const Index w,
if(nsupc >= 2) if(nsupc >= 2)
{ {
Index ldu = internal::first_multiple<Index>(u_rows, PacketSize); 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 // gather U
Index u_col = 0; 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); Index lda = glu.xlusup(fsupc+1) - glu.xlusup(fsupc);
no_zeros = (krep - u_rows + 1) - fsupc; no_zeros = (krep - u_rows + 1) - fsupc;
luptr += lda * no_zeros + no_zeros; 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); U = A.template triangularView<UnitLower>().solve(U);
// update // update
luptr += u_rows; 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); eigen_assert(tempv.size()>w*ldu + nrow*w + 1);
Index ldl = internal::first_multiple<Index>(nrow, PacketSize); Index ldl = internal::first_multiple<Index>(nrow, PacketSize);
Index offset = (PacketSize-internal::first_default_aligned(B.data(), PacketSize)) % PacketSize; Index offset = (PacketSize-internal::first_default_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(); L.setZero();
internal::sparselu_gemm<Scalar>(L.rows(), L.cols(), B.cols(), B.data(), B.outerStride(), U.data(), U.outerStride(), L.data(), L.outerStride()); internal::sparselu_gemm<Scalar>(L.rows(), L.cols(), B.cols(), B.data(), B.outerStride(), U.data(), U.outerStride(), L.data(), L.outerStride());