phynexis-deps/armadillo-code
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

speedup by transpose based multiplication

Browse Source
This commit is contained in:
conrad
2023-07-20 15:10:21 +10:00
parent 20db1e4eea
commit f8f4c82e1f
2 changed files with 47 additions and 28 deletions
Download Patch File Download Diff File Expand all files Collapse all files
include/armadillo_bits/sp_auxlib_bones.hpp
+1 -1
View File
@@ -128,7 +128,7 @@ class sp_auxlib
inline static void run_aupd_plain
(
const uword n_eigvals, char* which,
const SpMat<T>& X, const bool sym,
const SpMat<T>& X, const SpMat<T>& Xst, const bool sym,
blas_int& n, eT& tol, blas_int& maxiter,
podarray<T>& resid, blas_int& ncv, podarray<T>& v, blas_int& ldv,
podarray<blas_int>& iparam, podarray<blas_int>& ipntr,
include/armadillo_bits/sp_auxlib_meat.hpp
+46 -27
View File
@@ -479,7 +479,9 @@ sp_auxlib::eigs_sym_arpack(Col<eT>& eigval, Mat<eT>& eigvec, const SpMat<eT>& 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<eT> 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<T> >& eigval, Mat< std::complex<T>
}
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<T> 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<T> >& eigval, Mat< std::complex<T> >& 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<T> > 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<T>& X, const bool sym,
const SpMat<T>& X, const SpMat<T>& Xst, const bool sym,
blas_int& n, eT& tol, blas_int& maxiter,
podarray<T>& resid, blas_int& ncv, podarray<T>& v, blas_int& ldv,
podarray<blas_int>& iparam, podarray<blas_int>& 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<T> out(workd.memptr() + ipntr(1) - 1, n, false /* don't copy */);
// // Set the input to point at the right memory.
// Col<T> 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<T>::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<T> out(workd.memptr() + ipntr(1) - 1, n, false /* don't copy */);
// Set the input to point at the right memory.
Col<T> 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<T>::const_iterator X_it = X.begin();
typename SpMat<T>::const_iterator X_it_end = X.end();
Row<T> out(workd.memptr() + ipntr(1) - 1, n, false, true);
Row<T> 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<T> out(workd.memptr() + ipntr(1) - 1, n, false /* don't copy */);