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

2 files changed, 23 insertions, 8 deletions

diff --git a/Eigen/src/Core/products/TriangularMatrixMatrix.h b/Eigen/src/Core/products/TriangularMatrixMatrix.h
index be9362958..98305f993 100644
--- a/Eigen/src/Core/products/TriangularMatrixMatrix.h
+++ b/Eigen/src/Core/products/TriangularMatrixMatrix.h
@@ -75,7 +75,7 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,LhsIsTriangular,
     Scalar alpha)
   {
     ei_product_triangular_matrix_matrix<Scalar, Index,
-      (Mode&UnitDiag) | ((Mode&Upper) ? Lower : Upper),
+      (Mode&(UnitDiag|ZeroDiag)) | ((Mode&Upper) ? Lower : Upper),
       (!LhsIsTriangular),
       RhsStorageOrder==RowMajor ? ColMajor : RowMajor,
       ConjugateRhs,
@@ -111,7 +111,8 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,true,
     typedef ei_product_blocking_traits<Scalar> Blocking;
     enum {
       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
@@ -127,7 +128,10 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,true,
 
     Matrix<Scalar,SmallPanelWidth,SmallPanelWidth,LhsStorageOrder> triangularBuffer;
     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_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
           for (Index k=0;k<actualPanelWidth;++k)
           {
-            if (!(Mode&UnitDiag))
+            if (SetDiag)
               triangularBuffer.coeffRef(k,k) = lhs(startBlock+k,startBlock+k);
             for (Index i=IsLower ? k+1 : 0; IsLower ? i<actualPanelWidth : i<k; ++i)
               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;
     enum {
       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
@@ -252,7 +257,10 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,false,
 
     Matrix<Scalar,SmallPanelWidth,SmallPanelWidth,RhsStorageOrder> triangularBuffer;
     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_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
           for (Index j=0;j<actualPanelWidth;++j)
           {
-            if (!(Mode&UnitDiag))
+            if (SetDiag)
               triangularBuffer.coeffRef(j,j) = rhs(actual_j2+j,actual_j2+j);
             for (Index k=IsLower ? j+1 : 0; IsLower ? k<actualPanelWidth : k<j; ++k)
               triangularBuffer.coeffRef(k,j) = rhs(actual_j2+k,actual_j2+j);
diff --git a/test/product_trmm.cpp b/test/product_trmm.cpp
index e20b408c4..e4790fb66 100644
--- a/test/product_trmm.cpp
+++ b/test/product_trmm.cpp
@@ -35,7 +35,7 @@ template<typename Scalar> void trmm(int size,int /*othersize*/)
   DenseIndex cols = ei_random<DenseIndex>(1,size);
 
   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;
   MatrixRowMaj  rge3;
 
@@ -48,6 +48,8 @@ template<typename Scalar> void trmm(int size,int /*othersize*/)
   upTri = triH.template triangularView<Upper>();
   unitLoTri = triV.template triangularView<UnitLower>();
   unitUpTri = triH.template triangularView<UnitUpper>();
+  strictlyLoTri = triV.template triangularView<StrictlyLower>();
+  strictlyUpTri = triH.template triangularView<StrictlyUpper>();
   ge1.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( 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 = 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()