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

add support for strictly triangular matrix in trmm though it is not really useful

Browse Source
This commit is contained in:
Gael Guennebaud 2010-07-15 20:39:20 +02:00
parent cf9edd9958
commit 84fdbded4d
2 changed files with 23 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
Eigen/src/Core/products/TriangularMatrixMatrix.h
View File

@ -75,7 +75,7 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,LhsIsTriangular,
Scalar alpha) Scalar alpha)
{ {
ei_product_triangular_matrix_matrix<Scalar, Index, ei_product_triangular_matrix_matrix<Scalar, Index,
(Mode&UnitDiag) | ((Mode&Upper) ? Lower : Upper), (Mode&(UnitDiag|ZeroDiag)) | ((Mode&Upper) ? Lower : Upper),
(!LhsIsTriangular), (!LhsIsTriangular),
RhsStorageOrder==RowMajor ? ColMajor : RowMajor, RhsStorageOrder==RowMajor ? ColMajor : RowMajor,
ConjugateRhs, ConjugateRhs,
@ -111,7 +111,8 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,true,
typedef ei_product_blocking_traits<Scalar> Blocking; typedef ei_product_blocking_traits<Scalar> Blocking;
enum { enum {
SmallPanelWidth = EIGEN_PLAIN_ENUM_MAX(Blocking::mr,Blocking::nr), SmallPanelWidth = EIGEN_PLAIN_ENUM_MAX(Blocking::mr,Blocking::nr),
IsLower = (Mode&Lower) == Lower IsLower = (Mode&Lower) == Lower,
SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1
}; };
Index kc = depth; // cache block size along the K direction Index kc = depth; // cache block size along the K direction
@ -127,7 +128,10 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,true,
Matrix<Scalar,SmallPanelWidth,SmallPanelWidth,LhsStorageOrder> triangularBuffer; Matrix<Scalar,SmallPanelWidth,SmallPanelWidth,LhsStorageOrder> triangularBuffer;
triangularBuffer.setZero(); triangularBuffer.setZero();
triangularBuffer.diagonal().setOnes(); if((Mode&ZeroDiag)==ZeroDiag)
triangularBuffer.diagonal().setZero();
else
triangularBuffer.diagonal().setOnes();
ei_gebp_kernel<Scalar, Index, Blocking::mr, Blocking::nr, ConjugateLhs, ConjugateRhs> gebp_kernel; ei_gebp_kernel<Scalar, Index, Blocking::mr, Blocking::nr, ConjugateLhs, ConjugateRhs> gebp_kernel;
ei_gemm_pack_lhs<Scalar, Index, Blocking::mr,LhsStorageOrder> pack_lhs; ei_gemm_pack_lhs<Scalar, Index, Blocking::mr,LhsStorageOrder> pack_lhs;
@ -169,7 +173,7 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,true,
// To this end we do an extra triangular copy to a small temporary buffer // To this end we do an extra triangular copy to a small temporary buffer
for (Index k=0;k<actualPanelWidth;++k) for (Index k=0;k<actualPanelWidth;++k)
{ {
if (!(Mode&UnitDiag)) if (SetDiag)
triangularBuffer.coeffRef(k,k) = lhs(startBlock+k,startBlock+k); triangularBuffer.coeffRef(k,k) = lhs(startBlock+k,startBlock+k);
for (Index i=IsLower ? k+1 : 0; IsLower ? i<actualPanelWidth : i<k; ++i) for (Index i=IsLower ? k+1 : 0; IsLower ? i<actualPanelWidth : i<k; ++i)
triangularBuffer.coeffRef(i,k) = lhs(startBlock+i,startBlock+k); triangularBuffer.coeffRef(i,k) = lhs(startBlock+i,startBlock+k);
@ -237,7 +241,8 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,false,
typedef ei_product_blocking_traits<Scalar> Blocking; typedef ei_product_blocking_traits<Scalar> Blocking;
enum { enum {
SmallPanelWidth = EIGEN_PLAIN_ENUM_MAX(Blocking::mr,Blocking::nr), SmallPanelWidth = EIGEN_PLAIN_ENUM_MAX(Blocking::mr,Blocking::nr),
IsLower = (Mode&Lower) == Lower IsLower = (Mode&Lower) == Lower,
SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1
}; };
Index kc = depth; // cache block size along the K direction Index kc = depth; // cache block size along the K direction
@ -252,7 +257,10 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,false,
Matrix<Scalar,SmallPanelWidth,SmallPanelWidth,RhsStorageOrder> triangularBuffer; Matrix<Scalar,SmallPanelWidth,SmallPanelWidth,RhsStorageOrder> triangularBuffer;
triangularBuffer.setZero(); triangularBuffer.setZero();
triangularBuffer.diagonal().setOnes(); if((Mode&ZeroDiag)==ZeroDiag)
triangularBuffer.diagonal().setZero();
else
triangularBuffer.diagonal().setOnes();
ei_gebp_kernel<Scalar, Index, Blocking::mr, Blocking::nr, ConjugateLhs, ConjugateRhs> gebp_kernel; ei_gebp_kernel<Scalar, Index, Blocking::mr, Blocking::nr, ConjugateLhs, ConjugateRhs> gebp_kernel;
ei_gemm_pack_lhs<Scalar, Index, Blocking::mr,LhsStorageOrder> pack_lhs; ei_gemm_pack_lhs<Scalar, Index, Blocking::mr,LhsStorageOrder> pack_lhs;
@ -300,7 +308,7 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,false,
// append the triangular part via a temporary buffer // append the triangular part via a temporary buffer
for (Index j=0;j<actualPanelWidth;++j) for (Index j=0;j<actualPanelWidth;++j)
{ {
if (!(Mode&UnitDiag)) if (SetDiag)
triangularBuffer.coeffRef(j,j) = rhs(actual_j2+j,actual_j2+j); triangularBuffer.coeffRef(j,j) = rhs(actual_j2+j,actual_j2+j);
for (Index k=IsLower ? j+1 : 0; IsLower ? k<actualPanelWidth : k<j; ++k) for (Index k=IsLower ? j+1 : 0; IsLower ? k<actualPanelWidth : k<j; ++k)
triangularBuffer.coeffRef(k,j) = rhs(actual_j2+k,actual_j2+j); triangularBuffer.coeffRef(k,j) = rhs(actual_j2+k,actual_j2+j);

9
test/product_trmm.cpp
View File

@ -35,7 +35,7 @@ template<typename Scalar> void trmm(int size,int /*othersize*/)
DenseIndex cols = ei_random<DenseIndex>(1,size); DenseIndex cols = ei_random<DenseIndex>(1,size);
MatrixColMaj triV(rows,cols), triH(cols,rows), upTri(cols,rows), loTri(rows,cols), MatrixColMaj triV(rows,cols), triH(cols,rows), upTri(cols,rows), loTri(rows,cols),
unitUpTri(cols,rows), unitLoTri(rows,cols); unitUpTri(cols,rows), unitLoTri(rows,cols), strictlyUpTri(cols,rows), strictlyLoTri(rows,cols);
MatrixColMaj ge1(rows,cols), ge2(cols,rows), ge3; MatrixColMaj ge1(rows,cols), ge2(cols,rows), ge3;
MatrixRowMaj rge3; MatrixRowMaj rge3;
@ -48,6 +48,8 @@ template<typename Scalar> void trmm(int size,int /*othersize*/)
upTri = triH.template triangularView<Upper>(); upTri = triH.template triangularView<Upper>();
unitLoTri = triV.template triangularView<UnitLower>(); unitLoTri = triV.template triangularView<UnitLower>();
unitUpTri = triH.template triangularView<UnitUpper>(); unitUpTri = triH.template triangularView<UnitUpper>();
strictlyLoTri = triV.template triangularView<StrictlyLower>();
strictlyUpTri = triH.template triangularView<StrictlyUpper>();
ge1.setRandom(); ge1.setRandom();
ge2.setRandom(); ge2.setRandom();
@ -72,6 +74,11 @@ template<typename Scalar> void trmm(int size,int /*othersize*/)
VERIFY_IS_APPROX( rge3.noalias() = ge2 * triV.template triangularView<UnitLower>(), ge2 * unitLoTri); VERIFY_IS_APPROX( rge3.noalias() = ge2 * triV.template triangularView<UnitLower>(), ge2 * unitLoTri);
VERIFY_IS_APPROX( ge3 = ge2 * triV.template triangularView<UnitLower>(), ge2 * unitLoTri); VERIFY_IS_APPROX( ge3 = ge2 * triV.template triangularView<UnitLower>(), ge2 * unitLoTri);
VERIFY_IS_APPROX( ge3 = (s1*triV).adjoint().template triangularView<UnitUpper>() * ge2.adjoint(), ei_conj(s1) * unitLoTri.adjoint() * ge2.adjoint()); VERIFY_IS_APPROX( ge3 = (s1*triV).adjoint().template triangularView<UnitUpper>() * ge2.adjoint(), ei_conj(s1) * unitLoTri.adjoint() * ge2.adjoint());
VERIFY_IS_APPROX( ge3 = triV.template triangularView<StrictlyLower>() * ge2, strictlyLoTri * ge2);
VERIFY_IS_APPROX( rge3.noalias() = ge2 * triV.template triangularView<StrictlyLower>(), ge2 * strictlyLoTri);
VERIFY_IS_APPROX( ge3 = ge2 * triV.template triangularView<StrictlyLower>(), ge2 * strictlyLoTri);
VERIFY_IS_APPROX( ge3 = (s1*triV).adjoint().template triangularView<StrictlyUpper>() * ge2.adjoint(), ei_conj(s1) * strictlyLoTri.adjoint() * ge2.adjoint());
} }
void test_product_trmm() void test_product_trmm()