求助，矩阵相乘的优化？

hangover

大家有什么好的算法介绍一下，昨天看了一段code，区别不大阿，各位指点一下，谢谢

Matrix operator*(const Matrix& B)
{
        int B_M = B.getRow();
        int B_N = B.getColumn();
        if( N!= B_M)
        {
                throw "Dimension for multiplication do NOT match !" ;
        }

        Matrix C(M,N);
        /*-------------------------------------------------------------
        //  Standard implemetion of matrix multiplication formula

        for(int i=0; i<M; i++)
        for(int j=0;j<B_N; j++)
        {     
        C.data[i*N+j]= 0;
        for(int k=0; k<N;k++)
        C.data[i*N+j]+= data[i*N+k]* B.data[k*N+j];
        }
        //----------------------------------------------------------- */

        //-------------------------------------------------------------
        // An optimized way to speed up the multiplication calculation

        for(int i=0, e=0; i<M; i++)
                for(int j=0;j<B_N; j++, e++)
                {   
                        T val = data[i*N] * B.data[j];
                        for(int k=1, c = i * N+1, d = j + B_N; k<N;k++, c++, d += B_N)
                                val += data[c]* B.data[d];

                        C.data[e] = val;
                }               
                //---------------------------------------------------
                return C;
}

longshort

测试一下它们的处理器消耗，就知道了。

hangover

void dgemm(const Matrix &A, const Matrix &B, Matrix &C)
竟然比
Matrix operator % (const Matrix &a, const Matrix &b)
慢了好多。
operator overloading用的是最基本的算法,return 一个matrix，
而dgemm是pure pointer access但是多了一个loop.
怎么回事啊？好多人都说pure access要快一些，高手指点一下啊。



Matrix operator % (const Matrix &a, const Matrix &b)
{
#ifdef DEBUG__
        if (a.ncol_ != b.nrow_)
                errorReport("(Matrix operator %) : matrices do not conform in multiplication");
#endif

        Matrix c(a.nrow_, b.ncol_);
        double *ptr_a = a.data_, *ptr_b = b.data_, *ptr_c = c.data_;
        double sum;
        for (int i = 0; i < a.nrow_; i++)
                for (int j = 0; j < b.ncol_; j++)
                {
                        sum = 0.0;
                        for (int k = 0; k < a.ncol_; k++)
                                sum += ptr_a[i*a.ncol_+k]*ptr_b[k*b.ncol_+j];
                        ptr_c[i*c.ncol_+j]=sum;
                }

                return c;
}



void dgemm(const Matrix &A, const Matrix &B, Matrix &C)
// A[m][n] * B[n][l] = C[m][l]
{  
        int i, k, j;
        const double *ptr_a = A.data();
        int m = A.nrows(), n = A.ncols(), l = B.ncols();

        for ( i = 0; i < m; i++ )
        {  
                double *ptr_c = C.data() + i*l, Aik = *ptr_a++;
                const double *ptr_b = B.data();

                for ( j = 0; j < l; j++ )
                        *ptr_c++ = Aik * (*ptr_b++);

                for ( k = 1; k < n; k++ )
                {  
                        ptr_b  = B.data() + k*l;
                        ptr_c  = C.data() + i*l;
                        Aik = *ptr_a++;

                        for ( j = 0; j < l; j++ )
                        {  
                                *ptr_c++ += Aik * (*ptr_b++);  
                        }
                }
        }
}

sl_stanley

这个只有看编译后的汇编代码了,不过用数组下标访问和用指针,在现在的编译器里也没啥区别

tyc611

1. 两个矩阵相乘请参考斯特拉森(V.Strassen)乘法；
2. 多个矩阵连乘请参考矩阵连乘的动态规划算法（就是加括号确定一个最优计算顺序，不处理和处理的效率往往相差非常大！）
这些都可以在讲解算法的书中找到，另外可以参考：
线性代数与多项式的快速算法，游兆主编，上海科学技术出版社，1980（这本书我没看过，是算法书上在讲解矩阵乘法时推荐的，据说对矩阵相乘做了深入讨论）

ArXoR

要比较彻底一点的优化就用FFT吧