求助，矩阵相乘的优化？

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;
}

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++);  
                        }
                }
        }
}