Hello Xiss,

maybe the examples on sparse fill can help you with that. Take a look at: http://www.guwi17.de/ublas/matrix_sparse_usage.html. It is probable that using a generalized vector of coordinate matrices to intitally assemble the matrix (by basically changing the type you have) and then use that to fill (by using push_back) a compressed matrix, will make it a lot faster. If you don't need the structure of the compressed matrix you can neglect the last step.

If you are sort in memory, another idea that might work (although it is algorithmically harder and I haven't really tested  it) is to define a vector of matrix blocks (each matrix being maybe 2*band x 2*band), that you use them to assemble the stiffness matrix and then progressively push them back in the compressed matrix like,
[B1 C1  0  0 ]
[A2 B2 C2 0 ]
[0 A3 B3 C3]
[0  0  A4 B4]
neglecting the zero entries. "Bs" will be diagonals, "As "will be at least upper triangular and "Cs" will be at least lower triangular (the least means that their diagonal elements maybe zero). This would need some custom algorithms to navigate in the blocks though.

Best
Nasos

---

Date: Mon, 19 Apr 2010 00:14:36 -0300
From: xissburg@gmail.com
To: ublas@lists.boost.org
Subject: [ublas] Element-wise operations are really slow

In my algorithms I have to read/write from/to each individual element of  a matrix, and this is making my application really slow. More specifically, I'm assembling a stiffness matrix in Finite Element Method. The code is like this:

    for(int i=0; i<m_tetrahedrons.size(); ++i)
    {
        btTetrahedron* t = m_tetrahedrons[i];
        t->computeCorotatedStiffness();

        for(unsigned int j=0; j<4; ++j)
            for(unsigned int k=0; k<4; ++k)
            {
                unsigned int jj = t->getNodeIndex(j);
                unsigned int kk = t->getNodeIndex(k);

                for(unsigned int r=0; r<3; ++r)
                    for(unsigned int s=0; s<3; ++s)
                    {
                        m_RKR_1(3*jj+r, 3*kk+s) += t->getCorotatedStiffness0(3*j+r, 3*k+s);
                        m_RK(3*jj+r, 3*kk+s) += t->getCorotatedStiffness1(3*j+r, 3*k+s);
                    }
            }
    }

Where m_RKR_1 and m_RK are both compressed_matrix<float>, t->getCorotatedStiffness0/1 just returns the (i,j) element of a 12x12 compressed_matrix<float>. If I don't compute the co-rotated matrices the simulation still works but incorrectly (linear strain only), but very fast (basically, by commenting the element wise operations in that code). Whenever I turn the co-rotational stuff on, it gets damn slow, and those element wise operations are guilty.

What am I doing wrong? Is there any faster technique to do that? Well, there must be one...


Thanks in advance.

---

The New Busy think 9 to 5 is a cute idea. Combine multiple calendars with Hotmail. Get busy.