From f8f4c82e1feed186ac33c5e42e9b985ea1079e2f Mon Sep 17 00:00:00 2001 From: conrad Date: Thu, 20 Jul 2023 15:10:21 +1000 Subject: [PATCH] speedup by transpose based multiplication --- include/armadillo_bits/sp_auxlib_bones.hpp | 2 +- include/armadillo_bits/sp_auxlib_meat.hpp | 73 ++++++++++++++-------- 2 files changed, 47 insertions(+), 28 deletions(-) diff --git a/include/armadillo_bits/sp_auxlib_bones.hpp b/include/armadillo_bits/sp_auxlib_bones.hpp index 5a8e9ea8..4a06e74e 100644 --- a/include/armadillo_bits/sp_auxlib_bones.hpp +++ b/include/armadillo_bits/sp_auxlib_bones.hpp @@ -128,7 +128,7 @@ class sp_auxlib inline static void run_aupd_plain ( const uword n_eigvals, char* which, - const SpMat& X, const bool sym, + const SpMat& X, const SpMat& Xst, const bool sym, blas_int& n, eT& tol, blas_int& maxiter, podarray& resid, blas_int& ncv, podarray& v, blas_int& ldv, podarray& iparam, podarray& ipntr, diff --git a/include/armadillo_bits/sp_auxlib_meat.hpp b/include/armadillo_bits/sp_auxlib_meat.hpp index f8af9064..31e820ae 100644 --- a/include/armadillo_bits/sp_auxlib_meat.hpp +++ b/include/armadillo_bits/sp_auxlib_meat.hpp @@ -479,7 +479,9 @@ sp_auxlib::eigs_sym_arpack(Col& eigval, Mat& eigvec, const SpMat& X, } else { - run_aupd_plain(n_eigvals, which, X, true /* sym, not gen */, n, tol, maxiter, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); + const SpMat Xst = X.st(); + + run_aupd_plain(n_eigvals, which, X, Xst, true /* sym, not gen */, n, tol, maxiter, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); } if(info != 0) { return false; } @@ -863,7 +865,9 @@ sp_auxlib::eigs_gen_arpack(Col< std::complex >& eigval, Mat< std::complex } else { - run_aupd_plain(n_eigvals, which, X, false /* gen, not sym */, n, tol, maxiter, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); + const SpMat Xst = X.st(); + + run_aupd_plain(n_eigvals, which, X, Xst, false /* gen, not sym */, n, tol, maxiter, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); } if(info != 0) { return false; } @@ -1105,7 +1109,9 @@ sp_auxlib::eigs_gen(Col< std::complex >& eigval, Mat< std::complex >& eigv } else { - run_aupd_plain(n_eigvals, which, X, false /* gen, not sym */, n, tol, maxiter, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); + const SpMat< std::complex > Xst = X.st(); + + run_aupd_plain(n_eigvals, which, X, Xst, false /* gen, not sym */, n, tol, maxiter, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); } if(info != 0) { return false; } @@ -1893,7 +1899,7 @@ void sp_auxlib::run_aupd_plain ( const uword n_eigvals, char* which, - const SpMat& X, const bool sym, + const SpMat& X, const SpMat& Xst, const bool sym, blas_int& n, eT& tol, blas_int& maxiter, podarray& resid, blas_int& ncv, podarray& v, blas_int& ldv, podarray& iparam, podarray& ipntr, @@ -1978,30 +1984,46 @@ sp_auxlib::run_aupd_plain // where x is of length n and starts at workd(ipntr(0)), and y is of // length n and starts at workd(ipntr(1)). - // operator*(sp_mat, vec) doesn't properly put the result into the - // right place so we'll just reimplement it here for now... + // // OLD METHOD + // + // // operator*(sp_mat, vec) doesn't properly put the result into the + // // right place so we'll just reimplement it here for now... + // + // // Set the output to point at the right memory. We have to subtract + // // one from FORTRAN pointers... + // Col out(workd.memptr() + ipntr(1) - 1, n, false /* don't copy */); + // // Set the input to point at the right memory. + // Col in(workd.memptr() + ipntr(0) - 1, n, false /* don't copy */); + // + // out.zeros(); + // + // T* out_mem = out.memptr(); + // const T* in_mem = in.memptr(); + // + // typename SpMat::const_iterator X_it = X.begin(); + // + // const uword X_nnz = X.n_nonzero; + // + // for(uword count=0; count < X_nnz; ++count, ++X_it) + // { + // const eT X_it_val = (*X_it); + // const uword X_it_row = X_it.row(); + // const uword X_it_col = X_it.col(); + // + // out_mem[X_it_row] += X_it_val * in_mem[X_it_col]; + // } + // + // // No need to modify memory further since it was all done in-place. - // Set the output to point at the right memory. We have to subtract - // one from FORTRAN pointers... - Col out(workd.memptr() + ipntr(1) - 1, n, false /* don't copy */); - // Set the input to point at the right memory. - Col in(workd.memptr() + ipntr(0) - 1, n, false /* don't copy */); - out.zeros(); + // NEW METHOD + // + // both operator*(rowvec, sp_mat) and operator*(sp_mat, colvec) can now write to an existing object - typename SpMat::const_iterator X_it = X.begin(); - typename SpMat::const_iterator X_it_end = X.end(); + Row out(workd.memptr() + ipntr(1) - 1, n, false, true); + Row in(workd.memptr() + ipntr(0) - 1, n, false, true); - while(X_it != X_it_end) - { - const uword X_it_row = X_it.row(); - const uword X_it_col = X_it.col(); - - out[X_it_row] += (*X_it) * in[X_it_col]; - ++X_it; - } - - // No need to modify memory further since it was all done in-place. + out = in * Xst; break; } @@ -2234,9 +2256,6 @@ sp_auxlib::run_aupd_shiftinvert // where x is of length n and starts at workd(ipntr(0)), and y is of // length n and starts at workd(ipntr(1)). - // operator*(sp_mat, vec) doesn't properly put the result into the - // right place so we'll just reimplement it here for now... - // Set the output to point at the right memory. We have to subtract // one from FORTRAN pointers... Col out(workd.memptr() + ipntr(1) - 1, n, false /* don't copy */);