For sure, there has been many highly optimized library on this. But I am working on a matrix-free context since the problem size does not allow explicit storage of sparse matrix elements. I'd love have a general knowledge on how people handle this sort of problem, especially on how to alleviate the memory bottleneck.


As suggested by other commentators, I shall further elaborate on my situation.

My current implementation is below, which is just the basic one showing up in various textbooks.

void SpMVMtply(vector<det>& basis, ColVec& x, ColVec& b)

         "basis" is the orthonormal basis that spans the physical space 
         "x" ,"b" are just those in Ax = b
    int len = x.size(); 

    #pragma omp parallel
        #pragma omp for nowait 
        for (int i=0; i<len; i++)
            vector<det> pool_det; 
            vector<cdouble> pool_amp; 

               this function generates the non-zero off diagonal elements 
               with respect to the given basis vector. pool_det is the 
               collection of generated basis vectors, and pool_map collection 
               of the amplitude. 

            OffDiagGen(basis[i], pool_det, pool_amp);

            // this function generates the diagonal element
            cdouble sum = x(i)* DiagCal(basis[i]);

                Below, basis_pos is an unordered_map with basis vector as key 
                and that basis vector's position in "basis" array as value. 

            for (int j=0; j< pool_det.size(); j++) 
                sum += x(basis_pos[pool_det[j]]) * conj(pool_amp[j]); 

            b(i) = sum ;  

The size of vector could reach 1 billion, so the explicit storage of sparse matrix is impossible. A rough profiling shows that the non-sequential accessing of memory in x(basis_pos[pool_det[j]]) is the most time-consuming part. On a Xeon Gold 6242 with 12 cores in use, it takes 6.2 seconds to such part for a vector with 3 million elements, while it takes only 1.5 seconds when that part is replaced by a constant. That's why I call it a memory bottleneck. My main concern is that if such bottleneck is inevitable, since even when the sparse matrix is explicitly stored, such non-sequential memory access also shows up. Based on my knowledge I haven't figured out a way to avoid it.

  • 1
    $\begingroup$ It is a bit unclear what aspects you are interested in. By memory bottleneck, do you imply memory access or general "not enough memory" aspect? Any details on your matrix-free context and why do you still have memory bottleneck in the matrix-free solver? What paradigm of parallel computing are you interested in (distributed, shared, hybrid, GPU, FPGA)? Have you looked into any libraries that implement this functionality? Currently the question is extremely broad, I want to help you narrow it down. $\endgroup$ – Anton Menshov Apr 1 at 20:01
  • $\begingroup$ You want to read the papers by Martin Kronbichler! $\endgroup$ – Wolfgang Bangerth Apr 2 at 1:13
  • $\begingroup$ So you have a sparse matrix but not enough storage to hold the sparse matrix- you can only compute elements of the sparse matrix as needed, right? Can you iterate over all the nonzero entries in parallel in a reasonable way? Do you have enough storage to hold the vector input and the vector output? $\endgroup$ – Brian Borchers Apr 2 at 3:55
  • $\begingroup$ I have updated my question in order to address some ambiguities. $\endgroup$ – Izzy Vang Apr 2 at 22:02
  • $\begingroup$ @BrianBorchers Yes, only input and output vectors are stored explicitly, and there is a way to iterate over all non-zero elements. $\endgroup$ – Izzy Vang Apr 2 at 22:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.