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

improve trmm unit test and fix several bugs in trmm

Browse Source
This commit is contained in:
Gael Guennebaud 2010-06-15 23:38:21 +02:00
parent 2e792d1f42
commit 9726824f7c
2 changed files with 57 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -75,14 +75,14 @@ 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) | ((Mode&Upper) ? Lower : Upper),
(!LhsIsTriangular), (!LhsIsTriangular),
RhsStorageOrder==RowMajor ? ColMajor : RowMajor, RhsStorageOrder==RowMajor ? ColMajor : RowMajor,
ConjugateRhs, ConjugateRhs,
LhsStorageOrder==RowMajor ? ColMajor : RowMajor, LhsStorageOrder==RowMajor ? ColMajor : RowMajor,
ConjugateLhs, ConjugateLhs,
ColMajor> ColMajor>
::run(rows, cols, depth, rhs, rhsStride, lhs, lhsStride, res, resStride, alpha); ::run(cols, rows, depth, rhs, rhsStride, lhs, lhsStride, res, resStride, alpha);
} }
}; };
@ -138,6 +138,7 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,true,
Index actual_kc = std::min(IsLower ? k2 : depth-k2, kc); Index actual_kc = std::min(IsLower ? k2 : depth-k2, kc);
Index actual_k2 = IsLower ? k2-actual_kc : k2; Index actual_k2 = IsLower ? k2-actual_kc : k2;
// align blocks with the end of the triangular part for trapezoidal lhs
if((!IsLower)&&(k2<rows)&&(k2+actual_kc>rows)) if((!IsLower)&&(k2<rows)&&(k2+actual_kc>rows))
{ {
actual_kc = rows-k2; actual_kc = rows-k2;
@ -191,7 +192,7 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,true,
// the part below the diagonal => GEPP // the part below the diagonal => GEPP
{ {
Index start = IsLower ? k2 : 0; Index start = IsLower ? k2 : 0;
Index end = IsLower ? rows : actual_k2; Index end = IsLower ? rows : std::min(actual_k2,rows);
for(Index i2=start; i2<end; i2+=mc) for(Index i2=start; i2<end; i2+=mc)
{ {
const Index actual_mc = std::min(i2+mc,end)-i2; const Index actual_mc = std::min(i2+mc,end)-i2;
@ -258,14 +259,27 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,false,
IsLower ? k2<depth : k2>0; IsLower ? k2<depth : k2>0;
IsLower ? k2+=kc : k2-=kc) IsLower ? k2+=kc : k2-=kc)
{ {
const Index actual_kc = std::min(IsLower ? depth-k2 : k2, kc); Index actual_kc = std::min(IsLower ? depth-k2 : k2, kc);
Index actual_k2 = IsLower ? k2 : k2-actual_kc; Index actual_k2 = IsLower ? k2 : k2-actual_kc;
Index rs = IsLower ? actual_k2 : depth - k2;
Scalar* geb = blockB+actual_kc*actual_kc; // align blocks with the end of the triangular part for trapezoidal rhs
if(IsLower && (k2<cols) && (actual_k2+actual_kc>cols))
{
actual_kc = cols-k2;
k2 = actual_k2 + actual_kc - kc;
}
// remaining size
Index rs = IsLower ? std::min(cols,actual_k2) : cols - k2;
// size of the triangular part
Index ts = (IsLower && actual_k2>=cols) ? 0 : actual_kc;
Scalar* geb = blockB+ts*ts;
pack_rhs(geb, &rhs(actual_k2,IsLower ? 0 : k2), rhsStride, alpha, actual_kc, rs); pack_rhs(geb, &rhs(actual_k2,IsLower ? 0 : k2), rhsStride, alpha, actual_kc, rs);
// pack the triangular part of the rhs padding the unrolled blocks with zeros // pack the triangular part of the rhs padding the unrolled blocks with zeros
if(ts>0)
{ {
for (Index j2=0; j2<actual_kc; j2+=SmallPanelWidth) for (Index j2=0; j2<actual_kc; j2+=SmallPanelWidth)
{ {
@ -301,6 +315,7 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,false,
pack_lhs(blockA, &lhs(i2, actual_k2), lhsStride, actual_kc, actual_mc); pack_lhs(blockA, &lhs(i2, actual_k2), lhsStride, actual_kc, actual_mc);
// triangular kernel // triangular kernel
if(ts>0)
{ {
for (Index j2=0; j2<actual_kc; j2+=SmallPanelWidth) for (Index j2=0; j2<actual_kc; j2+=SmallPanelWidth)
{ {

66
test/product_trmm.cpp
View File

@ -24,48 +24,54 @@
#include "main.h" #include "main.h"
template<typename Scalar> void trmm(int size,int othersize) template<typename Scalar> void trmm(int size,int /*othersize*/)
{ {
typedef typename NumTraits<Scalar>::Real RealScalar; typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Matrix<Scalar,Dynamic,Dynamic,ColMajor> MatrixType; typedef Matrix<Scalar,Dynamic,Dynamic,ColMajor> MatrixColMaj;
typedef Matrix<Scalar,Dynamic,Dynamic,RowMajor> MatrixRowMaj;
MatrixType tri(size,size), upTri(size,size), loTri(size,size), int rows = size;
unitUpTri(size,size), unitLoTri(size,size); int cols = ei_random<int>(1,size);
MatrixType ge1(size,othersize), ge2(10,size), ge3;
Matrix<Scalar,Dynamic,Dynamic,RowMajor> rge3; MatrixColMaj triV(rows,cols), triH(cols,rows), upTri(cols,rows), loTri(rows,cols),
unitUpTri(cols,rows), unitLoTri(rows,cols);
MatrixColMaj ge1(rows,cols), ge2(cols,rows), ge3;
MatrixRowMaj rge3;
Scalar s1 = ei_random<Scalar>(), Scalar s1 = ei_random<Scalar>(),
s2 = ei_random<Scalar>(); s2 = ei_random<Scalar>();
tri.setRandom(); triV.setRandom();
loTri = tri.template triangularView<Lower>(); triH.setRandom();
upTri = tri.template triangularView<Upper>(); loTri = triV.template triangularView<Lower>();
unitLoTri = tri.template triangularView<UnitLower>(); upTri = triH.template triangularView<Upper>();
unitUpTri = tri.template triangularView<UnitUpper>(); unitLoTri = triV.template triangularView<UnitLower>();
unitUpTri = triH.template triangularView<UnitUpper>();
ge1.setRandom(); ge1.setRandom();
ge2.setRandom(); ge2.setRandom();
VERIFY_IS_APPROX( ge3 = tri.template triangularView<Lower>() * ge1, loTri * ge1); VERIFY_IS_APPROX( ge3 = triV.template triangularView<Lower>() * ge2, loTri * ge2);
VERIFY_IS_APPROX(rge3 = tri.template triangularView<Lower>() * ge1, loTri * ge1); VERIFY_IS_APPROX( ge3 = ge2 * triV.template triangularView<Lower>(), ge2 * loTri);
VERIFY_IS_APPROX( ge3 = tri.template triangularView<Upper>() * ge1, upTri * ge1); VERIFY_IS_APPROX( ge3 = triH.template triangularView<Upper>() * ge1, upTri * ge1);
VERIFY_IS_APPROX(rge3 = tri.template triangularView<Upper>() * ge1, upTri * ge1); VERIFY_IS_APPROX( ge3 = ge1 * triH.template triangularView<Upper>(), ge1 * upTri);
VERIFY_IS_APPROX( ge3 = (s1*tri.adjoint()).template triangularView<Upper>() * (s2*ge1), s1*loTri.adjoint() * (s2*ge1)); VERIFY_IS_APPROX( ge3 = (s1*triV.adjoint()).template triangularView<Upper>() * (s2*ge1), s1*loTri.adjoint() * (s2*ge1));
VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView<Upper>() * ge1, loTri.adjoint() * ge1); VERIFY_IS_APPROX( ge3 = ge1 * triV.adjoint().template triangularView<Upper>(), ge1 * loTri.adjoint());
VERIFY_IS_APPROX( ge3 = tri.adjoint().template triangularView<Lower>() * ge1, upTri.adjoint() * ge1); VERIFY_IS_APPROX( ge3 = triH.adjoint().template triangularView<Lower>() * ge2, upTri.adjoint() * ge2);
VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView<Lower>() * ge1, upTri.adjoint() * ge1); VERIFY_IS_APPROX( ge3 = ge2 * triH.adjoint().template triangularView<Lower>(), ge2 * upTri.adjoint());
VERIFY_IS_APPROX( ge3 = tri.template triangularView<Lower>() * ge2.adjoint(), loTri * ge2.adjoint()); VERIFY_IS_APPROX( ge3 = triV.template triangularView<Lower>() * ge1.adjoint(), loTri * ge1.adjoint());
VERIFY_IS_APPROX(rge3 = tri.template triangularView<Lower>() * ge2.adjoint(), loTri * ge2.adjoint()); VERIFY_IS_APPROX( ge3 = ge1.adjoint() * triV.template triangularView<Lower>(), ge1.adjoint() * loTri);
VERIFY_IS_APPROX( ge3 = tri.template triangularView<Upper>() * ge2.adjoint(), upTri * ge2.adjoint()); VERIFY_IS_APPROX( ge3 = triH.template triangularView<Upper>() * ge2.adjoint(), upTri * ge2.adjoint());
VERIFY_IS_APPROX(rge3 = tri.template triangularView<Upper>() * ge2.adjoint(), upTri * ge2.adjoint()); VERIFY_IS_APPROX(rge3.noalias() = triH.template triangularView<Upper>() * ge2.adjoint(), upTri * ge2.adjoint());
VERIFY_IS_APPROX( ge3 = (s1*tri).adjoint().template triangularView<Upper>() * ge2.adjoint(), ei_conj(s1) * loTri.adjoint() * ge2.adjoint()); VERIFY_IS_APPROX( ge3 = (s1*triV).adjoint().template triangularView<Upper>() * ge2.adjoint(), ei_conj(s1) * loTri.adjoint() * ge2.adjoint());
VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView<Upper>() * ge2.adjoint(), loTri.adjoint() * ge2.adjoint()); VERIFY_IS_APPROX(rge3.noalias() = triV.adjoint().template triangularView<Upper>() * ge2.adjoint(), loTri.adjoint() * ge2.adjoint());
VERIFY_IS_APPROX( ge3 = tri.adjoint().template triangularView<Lower>() * ge2.adjoint(), upTri.adjoint() * ge2.adjoint()); VERIFY_IS_APPROX( ge3 = triH.adjoint().template triangularView<Lower>() * ge1.adjoint(), upTri.adjoint() * ge1.adjoint());
VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView<Lower>() * ge2.adjoint(), upTri.adjoint() * ge2.adjoint()); VERIFY_IS_APPROX(rge3.noalias() = triH.adjoint().template triangularView<Lower>() * ge1.adjoint(), upTri.adjoint() * ge1.adjoint());
VERIFY_IS_APPROX( ge3 = tri.template triangularView<UnitLower>() * ge1, unitLoTri * ge1); VERIFY_IS_APPROX( ge3 = triV.template triangularView<UnitLower>() * ge2, unitLoTri * ge2);
VERIFY_IS_APPROX(rge3 = tri.template triangularView<UnitLower>() * ge1, unitLoTri * ge1); VERIFY_IS_APPROX( rge3.noalias() = ge2 * triV.template triangularView<UnitLower>(), ge2 * unitLoTri);
VERIFY_IS_APPROX( ge3 = (s1*tri).adjoint().template triangularView<UnitUpper>() * ge2.adjoint(), ei_conj(s1) * unitLoTri.adjoint() * ge2.adjoint()); 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());
} }
void test_product_trmm() void test_product_trmm()