From a556b45abf18f1bd509daaf63b66b7d55e9fd291 Mon Sep 17 00:00:00 2001 From: jjesswan Date: Mon, 22 Apr 2024 21:56:26 -0400 Subject: add engine version --- .../Eigen/src/SparseLU/SparseLU_gemm_kernel.h | 280 +++++++++++++++++++++ 1 file changed, 280 insertions(+) create mode 100644 engine-ocean/Eigen/src/SparseLU/SparseLU_gemm_kernel.h (limited to 'engine-ocean/Eigen/src/SparseLU/SparseLU_gemm_kernel.h') diff --git a/engine-ocean/Eigen/src/SparseLU/SparseLU_gemm_kernel.h b/engine-ocean/Eigen/src/SparseLU/SparseLU_gemm_kernel.h new file mode 100644 index 0000000..e37c2fe --- /dev/null +++ b/engine-ocean/Eigen/src/SparseLU/SparseLU_gemm_kernel.h @@ -0,0 +1,280 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Gael Guennebaud +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_SPARSELU_GEMM_KERNEL_H +#define EIGEN_SPARSELU_GEMM_KERNEL_H + +namespace Eigen { + +namespace internal { + + +/** \internal + * A general matrix-matrix product kernel optimized for the SparseLU factorization. + * - A, B, and C must be column major + * - lda and ldc must be multiples of the respective packet size + * - C must have the same alignment as A + */ +template +EIGEN_DONT_INLINE +void sparselu_gemm(Index m, Index n, Index d, const Scalar* A, Index lda, const Scalar* B, Index ldb, Scalar* C, Index ldc) +{ + using namespace Eigen::internal; + + typedef typename packet_traits::type Packet; + enum { + NumberOfRegisters = EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS, + PacketSize = packet_traits::size, + PM = 8, // peeling in M + RN = 2, // register blocking + RK = NumberOfRegisters>=16 ? 4 : 2, // register blocking + BM = 4096/sizeof(Scalar), // number of rows of A-C per chunk + SM = PM*PacketSize // step along M + }; + Index d_end = (d/RK)*RK; // number of columns of A (rows of B) suitable for full register blocking + Index n_end = (n/RN)*RN; // number of columns of B-C suitable for processing RN columns at once + Index i0 = internal::first_default_aligned(A,m); + + eigen_internal_assert(((lda%PacketSize)==0) && ((ldc%PacketSize)==0) && (i0==internal::first_default_aligned(C,m))); + + // handle the non aligned rows of A and C without any optimization: + for(Index i=0; i(BM, m-ib); // actual number of rows + Index actual_b_end1 = (actual_b/SM)*SM; // actual number of rows suitable for peeling + Index actual_b_end2 = (actual_b/PacketSize)*PacketSize; // actual number of rows suitable for vectorization + + // Let's process two columns of B-C at once + for(Index j=0; j(Bc0[0]); } + { b10 = pset1(Bc0[1]); } + if(RK==4) { b20 = pset1(Bc0[2]); } + if(RK==4) { b30 = pset1(Bc0[3]); } + { b01 = pset1(Bc1[0]); } + { b11 = pset1(Bc1[1]); } + if(RK==4) { b21 = pset1(Bc1[2]); } + if(RK==4) { b31 = pset1(Bc1[3]); } + + Packet a0, a1, a2, a3, c0, c1, t0, t1; + + const Scalar* A0 = A+ib+(k+0)*lda; + const Scalar* A1 = A+ib+(k+1)*lda; + const Scalar* A2 = A+ib+(k+2)*lda; + const Scalar* A3 = A+ib+(k+3)*lda; + + Scalar* C0 = C+ib+(j+0)*ldc; + Scalar* C1 = C+ib+(j+1)*ldc; + + a0 = pload(A0); + a1 = pload(A1); + if(RK==4) + { + a2 = pload(A2); + a3 = pload(A3); + } + else + { + // workaround "may be used uninitialized in this function" warning + a2 = a3 = a0; + } + +#define KMADD(c, a, b, tmp) {tmp = b; tmp = pmul(a,tmp); c = padd(c,tmp);} +#define WORK(I) \ + c0 = pload(C0+i+(I)*PacketSize); \ + c1 = pload(C1+i+(I)*PacketSize); \ + KMADD(c0, a0, b00, t0) \ + KMADD(c1, a0, b01, t1) \ + a0 = pload(A0+i+(I+1)*PacketSize); \ + KMADD(c0, a1, b10, t0) \ + KMADD(c1, a1, b11, t1) \ + a1 = pload(A1+i+(I+1)*PacketSize); \ + if(RK==4){ KMADD(c0, a2, b20, t0) }\ + if(RK==4){ KMADD(c1, a2, b21, t1) }\ + if(RK==4){ a2 = pload(A2+i+(I+1)*PacketSize); }\ + if(RK==4){ KMADD(c0, a3, b30, t0) }\ + if(RK==4){ KMADD(c1, a3, b31, t1) }\ + if(RK==4){ a3 = pload(A3+i+(I+1)*PacketSize); }\ + pstore(C0+i+(I)*PacketSize, c0); \ + pstore(C1+i+(I)*PacketSize, c1) + + // process rows of A' - C' with aggressive vectorization and peeling + for(Index i=0; i0) + { + const Scalar* Bc0 = B+(n-1)*ldb; + + for(Index k=0; k(Bc0[0]); + b10 = pset1(Bc0[1]); + if(RK==4) b20 = pset1(Bc0[2]); + if(RK==4) b30 = pset1(Bc0[3]); + + Packet a0, a1, a2, a3, c0, t0/*, t1*/; + + const Scalar* A0 = A+ib+(k+0)*lda; + const Scalar* A1 = A+ib+(k+1)*lda; + const Scalar* A2 = A+ib+(k+2)*lda; + const Scalar* A3 = A+ib+(k+3)*lda; + + Scalar* C0 = C+ib+(n_end)*ldc; + + a0 = pload(A0); + a1 = pload(A1); + if(RK==4) + { + a2 = pload(A2); + a3 = pload(A3); + } + else + { + // workaround "may be used uninitialized in this function" warning + a2 = a3 = a0; + } + +#define WORK(I) \ + c0 = pload(C0+i+(I)*PacketSize); \ + KMADD(c0, a0, b00, t0) \ + a0 = pload(A0+i+(I+1)*PacketSize); \ + KMADD(c0, a1, b10, t0) \ + a1 = pload(A1+i+(I+1)*PacketSize); \ + if(RK==4){ KMADD(c0, a2, b20, t0) }\ + if(RK==4){ a2 = pload(A2+i+(I+1)*PacketSize); }\ + if(RK==4){ KMADD(c0, a3, b30, t0) }\ + if(RK==4){ a3 = pload(A3+i+(I+1)*PacketSize); }\ + pstore(C0+i+(I)*PacketSize, c0); + + // aggressive vectorization and peeling + for(Index i=0; i0) + { + for(Index j=0; j1 ? Aligned : 0 + }; + typedef Map, Alignment > MapVector; + typedef Map, Alignment > ConstMapVector; + if(rd==1) MapVector(C+j*ldc+ib,actual_b) += B[0+d_end+j*ldb] * ConstMapVector(A+(d_end+0)*lda+ib, actual_b); + + else if(rd==2) MapVector(C+j*ldc+ib,actual_b) += B[0+d_end+j*ldb] * ConstMapVector(A+(d_end+0)*lda+ib, actual_b) + + B[1+d_end+j*ldb] * ConstMapVector(A+(d_end+1)*lda+ib, actual_b); + + else MapVector(C+j*ldc+ib,actual_b) += B[0+d_end+j*ldb] * ConstMapVector(A+(d_end+0)*lda+ib, actual_b) + + B[1+d_end+j*ldb] * ConstMapVector(A+(d_end+1)*lda+ib, actual_b) + + B[2+d_end+j*ldb] * ConstMapVector(A+(d_end+2)*lda+ib, actual_b); + } + } + + } // blocking on the rows of A and C +} +#undef KMADD + +} // namespace internal + +} // namespace Eigen + +#endif // EIGEN_SPARSELU_GEMM_KERNEL_H -- cgit v1.2.3-70-g09d2