From 35927e78c21b11113889f93325b980fc13349af9 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 24 Jul 2009 16:21:52 +0200 Subject: [PATCH] add WIP trsm --- Eigen/src/Core/products/GeneralMatrixMatrix.h | 20 +- .../Core/products/TriangularSolverMatrix.h | 203 ++++++++++++++++++ 2 files changed, 214 insertions(+), 9 deletions(-) create mode 100644 Eigen/src/Core/products/TriangularSolverMatrix.h diff --git a/Eigen/src/Core/products/GeneralMatrixMatrix.h b/Eigen/src/Core/products/GeneralMatrixMatrix.h index 776b7c5d6..a1be8ea50 100644 --- a/Eigen/src/Core/products/GeneralMatrixMatrix.h +++ b/Eigen/src/Core/products/GeneralMatrixMatrix.h @@ -110,10 +110,12 @@ static void run(int rows, int cols, int depth, template struct ei_gebp_kernel { - void operator()(Scalar* res, int resStride, const Scalar* blockA, const Scalar* blockB, int rows, int depth, int cols) + void operator()(Scalar* res, int resStride, const Scalar* blockA, const Scalar* blockB, int rows, int depth, int cols, int strideA=-1, int strideB=-1, int offsetA=0, int offsetB=0) { typedef typename ei_packet_traits::type PacketType; enum { PacketSize = ei_packet_traits::size }; + if(strideA==-1) strideA = depth; + if(strideB==-1) strideB = depth; Conj cj; int packet_cols = (cols/nr) * nr; const int peeled_mc = (rows/mr)*mr; @@ -123,7 +125,7 @@ struct ei_gebp_kernel // loops on each register blocking of lhs/res for(int i=0; i +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_TRIANGULAR_SOLVER_MATRIX_H +#define EIGEN_TRIANGULAR_SOLVER_MATRIX_H + +/* Optimized triangular solver with multiple right hand side (_TRSM) + */ +template +struct ei_triangular_solve_matrix// +{ + + static EIGEN_DONT_INLINE void run( + int size, int cols, + const Scalar* _lhs, int lhsStride, + Scalar* _rhs, int rhsStride) + { + Map > lhs(_lhs, size, size); + Map > rhs(_rhs, size, cols); + //ei_const_blas_data_mapper lhs(_lhs,lhsStride); + //ei_const_blas_data_mapper rhs(_rhs,rhsStride); + + typedef ei_product_blocking_traits Blocking; + enum { + SmallPanelWidth = EIGEN_ENUM_MAX(Blocking::mr,Blocking::nr), + IsLowerTriangular = (Mode&LowerTriangular) == LowerTriangular + }; + + int kc = 8;//std::min(Blocking::Max_kc,size); // cache block size along the K direction + int mc = 8;//std::min(Blocking::Max_mc,size); // cache block size along the M direction + + Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc); + Scalar* blockB = ei_aligned_stack_new(Scalar, kc*cols*Blocking::PacketSize); + + ei_gebp_kernel > gebp_kernel; + + for(int k2=0; k2 general block copy + // - R1 = L1^-1 B => tricky part + // - update B from the new R1 => actually this has to performed continuously during the above step + // - R2 = L2 * B => GEPP + + // B = R1 + ei_gemm_pack_rhs() + (blockB, &rhs(k2,0), rhsStride, -1, actual_kc, cols); + + Map(blockB,Blocking::PacketSize*Blocking::nr*actual_kc, cols/Blocking::nr+(cols%Blocking::nr)).setZero(); + + // The tricky part: R1 = L1^-1 B while updating B from R1 + // The idea is to split L1 into multiple small vertical panels. + // Each panel can be split into a small triangular part A1 which is processed without optimization, + // and the remaining small part A2 which is processed using gebp with appropriate block strides + { + // pack L1 +// ei_gemm_pack_lhs() +// (blockA, &lhs(k2, k2), lhsStride, actual_kc, actual_kc); + + // for each small vertical panels of lhs + for (int k1=0; k1(SmallPanelWidth,actual_kc-k1); + // tr solve + for (int k=0; k0) + { + rhs.block(i+1,0,rs,cols) -= + lhs.col(i).segment(IsLowerTriangular ? i+1 : i-rs, rs) * rhs.row(i); + } + } + // update the respective row of B from rhs + { + const Scalar* lr = _rhs+k2+k1; + int packet_cols = (cols/Blocking::nr) * Blocking::nr; + int count = 0; + for(int j2=0; j2(blockB,Blocking::PacketSize*Blocking::nr*actual_kc, cols/Blocking::nr+(cols%Blocking::nr)) << "\n\n"; + +// MatrixXf aux(Blocking::PacketSize*Blocking::nr*actual_kc, cols/Blocking::nr+(cols%Blocking::nr)); +// aux.setZero(); + +// ei_gemm_pack_rhs() +// (aux.data(), &rhs(k2,0), rhsStride, -1, actual_kc, cols); + +// std::cerr << Map(blockB,Blocking::PacketSize*Blocking::nr*actual_kc, cols/Blocking::nr+(cols%Blocking::nr)) - aux << "\n\n"; + + + // gebp + int i = k1+actualPanelWidth; + int rs = actual_kc-i; + +// ei_gemm_pack_rhs() +// (blockB, &rhs(k1,0), rhsStride, -1, actualPanelWidth, cols); + + ei_gemm_pack_lhs() + (blockA, &lhs(k2+i, k2+k1), lhsStride, actualPanelWidth, rs); + + if (rs>0) + rhs.block(i,0,actual_kc-i,cols) -= lhs.block(i,k1,rs,actualPanelWidth) * rhs.block(k1,0,actualPanelWidth,cols); + +// gebp_kernel(_rhs+i+k2, rhsStride, +// blockA/*+actual_kc*i+k1*rs*/, blockB/*+k1*Blocking::PacketSize*Blocking::nr*/, rs, actualPanelWidth, cols, actualPanelWidth/*actual_kc*/, actual_kc, 0, k1*Blocking::PacketSize); + +// gebp_kernel(_rhs+i, rhsStride, +// blockA+actual_kc*i+k1*rs, blockB+k1*Blocking::PacketSize*Blocking::nr, rs, actualPanelWidth, cols, actual_kc, actual_kc); + +// gebp_kernel(_rhs+k2+i, rhsStride, +// blockA+actual_kc*i+k1, blockB+k1*Blocking::PacketSize, actual_kc-i, actualPanelWidth, cols, actual_kc, actual_kc); + } + } + + // - R2 = A2 * B => GEPP + for(int i2=k2+kc; i2() + (blockA, &lhs(k2, i2), lhsStride, actual_kc, actual_mc); + + std::cerr << i2 << " sur " << actual_mc << " -= " << i2 << "x" << k2 << "+" << actual_mc<<"," <